# Path-integral simulations with fermionic and bosonic reservoirs:
Transport and dissipation in molecular electronic junctions

Lena Simine  Chemical Physics Theory Group, Department of Chemistry,
University of Toronto, 80 Saint George St. Toronto, Ontario, Canada M5S 3H6
Dvira Segal  Chemical Physics Theory Group, Department of Chemistry,
University of Toronto, 80 Saint George St. Toronto, Ontario, Canada M5S 3H6

###### Abstract

We expand iterative numerically-exact influence functional path-integral tools
and present a method capable of following the nonequilibrium time evolution of
subsystems coupled to multiple bosonic and fermionic reservoirs
simultaneously. Using this method, we study the real-time dynamics of charge
transfer and vibrational mode excitation in an electron conducting molecular
junction. We focus on nonequilibrium vibrational effects, particularly, the
development of vibrational instability in a current-rectifying junction. Our
simulations are performed by assuming large molecular vibrational
anharmonicity (or low temperature). This allows us to truncate the molecular
vibrational mode to include only a two-state system. Exact numerical results
are compared to perturbative Master equation calculations demonstrating an
excellent agreement in the weak electron-phonon coupling regime. Significant
deviations take place only at strong coupling. Our simulations allow us to
quantify the contribution of different transport mechanisms, coherent dynamics
and inelastic transport, in the overall charge current. This is done by
studying two model variants: The first admits inelastic electron transmission
only, while the second one allows for both coherent and incoherent pathways.

##  I Introduction

Following the quantum dynamics of an open-dissipative many-body system with
multiple bosonic and fermionic reservoirs in a nonequilibrium state, beyond
the linear response regime, is a significant theoretical and computational
challenge. In the realm of molecular conducting junctions, we should describe
the out-of-equilibrium dynamics of the molecular unit while handling both
electrons and molecular vibrations, accounting for many-body effects such as
electron-electron, phonon-phonon and electron-phonon interactions. Given this
complexity, studies in this field are mostly focused on steady-state
properties, using e.g., scattering theory Bonca and Trugmann (1995); Ness  _et
al._ (2001); Cízek  _et al._ (2004), while ignoring vibrational nonequilibrium
effects. Perturbative treatments (in either the molecule-leads coupling
parameter or the electron-phonon interaction energy) are commonly used,
including the nonequilibrium Green’s function technique Galperin  _et al._
(2007a); Mitra  _et al._ (2004); Galperin  _et al._ (2007b, 2009); Fransson
and Galperin (2010) and Master equation approaches Mitra  _et al._ (2004);
Segal and Nitzan (2002); Leijnse and Wegewijs (2008); Hartle and Thoss (2011a,
b); Volkovich  _et al._ (2010). For following the real-time dynamics of such
systems, involved methods have been recently developed, e.g., semiclassical
approaches Swenson  _et al._ (2011, 2012).

Figure 1: Left Panel: Generic setup considered in this work, including a
subsystem ($S$) coupled to multiple fermionic ($F$) and bosonic ($B$)
reservoirs. Right panel: Molecular electronic realization with two metals, $L$
and $R$, connected by two electronic levels, $D$ and $A$. Electronic
transitions in this junction are coupled to excitation/de-excitation processes
of a particular, anharmonic, vibrational mode that plays the role of the
“subsystem”. This mode may dissipate its excess energy to a secondary phonon
bath $B$.

In this work, we extend numerically-exact path-integral methods, and follow
the dynamics of a subsystem coupled to multiple out-of-equilibrium bosonic and
fermionic reservoirs. The technique is then applied on a molecular junction
realization, with the motivation to address basic problems in the field of
molecular electronics. Particularly, in this work we consider the dynamics and
steady-state properties of a conducting molecular junction acting as a charge
rectifier. A scheme of the generic setup and a particular molecular junction
realization are depicted in Fig. 1.

The time evolution scheme developed in this paper treats both bosonic and
fermionic reservoirs. This is achieved by combining two related iterative
path-integral methods: (i) The quasi-adiabatic path-integral approach (QUAPI)
of Makri et al. Makri and Makarov (1995a); *QUAPI2, applicable for the study
of subsystem-boson models, and (i) the recently developed influence-functional
path-integral (INFPI) technique Segal  _et al._ (2010); *IF2, able to produce
the dynamics of subsystems in contact with multiple fermi baths. The latter
method (INFPI) essentially generalizes QUAPI. It relies on the observation
that in out-of-equilibrium (and/or finite temperature) situations bath
correlations have a finite range, allowing for their truncation beyond a
memory time dictated by the voltage-bias and the temperature. Taking advantage
of this fact, an iterative-deterministic time-evolution scheme can be
developed, where convergence with respect to the memory length can in
principle be reached.

The principles of the INFPI approach have been detailed in Segal  _et al._
(2010); *IF2, where it has been adopted for investigating dissipation effects
in the nonequilibrium spin-fermion model and charge occupation dynamics in
correlated quantum dots. Recently, it was further utilized for examining the
effect of a magnetic flux on the intrinsic coherence dynamics in a double
quantum dot system Bedkihal and Segal (2012), and for studying relaxation and
equilibration dynamics in finite metal grains Kulkarni  _et al._ (2012);
*Kunal2.

Numerically-exact methodologies are typically limited to simple models;
analytic results are further restricted to specific parameters. The Anderson-
Holstein (AH) model has been studied extensively in this context. In this
model the electronic structure of the molecule is represented by a single
spinless electronic level, with electron occupation on the dot coupled to the
displacement of a single oscillator mode, representing an internal vibration.
This vibration may connect with a secondary phonon bath, representing a larger
phononic environment (internal modes, solvent). The AH model has been
simulated exactly with the secondary phonon bath, using a a real-time path-
integral Monte Carlo approach Muhlbacher and Rabani (2008), and by extending
the multilayer multiconfiguration time-dependent Hartree method to include
fermionic degrees of freedom Wang  _et al._ (2011). More recently, the model
has been simulated by adopting the iterative-summation of path-integral
approach Weiss  _et al._ (2008); Eckel  _et al._ (2010); Hutzen  _et al._
(2012).

In this paper, we examine a variant of the AH model, the Donor (D)-Acceptor
(A) electronic rectifier model Aviram and Ratner (1974). This model
incorporates nonlocal electron-vibration interactions: electronic transitions
between the two molecular states, A and D, are coupled to a particular
internal molecular vibrational mode. Within this simple system, we are
concerned with the development of vibrational instability: Significant
molecular heating can take place once the D level is lifted above the A level,
as the excess electronic energy is used to excite the vibrational mode. This
process may ultimately lead to junction instability and breakdown Lu  _et al._
(2011). We have recently studied a variant of this model (excluding direct D-A
tunneling element), using a Master equation method, by working in the weak
electron-phonon coupling limit. Simine and Segal (2012). An important
observation in that work has been that since the development of this type of
instability is directly linked to the breakdown of the detailed balance
relation above a certain bias (resulting in an enhanced vibrational excitation
rate constant, over relaxation), it suffices to describe the vibrational mode
as a truncated two-level system. In this picture, population inversion in the
two-state system evinces on the development of vibrational instability.

Our objectives here are threefold: (i) To present a numerically-exact
iterative scheme for following the dynamics of a quantum system driven to a
nonequilibrium steady-state due to its coupling to multiple bosonic and
fermionic reservoirs. (ii) To demonstrate the applicability of the method in
the field of molecular electronics. Particularly, to explore the development
of vibrational instability in conducting molecules. (iii) To evaluate the
performance and accuracy of standard-perturbative Master equation treatments,
by comparing their predictions to exact results. Since Master equation
techniques are extensively used for explaining charge transfer phenomenology,
scrutinizing their validity and accuracy is an important task.

The plan of the paper is as follows. In Sec. II we introduce the path-integral
formalism. We describe the iterative time evolution scheme in Sec. III, by
exemplifying it to the case of a spin subsystem. Sec. IV describes a molecular
electronics application, and we follow both electrons and vibrational dynamics
in a dissipative molecular rectifier. Sec. V concludes. For simplicity, we use
the conventions $\hbar\equiv 1$, electron charge $e\equiv 1$, and Boltzmann
constant $k_{B}=1$.

##  II Path-integral formulation

We consider a multi-level subsystem, with the Hamiltonian $H_{S}$, coupled to
multiple bosonic ($B$) and fermionic ($F$) reservoirs that are prepared in an
out-of-equilibrium initial state. The total Hamiltonian $H$ is written as

$\displaystyle H=H_{S}+H_{B}+H_{F}+V_{SB}+V_{SF}.$ (1)

In the energy representation of the isolated subsystem, its Hamiltonian can be
written as

$\displaystyle H_{S}=\sum_{s}\epsilon_{s}|s\rangle\langle s|+\sum_{s\neq
s^{\prime}}v_{s,s^{\prime}}|s\rangle\langle s^{\prime}|.$ (2)

The Hamiltonian $H_{F}$ may comprise of multiple fermionic baths, and
similarly, $H_{B}$ may contain more than a single bosonic reservoir. The terms
$V_{SF}$ and $V_{SB}$ include the coupling of the subsystem to the fermionic
and bosonic environments, respectively. Coupling terms which directly link the
subsystem to both bosonic and fermionic degrees of freedom are not included.
However, $V_{SB}$ and $V_{SF}$ may contain non-additive contributions with
their own set of reservoirs. For example, $V_{SF}$ may admit subsystem
assisted tunneling terms, between separate fermionic baths (metals), see Fig.
1.

We are interested in the time evolution of the reduced density matrix
$\rho_{S}(t)$. This quantity is obtained by tracing the total density matrix
$\rho$ over the bosonic and fermionic reservoirs’ degrees of freedom

$\displaystyle\rho_{S}(t)={\rm Tr}_{B}{\rm
Tr}_{F}\left[e^{-iHt}\rho(0)e^{iHt}\right].$ (3)

We also study the dynamics of certain expectation values, for example, charge
current and energy current. The time evolution of an operator $A$ can be
calculated using the relation

$\displaystyle\langle A(t)\rangle$ $\displaystyle=$ $\displaystyle{\rm
Tr}[\rho(0)A(t)]$ (4) $\displaystyle=$ $\displaystyle\lim_{\lambda\rightarrow
0}\frac{\partial}{\partial\lambda}{\rm Tr}\big{[}\rho(0)e^{iHt}e^{\lambda
A}e^{-iHt}\big{]}.$

Here $\lambda$ is a real number, taken to vanish at the end of the
calculation. When unspecified, the trace is performed over the subsystem
states and all the environmental degrees of freedom. In what follows, we
detail the path-integral approach for the calculation of the reduced density
matrix. Section III.5 presents expressions useful for time-evolving
expectation values of operators.

As in standard path-integral approaches, we decompose the time evolution
operator into a product of $N$ exponentials, $e^{iHt}=\left(e^{iH\delta
t}\right)^{N}$ where $t=N\delta t$, and define the discrete time evolution
operator $\mathcal{G}\equiv e^{iH\delta t}$. Using the Trotter decomposition,
we approximate $\mathcal{G}$ by

$\displaystyle\mathcal{G}\sim\mathcal{G_{F}}\mathcal{G_{B}}\mathcal{G_{S}}\mathcal{G_{B}}\mathcal{G_{F}},$
(5)

where we define

$\displaystyle\mathcal{G_{F}}$ $\displaystyle\equiv$ $\displaystyle
e^{i(H_{F}+V_{SF})\delta t/2},\,\,\,\,\mathcal{G_{B}}\equiv
e^{i(H_{B}+V_{SB})\delta t/2}$ $\displaystyle\mathcal{G_{S}}$
$\displaystyle\equiv$ $\displaystyle e^{iH_{S}\delta t}$ (6)

Note that the breakup of the subsystem-bath term,
$e^{i(H_{F}+V_{SF}+H_{B}+H_{SB})\delta
t/2}\sim\mathcal{G_{B}}\mathcal{G_{F}}$, is exact if the commutator
$[V_{SB},V_{SF}]$ vanishes. This fact allows for an exact separation between
the bosonic and fermionic influence functionals, as we explain below. This
commutator nullifies if the fermionic and bosonic baths couple to commuting
subsystem degrees of freedom, for example, $V_{SB}\propto|s\rangle\langle s|$
and $V_{SF}\propto|s^{\prime}\rangle\langle s^{\prime}|$.

As an initial condition, we assume that at time $t=0$ the subsystem and the
baths are decoupled, $\rho(0)=\rho_{S}(0)\otimes\rho_{B}\otimes\rho_{F}$, and
the baths are prepared in a nonequilibrium (biased) state. For example, we may
include in $H_{F}$ two Fermi seas that are prepared each in a grand canonical
state with different chemical potentials and temperatures. The overall time
evolution can be represented by a path-integral over the subsystem states,

$\displaystyle\langle s_{N}^{+}|\rho_{S}(t)|s_{N}^{-}\rangle$
$\displaystyle=\sum_{s^{\pm}_{0}}\sum_{s_{1}^{\pm}}...\sum_{s_{N-1}^{\pm}}{\rm
Tr}_{B}{\rm Tr}_{F}\Big{[}\langle
s_{N}^{+}|\mathcal{G}^{\dagger}|s^{+}_{N-1}\rangle\langle
s_{N-1}^{+}|\mathcal{G}^{\dagger}|s^{+}_{N-2}\rangle...\langle
s_{0}^{+}|\rho(0)|s^{-}_{0}\rangle...\langle
s_{N-2}^{-}|\mathcal{G}|s^{-}_{N-1}\rangle\langle
s_{N-1}^{-}|\mathcal{G}|s_{N}^{-}\rangle\Big{]}.$ (7)

Here $s_{k}^{\pm}$ represents the discrete path on the forward ($+$) and
backward ($-$) contour. The calculation of each discrete term is done by
introducing four additional summations, e.g.,

$\displaystyle\langle
s_{k}^{-}|{\mathcal{G}}|s_{k+1}^{-}\rangle=\sum_{f^{-}_{k}}\sum_{g^{-}_{k}}\sum_{m^{-}_{k}}\sum_{n^{-}_{k}}\langle
s_{k}^{-}|\mathcal{G_{F}}|f_{k}^{-}\rangle\langle
f_{k}^{-}|\mathcal{G_{B}}|m_{k}^{-}\rangle\langle
m_{k}^{-}|\mathcal{G_{S}}|n_{k}^{-}\rangle\langle
n_{k}^{-}|\mathcal{G_{B}}|g_{k}^{-}\rangle\langle
g_{k}^{-}|\mathcal{G_{F}}|s_{k+1}^{-}\rangle.$ (8)

We substitute Eq. (8) into Eq. (7), further utilizing the factorized
subsystem-reservoirs initial condition as mentioned above, and find that the
function under the sum can be written as a product of separate terms,

$\displaystyle\langle s_{N}^{+}|\rho_{S}(t)|s_{N}^{-}\rangle=\sum_{\bf
s^{\pm}}\sum_{\bf f^{\pm}}\sum_{\bf g^{\pm}}\sum_{\bf m^{\pm}}\sum_{\bf
n^{\pm}}I_{S}({\bf m^{\pm}},{\bf n^{\pm}},s^{\pm}_{0})I_{F}({\bf
s^{\prime\pm}},{\bf f^{\pm}},{\bf g^{\pm}})I_{B}({\bf f^{\pm}},{\bf
m^{\pm}},{\bf n^{\pm}},{\bf g^{\pm}}).$ (9)

Here $I_{S}$ follows the subsystem ($H_{S}$) free evolution. The term $I_{F}$
is referred to as a fermionic “influence functional” (IF), and it contains the
effect of the fermionic degrees of freedom on the subsystem dynamics.
Similarly, $I_{B}$, the bosonic IF, describes how the bosonic degrees of
freedom affect the subsystem. Bold letters correspond to a path, for example,
${\bf m^{\pm}}=\\{m_{0}^{\pm},m_{1}^{\pm},...,m_{N-1}^{\pm}\\}$. We also
define the path ${\bf
s^{\pm}}=\\{s_{0}^{\pm},s_{1}^{\pm},...,s_{N-1}^{\pm}\\}$, and the associate
path which covers $N+1$ points, ${\bf
s^{\prime\pm}}=\\{s_{0}^{\pm},s_{1}^{\pm},...,s_{N-1}^{\pm},s_{N}^{\pm}\\}$.
Given the product structure of Eq. (9), the subsystem, bosonic and the
fermionic terms can be independently evaluated, while coordinating their path.
Explicitly, the elements in Eq. (9) are given by

$\displaystyle I_{S}$ $\displaystyle=$ $\displaystyle\langle
s_{0}^{+}|\rho_{S}(0)|s_{0}^{-}\rangle\Pi_{k=0,...,N-1}\langle
m_{k}^{-}|{\mathcal{G}_{S}}|n_{k}^{-}\rangle\langle
n_{k}^{+}|{\mathcal{G}_{S}^{\dagger}}|m_{k}^{+}\rangle$ $\displaystyle I_{F}$
$\displaystyle=$ $\displaystyle{\rm Tr}_{F}\Big{[}\langle
s_{N}^{+}|{\mathcal{G}_{F}^{\dagger}}|g_{N-1}^{+}\rangle\langle
f_{N-1}^{+}|{\mathcal{G}_{F}^{\dagger}}|s_{N-1}^{+}\rangle...$
$\displaystyle\times$ $\displaystyle\langle
s_{1}^{+}|{\mathcal{G}_{F}^{\dagger}}|g_{0}^{+}\rangle\langle
f_{0}^{+}|{\mathcal{G}_{F}^{\dagger}}|s_{0}^{+}\rangle\rho_{F}\langle
s_{0}^{-}|{\mathcal{G}_{F}}|f_{0}^{-}\rangle\langle
g_{0}^{-}|{\mathcal{G}_{F}}|s_{1}^{-}\rangle...$ $\displaystyle\times$
$\displaystyle\langle s_{N-1}^{-}|{\mathcal{G}_{F}}|f_{N-1}^{-}\rangle\langle
g_{N-1}^{-}|{\mathcal{G}_{F}}|s_{N}^{-}\rangle\Big{]}$ $\displaystyle I_{B}$
$\displaystyle=$ $\displaystyle{\rm Tr}_{B}\Big{[}\langle
g_{N-1}^{+}|{\mathcal{G}_{B}^{\dagger}}|n_{N-1}^{+}\rangle\langle
m_{N-1}^{+}|{\mathcal{G}_{B}^{\dagger}}|f_{N-1}^{+}\rangle...$ (10)
$\displaystyle\times$ $\displaystyle\langle
g_{0}^{+}|{\mathcal{G}_{B}^{\dagger}}|n_{0}^{+}\rangle\langle
m_{0}^{+}|{\mathcal{G}_{B}^{\dagger}}|f_{0}^{+}\rangle\rho_{B}\langle
f_{0}^{-}|{\mathcal{G}_{B}}|m_{0}^{-}\rangle\langle
n_{0}^{-}|{\mathcal{G}_{B}}|g_{0}^{-}\rangle...$ $\displaystyle\times$
$\displaystyle\langle f_{N-1}^{-}|{\mathcal{G}_{B}}|m_{N-1}^{-}\rangle\langle
n_{N-1}^{-}|{\mathcal{G}_{B}}|g_{N-1}^{-}\rangle\Big{]}.$

The dynamics in Eq. (9) can be retrieved by following an iterative scheme, by
using the principles of the INFPI approach Segal  _et al._ (2010); *IF2. In
the next section we illustrate this evolution with a spin subsystem.

##  III Iterative time evolution scheme

We consider here the spin-boson-fermion model. It includes a two-state
subsystem that is coupled through its polarization to bosonic and fermionic
reservoirs. With this relatively simple model, we exemplify the iterative
propagation technique, see Secs. III.1-III.5. Relevant expressions for a
multi-level subsystem and general interaction form are included in Sec. III.6.

###  III.1 spin-boson-fermion model

The spin-fermion model, with a qubit, spin, coupled to a fermionic bath is
kindred to the eminent spin-boson model, describing a qubit interacting with
bosonic environment. It is also related to the Kondo model Kondo (1964), only
lacking direct coupling of the reservoir degrees of freedom to spin-flip
processes. It provides a minimal setting for the study of dissipation and
decoherence effects in the presence of nonequilibrium reservoirs Mitra and
Millis (2005, 2007); Segal  _et al._ (2007); Lutchyn  _et al._ (2008). Here we
put together the spin-boson and the spin-fermion models, and present it in the
general form,

$\displaystyle H_{S}$ $\displaystyle=$
$\displaystyle\Delta\sigma_{x}+B\sigma_{z},\,\,\,\,\,$ $\displaystyle H_{F}$
$\displaystyle=$
$\displaystyle\sum_{j}\epsilon_{j}c_{j}^{\dagger}c_{j}+\sum_{j\neq
j^{\prime}}v_{j,j^{\prime}}^{F}c_{j}^{\dagger}c_{j^{\prime}}$ $\displaystyle
V_{SF}$ $\displaystyle=$
$\displaystyle\sigma_{z}\sum_{j,j^{\prime}}\xi^{F}_{j,j^{\prime}}c_{j}^{\dagger}c_{j^{\prime}}.$
$\displaystyle H_{B}$ $\displaystyle=$
$\displaystyle\sum_{p}\omega_{p}b_{p}^{\dagger}b_{p}+\sum_{p,p^{\prime}}v_{p,p^{\prime}}^{B}b_{p}^{\dagger}b_{p^{\prime}},$
$\displaystyle V_{SB}$ $\displaystyle=$
$\displaystyle\sigma_{z}\sum_{p}\xi_{p}^{B}\left(b_{p}^{\dagger}+b_{p}\right)+\sigma_{z}\sum_{p,p^{\prime}}\zeta_{p,p^{\prime}}^{B}b_{p}^{\dagger}b_{p^{\prime}}.$
(11)

The subsystem includes only two states, with an energy gap $2B$ and a
tunneling splitting $2\Delta$. This minimal subsystem is coupled here through
its polarization to a set of boson and fermion degrees of freedom, where
$\sigma_{z}$ and $\sigma_{x}$ denote the $z$ and $x$ Pauli matrices for a two-
state subsystem, respectively. $b_{p}$ stands for a bosonic operator, to
destroy a mode of frequency $\omega_{p}$, similarly, $c_{j}$ is a fermionic
operator, to annihilate an electron of energy $\epsilon_{j}$ (we assume later
a linear dispersion relation). In this model, spin polarization couples to
harmonic displacements, to scattering events between electronic states in the
metals (fermi reservoirs), and to scattering evens between different modes in
the harmonic bath. Since the commutator between the interaction terms vanish,
$[V_{SF},V_{SB}]=0$, the separation between the bosonic and fermionic IFs is
exact. Moreover, since the fermionic and bosonic operators couple both to
$\sigma_{z}$, we immediately note that $f_{k}^{\pm}=s_{k}^{\pm}$,
$m_{k}^{\pm}=f_{k}^{\pm}$, $n_{k}^{\pm}=g_{k}^{\pm}$ and
$g_{k}^{\pm}=s_{k+1}^{\pm}$. Eq. (9) then simplifies to

$\displaystyle\langle s_{N}^{+}|\rho_{S}(t)|s_{N}^{-}\rangle=\sum_{\bf
s^{\pm}}I_{S}({\bf s^{\prime\pm}})I_{F}({\bf s^{\prime\pm}})I_{B}({\bf
s^{\prime\pm}}),$ (12)

where we recall the definitions of the paths ${\bf
s^{\pm}}=\\{s_{0}^{\pm},s_{1}^{\pm},...,s_{N-1}^{\pm}\\}$ and ${\bf
s^{\prime\pm}}=\\{s_{0}^{\pm},s_{1}^{\pm},...,s_{N-1}^{\pm},s_{N}^{\pm}\\}$.
The subsystem evolution and the IFs are now given by

$\displaystyle I_{S}({\bf s^{\prime\pm}})=\langle
s_{0}^{+}|\rho_{S}(0)|s_{0}^{-}\rangle
K(s_{N}^{\pm},s_{N-1}^{\pm})...K(s_{2}^{\pm},s_{1}^{\pm})K(s_{1}^{\pm},s_{0}^{\pm})$
$\displaystyle I_{B}({\bf s^{\prime\pm}})={\rm
Tr}_{B}\Big{[}e^{-iW_{B}(s_{N}^{+})\delta t/2}e^{-iW_{B}(s_{N-1}^{+})\delta
t}...e^{-iW_{B}(s_{0}^{+})\delta t/2}\rho_{B}e^{iW_{B}(s_{0}^{-})\delta
t/2}....e^{iW_{B}(s_{N-1}^{-})\delta t}e^{iW_{B}(s_{N}^{-})\delta
t/2}\Big{]}.$ $\displaystyle I_{F}({\bf s^{\prime\pm}})={\rm
Tr}_{F}\Big{[}e^{-iW_{F}(s_{N}^{+})\delta t/2}e^{-iW_{F}(s_{N-1}^{+})\delta
t}...e^{-iW_{F}(s_{0}^{+})\delta t/2}\rho_{F}e^{iW_{F}(s_{0}^{-})\delta
t/2}....e^{iW_{F}(s_{N-1}^{-})\delta t}e^{iW_{F}(s_{N}^{-})\delta
t/2}\Big{]},$ (13)

where

$\displaystyle K(s_{k+1}^{\pm},s_{k}^{\pm})=\langle
s_{k+1}^{+}|e^{-iH_{S}\delta t}|s_{k}^{+}\rangle\langle
s_{k}^{-}|e^{iH_{S}\delta t}|s_{k+1}^{-}\rangle$ (14)

is the propagator matrix for the subsystem. We have also used the short
notation $W$ for bath operators that are evaluated along the path,

$\displaystyle W_{F}(s)$ $\displaystyle=$ $\displaystyle H_{F}+\langle
s|V_{SF}|s\rangle,$ $\displaystyle W_{B}(s)$ $\displaystyle=$ $\displaystyle
H_{B}+\langle s|V_{SB}|s\rangle.$ (15)

In the next sections we explain how we compute the bosonic and fermionic IFs.
The former has a closed analytic form in certain situations. The latter is
computed only numerically.

###  III.2 Bosonic IF

We present the structure of the bosonic IF in two separate models,
corresponding to different types of subsystem-boson bath interactions. In both
cases the bosonic bath is prepared in a canonical state of inverse temperature
$\beta_{ph}=1/T_{ph}$,

$\displaystyle\rho_{B}=e^{-\beta_{ph}H_{B}}/{\rm
Tr}_{B}[e^{-\beta_{ph}H_{B}}].$ (16)

Displacement interaction model, $v_{p,p^{\prime}}^{B}=0$ and
$\zeta_{p,p^{\prime}}^{B}=0$. Given the remaining linear displacement-
polarization interaction, an analytic form for the bosonic IF can be written,
the so-called “Feynman-Vernon” influence functional (FV IF) Feynman and Hibbs
(1965). In its time-discrete form, the bosonic IF is given by an exponent with
pairwise interactions along the path Makri and Makarov (1995a); *QUAPI2

$\displaystyle
I_{B}(s_{0}^{\pm},...,s_{N}^{\pm})=\exp\left[-\sum_{k=0}^{N}\sum_{k^{\prime}=0}^{k}(s_{k}^{+}-s_{k}^{-})(\eta_{k,k^{\prime}}s_{k^{\prime}}^{+}-\eta_{k,k^{\prime}}^{*}s_{k^{\prime}}^{-})\right].$
(17)

The coefficients $\eta_{k,k^{\prime}}$ are additive in the number of thermal
baths, and they depend on these baths’ spectral functions and initial
temperatures Makri and Makarov (1995a); *QUAPI2. For completeness, these
coefficients are included in Appendix A.

Boson scattering model, $\xi_{p}^{B}=0$. The bosonic IF can now be computed
numerically, by using the trace formula for bosons Klich (2003)

$\displaystyle{\rm
Tr}_{B}[e^{M_{1}}e^{M_{2}}...e^{M_{k}}]=\det[1-e^{m_{1}}e^{m_{2}}...e^{m_{k}}]^{-1}.$
(18)

Here $m_{k}$ is a single particle operator corresponding to a quadratic
bosonic operator
$M_{k}=\sum_{p,p^{\prime}}(m_{k})_{p,p^{\prime}}b_{p}^{\dagger}b_{p^{\prime}}$.
Application of the trace formula to the bosonic IF (13) leads to

$\displaystyle I_{B}$ $\displaystyle=$ $\displaystyle{\rm
Tr}_{B}[e^{M_{1}}e^{M_{2}}...e^{M_{k}}\rho_{B}]$ (19) $\displaystyle=$
$\displaystyle{\rm
det}\Big{\\{}[\hat{I}_{B}+f_{B}]-e^{m_{1}}e^{m_{2}}...e^{m_{k}}f_{B}\Big{\\}}^{-1}.$

The matrix $\hat{I}_{B}$ is an identity matrix, and the function $f_{B}$
stands for the Bose-Einstein distribution,
$f_{B}=[e^{\beta_{ph}\omega}-1]^{-1}$. The determinant in Eq. (19) can be
evaluated numerically by taking into account $L_{B}$ modes for the boson bath.
This discretization implies a numerical error. Generalizations, to include
more that one bosonic baths, are immediate.

###  III.3 Fermionic IF

The fermionic IF is computed numerically since an exact analytic form is not
known in the general strong coupling limit Mitra and Millis (2005, 2007);
Segal  _et al._ (2007). It is calculated by using the trace formula for
fermions Klich (2003)

$\displaystyle{\rm
Tr}_{F}[e^{M_{1}}e^{M_{2}}...e^{M_{k}}]=\det[1+e^{m_{1}}e^{m_{2}}...e^{m_{k}}].$
(20)

Here $m_{k}$ is a single particle operator corresponding to a quadratic
operator $M_{k}=\sum_{i,j}(m_{k})_{i,j}c_{i}^{\dagger}c_{j}$. In the next
section we consider a model with two Fermi seas, $H_{F}=H_{L}+H_{R}$, prepared
in a factorized state of distinct grand canonical states,
$\rho_{F}=\rho_{L}\otimes\rho_{R}$, with

$\displaystyle\rho_{\nu}=e^{-\beta_{\nu}(H_{\nu}-\mu_{\nu}N_{\nu})}/{\rm
Tr}_{F}[e^{-\beta_{\nu}(H_{\nu}-\mu_{\nu}N_{\nu})}],\,\,\,\,\ \nu=L,R$ (21)

Here $\beta_{\nu}=1/T_{\nu}$ stands for an inverse temperature, and
$\mu_{\nu}$ denotes the chemical potential of the $\nu$ bath. Application of
the trace formula to the fermionic IF in Eq. (13) leads to

$\displaystyle I_{F}$ $\displaystyle=$ $\displaystyle{\rm
Tr}_{F}[e^{M_{1}}e^{M_{2}}...e^{M_{k}}\rho_{F}]$ (22) $\displaystyle=$
$\displaystyle{\rm
det}\Big{\\{}[\hat{I}_{L}-f_{L}]\otimes[\hat{I}_{R}-f_{R}]+e^{m_{1}}e^{m_{2}}...e^{m_{k}}[f_{L}\otimes
f_{R}]\Big{\\}}.$

The matrices $\hat{I}_{\nu}$ are the identity matrices for the $\nu=L,R$
space. The functions $f_{L}$ and $f_{R}$ are the bands electrons’ energy
distribution, $f_{\nu}=[e^{\beta_{\nu}(\epsilon-\mu_{\nu})}+1]^{-1}$. The
determinant in Eq. (22) can be evaluated numerically by taking into account
$L_{s}$ electronic states for each metal. This discretization implies a
numerical error.

###  III.4 The iterative scheme

The dynamics described by Equation (12) includes long-range interactions along
the path, limiting brute force direct numerical simulations to very short
times. The iterative scheme, developed in Ref. Segal  _et al._ (2010); *IF2,
is based on the observation that in standard nonequilibrium situations and at
finite temperatures bath correlations exponentially die Segal  _et al._
(2007); Hutzen  _et al._ (2012), thus the IF can be truncated beyond a memory
time $\tau_{c}=N_{s}\delta t$, corresponding to the time where bath
correlations sustain. Here $N_{s}$ is an integer, $\delta t$ is the
discretized time step, and the correlation time $\tau_{c}$ is dictated by the
bias and temperature. Roughly, for a system under a potential bias $\Delta\mu$
and a temperature $T$, $\tau_{c}\sim\max\\{1/T,\,\,1/\Delta\mu\\}$ Segal  _et
al._ (2010); *IF2. By recursively breaking the IF to include terms only within
$\tau_{c}$, we reach the following (non-unique) structure for the $\alpha=B,F$
influence functional,

$\displaystyle
I_{\alpha}(s_{0}^{\pm},s_{1}^{\pm},s_{2}^{\pm},...,s_{N}^{\pm})\approx
I_{\alpha}(s_{0}^{\pm},s_{1}^{\pm},...,s_{N_{s}}^{\pm})I_{\alpha}^{(N_{s})}(s_{1}^{\pm},s_{2}^{\pm},...,s_{N_{s}+1}^{\pm})I_{\alpha}^{(N_{s})}(s_{2}^{\pm},s_{3}^{\pm},...,s_{N_{s}+2}^{\pm})...$
$\displaystyle\times
I_{\alpha}^{(N_{s})}(s_{N-N_{s}}^{\pm},s_{N-N_{s}+1}^{\pm},...,s_{N}^{\pm}),$
(23)

where we identify the “truncated IF”, $I_{\alpha}^{(N_{s})}$, as the ratio
between two IFs, with the numerator calculated with an additional time step,

$\displaystyle
I_{\alpha}^{(N_{s})}(s_{k},s_{k+1},...,s_{k+N_{s}})=\frac{I_{\alpha}(s_{k}^{\pm},s_{k+1}^{\pm},...,s_{k+N_{s}}^{\pm})}{I_{\alpha}(s_{k}^{\pm},s_{k+1}^{\pm},...,s_{k+N_{s}-1}^{\pm})}.$
(24)

The truncated IF is the central object in our calculations. For fermions, its
numerator and denominator are separately computed using Eq. (22). The bosonic
IF is similarly computed with the help of Eq. (19) when $\xi_{p}^{B}=0$. In
the complementary case, $\zeta_{p,p^{\prime}}^{B}=0$ and
$v_{p,p^{\prime}}^{B}=0$, the truncated-bosonic IF has a closed analytic form:
Using Eq. (17) we find that it comprises only two-body interactions, of
$s_{k+N_{s}}$ with the preceding spins, down to $s_{k}$,

$\displaystyle
I_{B}^{(N_{s})}(s_{k},s_{k+1},...,s_{k+N_{s}})=\exp\left[-\sum_{k^{\prime}=k}^{k+N_{s}}(s_{k+N_{s}}^{+}-s_{k+N_{s}}^{-})(\eta_{k+N_{s},k^{\prime}}s_{k^{\prime}}^{+}-\eta_{k+N_{s},k^{\prime}}^{*}s_{k^{\prime}}^{-})\right].$

Based on the decompositions (24) and (III.4), we time-evolve Eq. (12)
iteratively, by defining a multi-time reduced density matrix
$\tilde{\rho}_{S}(s_{k},s_{k+1},..,s_{k+N_{s}-1})$. Its initial value is given
by

$\displaystyle\tilde{\rho}_{S}(s_{0}^{\pm},...,s_{N_{s}}^{\pm})=I_{S}(s_{0}^{\pm},...,s_{N_{s}}^{\pm})I_{B}(s_{0}^{\pm},...,s_{N_{s}}^{\pm})I_{F}(s_{0}^{\pm},...,s_{N_{s}}^{\pm}).$
(26)

Its evolution is dictated by

$\displaystyle\tilde{\rho}_{S}(s_{k+1}^{\pm},...,s_{k+N_{s}}^{\pm})=\sum_{s_{k}^{\pm}}\tilde{\rho}_{S}(s_{k}^{\pm},...,s_{k+N_{s}-1}^{\pm})K(s_{k+N_{s}}^{\pm},s_{k+N_{s}-1}^{\pm})$
$\displaystyle\times
I_{F}^{(N_{s})}(s_{k}^{\pm},...,s_{k+N_{s}}^{\pm})I_{B}^{(N_{s})}(s_{k}^{\pm},...,s_{k+N_{s}}^{\pm}).$
(27)

The time-local ($t_{k}=k\delta t$) reduced density matrix, describing the
state of the subsystem at a certain time, is reached by summing over all
intermediate states,

$\displaystyle\rho_{S}(t_{k})=\sum_{s_{k-1}^{\pm}...s_{k-N_{s}+1}^{\pm}}\tilde{\rho}_{S}(s_{k-N_{s}+1}^{\pm},...,s_{k}^{\pm}).$
(28)

The bosonic and fermionic IFs may be (and often this is the case)
characterized by different memory time. Thus, in principle we could truncate
the fermionic IF to include $N_{s}^{F}$ terms, and the bosonic IF to include
$N_{s}^{B}$ elements. However, the efficiency of the computation is dictated
by the longest memory time, thus, for convenience, we truncate both IFs using
the largest value, identified by $N_{s}$.

By construction, this iterative approach conserves the trace of the reduced
density matrix, ensuring the stability of the iterative algorithm to long
times Makri and Makarov (1995a); *QUAPI2. This property can be inferred from
Eqs. (12) and (13), by using the formal expressions for the truncated IFs, Eq.
(23) and (24). To prover this property, we trace over the reduced density
matrix at time $t$, identifying $s_{N}=s_{N}^{+}=s_{N}^{-}$,

$\displaystyle{\rm Tr}_{S}[\rho_{S}(t)]$ $\displaystyle\equiv$
$\displaystyle\sum_{s_{N}}\langle s_{N}|\rho_{S}(t)|s_{N}\rangle$
$\displaystyle=$ $\displaystyle\sum_{{\bf s^{\prime\pm}}}I_{S}({\bf
s^{\prime\pm}})I_{F}({\bf s^{\prime\pm}})I_{B}({\bf
s^{\prime\pm}})\delta(s_{N}^{+}-s_{N}^{-})$

Using the cyclic property of the trace, we note that both the fermionic and
bosonic IFs are independent of $s_{N}$, when $s_{N}^{+}=s_{N}^{-}$. Therefore,
the summation over the $s_{N}$ coordinate reduces to a simple sum which can be
performed using the completeness relation for the subsystem states, resulting
in

$\displaystyle\sum_{s_{N}}\langle s_{N}|e^{-iH_{S}\delta
t}|s_{N-1}^{+}\rangle\langle s_{N-1}^{-}|e^{iH_{S}\delta
t}|s_{N}\rangle=\delta(s_{N-1}^{+}-s_{N-1}^{-}).$ (29)

Iterating in this manner we conclude that

$\displaystyle{\rm Tr}_{S}[\rho_{S}(t)]$ $\displaystyle\equiv$
$\displaystyle\sum_{s_{N}}\langle s_{N}|\rho_{S}(t)|s_{N}\rangle$ (30)
$\displaystyle=$ $\displaystyle\sum_{{\bf s^{\prime\pm}}}I_{S}({\bf
s^{\prime\pm}})I_{F}({\bf s^{\prime\pm}})I_{B}({\bf
s^{\prime\pm}})\delta(s_{N}^{+}-s_{N}^{-})\delta(s_{N-1}^{+}-s_{N-1}^{-})...\delta(s_{1}^{+}-s_{1}^{-})\delta(s_{0}^{+}-s_{0}^{-})$
$\displaystyle=$ $\displaystyle\sum_{s_{0}}\langle
s_{0}|\rho_{S}(0)|s_{0}\rangle={\rm Tr}_{S}[\rho_{S}(0)]$

We emphasize that the trace conservation is maintained even with the use of
the truncated form for the IFs. Moreover, it holds irrespective of the details
of the bath and the system-bath interaction form. It is also obeyed in the
more general case, Eq. (9). Equation (27) [and its generalized form, Eq. (34)
below], describe a linear map. Its fixed points are stable if the eigenvalues
of the map have modulus less than one, which is the case here. Thus, our
scheme is expected to approach a stationary-state in the long time limit.

###  III.5 Expectation values for operators

Besides the reduced density matrix, we can also acquire the time evolution of
several expectation values. Adopting the Hamiltonian (11), we illustrate next
how we achieve the charge current behavior. For simplicity, we consider the
case with only two fermionic reservoirs, $\nu=L,R$. The current operator,
e.g., at the $L$ bath is defined as the time derivative of the number
operator. The expectation value of this current is given by

$\displaystyle j_{L}=-\frac{d}{dt}{\rm Tr}[\rho
N_{L}],\,\,\,\,\,N_{L}\equiv\sum_{j\in L}c_{j}^{\dagger}c_{j}$ (31)

We consider the time evolution of the related exponential operator $e^{\lambda
N_{L}}$, with $\lambda$ a real number that is taken to vanish at the end of
the calculation,

$\displaystyle\langle N_{L}(t)\rangle$ $\displaystyle\equiv$
$\displaystyle{\rm Tr}\left[\rho N_{L}(t)\right]$ (32) $\displaystyle=$
$\displaystyle\lim_{\lambda\rightarrow 0}\frac{\partial}{\partial\lambda}{\rm
Tr}\big{[}\rho(0)e^{iHt}e^{\lambda N_{L}}e^{-iHt}\big{]}.$

As before, the initial condition is factorized at $t=0,$
$\rho(0)=\rho_{S}(0)\otimes\rho_{B}\otimes\rho_{F}$. The trace is performed
over subsystem and reservoirs degrees of freedom. By following the same steps
as in Eqs. (3)-(7), we reach the path-integral expression

$\displaystyle\langle e^{\lambda
N_{L}(t)}\rangle=\sum_{s_{0}^{\pm}}\sum_{s_{1}^{\pm}}...\sum_{s_{N-1}^{\pm}}\sum_{s_{N}}{\rm
Tr}_{B}{\rm Tr}_{F}\Big{[}e^{\lambda N_{L}}\langle
s_{N}|\mathcal{G}^{\dagger}|s^{+}_{N-1}\rangle\langle
s_{N-1}^{+}|\mathcal{G}^{\dagger}|s^{+}_{N-2}\rangle...$ (33)
$\displaystyle\times$ $\displaystyle\langle
s_{0}^{+}|\rho(0)|s_{0}^{-}\rangle...\langle
s^{-}_{N-2}|\mathcal{G}|s^{-}_{N-1}\rangle\langle
s_{N-1}^{-}|\mathcal{G}|s_{N}\rangle\Big{]}.$

Factorizing the time evolution operators using Eq. (5), we accomplish the
compact form

$\displaystyle\langle e^{\lambda N_{L}(t)}\rangle$ $\displaystyle=$
$\displaystyle\sum_{\bf s^{\prime\pm}}I_{S}({\bf s^{\prime\pm}})I_{B}({\bf
s^{\prime\pm}})\tilde{I}_{F}({\bf s^{\prime\pm}})\delta(s_{N}^{+}-s_{N}^{-}).$

The terms $I_{S}$ and $I_{B}$ are given in Eq. (13). The fermionic IF
accommodates an additional exponent,

$\displaystyle\tilde{I}_{F}({\bf s^{\prime\pm}})={\rm Tr}_{F}\Big{[}e^{\lambda
N_{L}}e^{-iW_{F}(s_{N}^{+})\delta t/2}e^{-iW_{F}(s_{N-1}^{+})\delta
t}...e^{-iW_{F}(s_{0}^{+})\delta t/2}\rho_{F}e^{iW_{F}(s_{0}^{-})\delta
t/2}....e^{iW_{F}(s_{N-1}^{-})\delta t}e^{iW_{F}(s_{N}^{-})\delta
t/2}\Big{]}.$

We can time evolve the operator $\langle e^{\lambda N_{L}}\rangle$ by using
the iterative scheme of Sec. III.D, by truncating the bosonic and fermionic
IFs up to the memory time $\tau_{c}=N_{s}\delta t$, for several values of
$\lambda$. We then take the numerical derivative with respect to $\lambda$ and
$t$, to attain the charge current itself.

The approach explained here could be used to explore several fermionic
operators, for example, the averaged current $j_{av}=(j_{L}-j_{R})/2$. The
minus sign in front of $j_{R}$ originates from the sign notation, with the
current defined positive when flowing $L$ to $R$. The implementation of a heat
current operator, describing the heat current flowing between two bosonic
reservoirs, requires first the derivation of an analytic form for the bosonic
IF, an expression analogous to the FV IF, and the subsequent time
discretization of this IF, to reach an expression analogous to (17).

###  III.6 Expression for multilevel subsystems and general interactions

So far we have detailed the iterative time evolution scheme for the spin-
boson-fermion model (11). The procedure can be extended, to treat more complex
cases. Based on the general principles outlined in Sec. III.4, one notes that
the path-integral expression (9) can be evaluated iteratively by generalizing
Eq. (27) to the form

$\displaystyle\tilde{\rho}_{S}(v_{k+1}^{\pm},...,v_{k+N_{s}}^{\pm})=$
$\displaystyle\sum_{v_{k}^{\pm}}\tilde{\rho}_{S}(v_{k}^{\pm},...,v_{k+N_{s}-1}^{\pm})K(m_{k+N_{s}}^{\pm},n_{k+N_{s}}^{\pm})I_{F}^{(N_{S})}(s_{k}^{\pm},f_{k}^{\pm},g_{k}^{\pm},...,s_{k+N_{s}}^{\pm},f_{k+N_{s}}^{\pm},g_{k+N_{s}}^{\pm})$
$\displaystyle\times
I_{B}^{(N_{s})}(f_{k}^{\pm},g_{k}^{\pm},m_{k}^{\pm},n_{k}^{\pm},...,f_{k+N_{s}}^{\pm},g_{k+N_{s}}^{\pm},m_{k+N_{s}}^{\pm},n_{k+N_{s}}^{\pm})$
(34)

where we compact several variables,
$v_{k}^{\pm}=\\{s_{k}^{\pm},f_{k}^{\pm},g_{k}^{\pm},m_{k}^{\pm},n_{k}^{\pm}\\}$.
It should be noted that in cases when the IF is time invariant, as in the
molecular electronics case discussed below, one needs to evaluate
$I_{B}^{(N_{s})}$ and $I_{F}^{(N_{s})}$ only once, then use the saved array to
time-evolve the auxiliary density matrix.

Figure 2: Molecular electronic rectifier setup. A biased donor-acceptor
electronic junction is coupled to an anharmonic mode, represented by the two-
state system with vibrational levels $|0\rangle$ and $|1\rangle$. This
molecular vibrational mode may further relax its energy to a phononic thermal
reservoir. This process is represented by a dashed arrow. Direct electron
tunneling element between D and A is depicted by a dotted double arrow. Top:
$\Delta\mu>0$. In our construction both molecular electronic levels are placed
within the bias window at large positive bias, resulting in a large (resonant)
current. Bottom: At negative bias the energy of A is placed outside the bias
window, thus the total charge current is small.

##  IV Application: Molecular Rectifier

The functionality and stability of electron-conducting molecular junctions are
directly linked to heating and cooling effects experienced by molecular
vibrational modes in biased situations Galperin  _et al._ (2007a); Yu.  _et
al._ (2004); Djukic  _et al._ (2005); Kumar  _et al._ (2012); Pasupathy  _et
al._ (2005); Ioffe  _et al._ (2008); Ward  _et al._ (2011); Huang  _et al._
(2007). In particular, junction heating and breakdown may occur once the bias
voltage exceeds typical molecular vibrational frequencies, when the electronic
levels are situated within the bias window, if energy dissipation from the
molecule to its environment is not efficient.

In this section we study the dynamics and steady-state behavior of electrons
and a specific vibrational mode in a molecular conducting junction serving as
an electrical rectifier. The rectifier Hamiltonian is detailed in Sec. IV.1.
In Sec. IV.2 we show that this model can be mapped onto the spin-boson-fermion
Hamiltonian (11). This allows us to employ the path-integral technique of Sec.
III for simulating the rectifier dynamics. The rectification mechanism is
explained in Sec. IV.3. Relevant expressions of a (perturbative) Master
equation method are described in Sec. IV.4, to be compared to our path-
integral based results in Sec. IV.5. Convergence issues and computational
aspects are discussed in Sec. IV.6.

###  IV.1 Rectifier Hamiltonian

The D-A rectifier model includes a biased molecular electronic junction and a
selected (generally anharmonic) internal vibrational mode which is coupled to
an electronic transition in the junction and to a secondary phonon bath,
representing other molecular and environmental degrees of freedom. In the
present study we model the anharmonic mode by a two-state system, and this
model can already capture the essence of the vibrational instability effect
Simine and Segal (2012). For a schematic representation, see Fig. 2. This
model allows us to investigate the exchange of electronic energy with
molecular vibrational heating, and the competition between elastic and
inelastic transport mechanisms. Its close variant has been adopted in Refs.
Entin-Wohlman  _et al._ (2010); Jiang  _et al._ (2012); Entin-Wohlman and
Aharony (2012) for studying the thermopower and thermal transport of electrons
in molecular junctions with electron-phonon interactions, within the linear
response regime.

We assume that the D molecular group is strongly attached to the neighboring
$L$ metal surface, and that this unit is overall characterized by the chemical
potential $\mu_{L}$. Similarly, the A group is connected to the metal $R$,
characterized by $\mu_{R}$. At time $t=0$ the D and A states are put into
contact. Experimentally, the $R$ metal may stand for an STM tip decorated by a
molecular group. This tip is approaching the D site which is attached to the
metal surface $L$. Once the D and A molecular groups are put into contact,
electrons can flow across the junction in two parallel pathways: (i) through a
direct D-A tunneling mechanism, and (ii) inelastically, assisted by a
vibration: excess electron energy goes to excite the D-A vibrational motion,
and vice versa.

The rectifier (rec) Hamiltonian includes the electronic Hamiltonian $H_{el}$
with decoupled D and A states, the vibrational, two-state subsystem $H_{vib}$,
electronic-vibrational coupling $H_{I}$, a free phonon Hamiltonian $H_{ph}$,
and the coupling of this secondary phonon bath to the selected vibration,

$\displaystyle\bar{H}_{rec}=H_{el}+H_{vib}+H_{I}+H_{ph}+H_{vib-ph}.$ (35)

The electronic (fermionic) contribution $H_{el}$ attends for all fermionic
terms besides the direct D and A tunneling term, which for convenience is
included in $H_{I}$,

$\displaystyle H_{el}$ $\displaystyle=$ $\displaystyle
H_{M}+H_{L}^{0}+H_{R}^{0}+H_{C}$ $\displaystyle H_{M}$ $\displaystyle=$
$\displaystyle\epsilon_{d}c_{d}^{\dagger}c_{d}+\epsilon_{a}c_{a}^{\dagger}c_{a}$
$\displaystyle H_{L}^{0}$ $\displaystyle=$ $\displaystyle\sum_{l\in
L}\epsilon_{l}c_{l}^{\dagger}c_{l};\,\,\,\,\,\,\,H_{R}^{0}=\sum_{r\in
R}\epsilon_{r}c_{r}^{\dagger}c_{r}.$ $\displaystyle H_{C}$ $\displaystyle=$
$\displaystyle\sum_{l}v_{l}\left(c_{l}^{\dagger}c_{d}+c_{d}^{\dagger}c_{l}\right)+\sum_{r}v_{r}\left(c_{r}^{\dagger}c_{a}+c_{a}^{\dagger}c_{r}\right).$
(36)

$H_{M}$ stands for the molecular electronic part including two electronic
states, a donor D and an acceptor A. $c_{d/a}^{\dagger}$ ($c_{d/a}$) is a
fermionic creation (annihilation) operator of an electron on the D or A sites,
of energies $\epsilon_{d,a}$. The two metals, $H_{\nu}^{0}$, $\nu=L,R$, are
each composed of a collection of noninteracting electrons. The hybridization
of the D state to the left ($L$) bath, and similarly, the coupling of the A
site to the right ($R$) metal, are described by $H_{C}$. The vibrational
Hamiltonian includes a special nuclear anharmonic vibrational mode of
frequency $\omega_{0}$,

$\displaystyle H_{vib}=\frac{\omega_{0}}{2}\sigma_{z}.$ (37)

The displacement of this mode from equilibrium is coupled to an electron
transition in the system, with an energy cost $\kappa$, resulting in heating
and/or cooling effects,

$\displaystyle
H_{I}=\left(\kappa\sigma_{x}+v_{da}\right)\left(c_{d}^{\dagger}c_{a}+c_{a}^{\dagger}c_{d}\right).$
(38)

Besides the electron-vibration coupling term, $H_{I}$ further includes a
direct electron tunneling element between the D and the A states, of strength
$v_{da}$. Electron transfer between the two metals can therefore proceed
through two mechanisms: coherent tunneling and vibrational-assisted inelastic
transport.

The selected vibrational mode may couple to many other phonons, either
internal to the molecules or external, grouped into a harmonic reservoir,

$\displaystyle H_{ph}$ $\displaystyle=$
$\displaystyle\sum_{p}\omega_{p}b_{p}^{\dagger}b_{p}$ $\displaystyle H_{vib-
ph}$ $\displaystyle=$
$\displaystyle\sigma_{x}\sum_{p}\xi_{p}^{B}\left(b_{p}^{\dagger}+b_{p}\right)$
(39)

The Hamiltonian $H_{vib-ph}$ corresponds to a displacement-displacement
interaction type.

The motivation behind the choice of the two-level system (TLS) mode is
twofold. First, as we showed in Ref. Simine and Segal (2012), the development
of vibrational instability in the D-A rectifier does not depend on the mode
harmonicity, at least in the weak electron-phonon coupling limit. Since it is
easier to simulate a truncated mode with our approach, rather than a harmonic
mode, we settle on the TLS model. Second, while there are many studies where a
perfectly harmonic mode is assumed, for example, see Refs. Muhlbacher and
Rabani (2008); Wang  _et al._ (2011); Hutzen  _et al._ (2012), to the best of
our knowledge our work is the first to explore electron conduction in the
limit of strong vibrational anharmonicity.

###  IV.2 Mapping to the spin-boson-fermion model

We diagonalize the electronic part of the Hamiltonian $H_{el}$ to acquire,
separately, the exact eigenstates for the $L$-half and $R$-half ends of
$H_{el}$,

$\displaystyle H_{el}$ $\displaystyle=$ $\displaystyle H_{L}+H_{R}$
$\displaystyle H_{L}$ $\displaystyle=$
$\displaystyle\sum_{l}\epsilon_{l}a_{l}^{\dagger}a_{l},\,\,\,H_{R}=\sum_{r}\epsilon_{r}a_{r}^{\dagger}a_{r}.$
(40)

Assuming that the reservoirs are dense, their new operators are assigned
energies that are the same as those before diagonalization. The D and A (new)
energies are assumed to be placed within a band of continuous states,
excluding the existence of bound states. The old operators are related to the
new ones by Mahan (2000)

$\displaystyle c_{d}$ $\displaystyle=$
$\displaystyle\sum_{l}\lambda_{l}a_{l},\,\,\,\,\,\,\
c_{l}=\sum_{l^{\prime}}\eta_{l,l^{\prime}}a_{l^{\prime}}$ $\displaystyle
c_{a}$ $\displaystyle=$ $\displaystyle\sum_{r}\lambda_{r}a_{r},\,\,\,\,\,\,\
c_{r}=\sum_{r^{\prime}}\eta_{r,r^{\prime}}a_{r^{\prime}},$ (41)

where the coefficients, e.g., for the $L$ set, are given by

$\displaystyle\lambda_{l}$ $\displaystyle=$
$\displaystyle\frac{v_{l}}{\epsilon_{l}-\epsilon_{d}-\sum_{l^{\prime}}\frac{v_{l^{\prime}}^{2}}{\epsilon_{l}-\epsilon_{l^{\prime}}+i\delta}}$
$\displaystyle\eta_{l,l^{\prime}}$ $\displaystyle=$
$\displaystyle\delta_{l,l^{\prime}}-\frac{v_{l}\lambda_{l^{\prime}}}{\epsilon_{l}-\epsilon_{l^{\prime}}+i\delta}.$
(42)

Similar expressions hold for the $R$ set. It is easy to derive the following
relation,

$\displaystyle\sum_{l^{\prime}}\frac{v_{l^{\prime}}^{2}}{\epsilon_{l}-\epsilon_{l^{\prime}}+i\delta}=PP\sum_{l^{\prime}}\frac{v_{l^{\prime}}^{2}}{\epsilon_{l}-\epsilon_{l^{\prime}}}-i\Gamma_{L}(\epsilon_{l})/2,$
(43)

with the hybridization strength ($v_{j}$ is assumed real),

$\displaystyle\Gamma_{L}(\epsilon)=2\pi\sum_{l}v_{l}^{2}\delta(\epsilon-\epsilon_{l}).$
(44)

With the new operators, the Hamiltonian (35) can be rewritten as

$\displaystyle\bar{H}_{rec}$ $\displaystyle=$
$\displaystyle\sum_{l}\epsilon_{l}a_{l}^{\dagger}a_{l}+\sum_{r}\epsilon_{r}a_{r}^{\dagger}a_{r}+\frac{\omega_{0}}{2}\sigma_{z}$
(45) $\displaystyle+$
$\displaystyle\left(\kappa\sigma_{x}+v_{da}\right)\sum_{l,r}\left[\lambda_{l}^{*}\lambda_{r}a_{l}^{\dagger}a_{r}+\lambda_{r}^{*}\lambda_{l}a_{r}^{\dagger}a_{l}\right]$
$\displaystyle+$
$\displaystyle\sum_{p}\omega_{p}b_{p}^{\dagger}b_{p}+\sigma_{x}\sum_{p}\xi_{p}^{B}\left(b_{p}^{\dagger}+b_{p}\right).$

This Hamiltonian can be transformed into the spin-boson-fermion model of zero
energy spacing, using the unitary transformation

$\displaystyle
U^{\dagger}\sigma_{z}U=\sigma_{x},\,\,\,\,\,U^{\dagger}\sigma_{x}U=\sigma_{z},$
(46)

with $U=\frac{1}{\sqrt{2}}(\sigma_{x}+\sigma_{z})$. The transformed
Hamiltonian $H_{rec}=U^{\dagger}\bar{H}_{rec}U$ includes a $\sigma_{z}$-type
electron-vibration coupling,

$\displaystyle H_{rec}$ $\displaystyle=$
$\displaystyle\sum_{l}\epsilon_{l}a_{l}^{\dagger}a_{l}+\sum_{r}\epsilon_{r}a_{r}^{\dagger}a_{r}+\frac{\omega_{0}}{2}\sigma_{x}$
(47) $\displaystyle+$
$\displaystyle\left(\kappa\sigma_{z}+v_{da}\right)\sum_{l,r}\left[\lambda_{l}^{*}\lambda_{r}a_{l}^{\dagger}a_{r}+\lambda_{r}^{*}\lambda_{l}a_{r}^{\dagger}a_{l}\right]$
$\displaystyle+$
$\displaystyle\sum_{p}\omega_{p}b_{p}^{\dagger}b_{p}+\sigma_{z}\sum_{p}\xi_{p}^{B}\left(b_{p}^{\dagger}+b_{p}\right).$

It describes a spin (TLS) coupled diagonally to two fermionic environments and
to a single boson bath. One can immediately confirm that this Hamiltonian is
accounted for by Eq. (11). To simplify our notation, we further identify the
electronic-vibration effective coupling parameter

$\displaystyle\xi_{l,r}^{F}=\kappa\lambda_{l}^{*}\lambda_{r}.$ (48)

For later use we also define the spectral function of the secondary phonon
bath as

$\displaystyle
J_{ph}(\omega)=\pi\sum_{p}(\xi^{B}_{p})^{2}\delta(\omega-\omega_{p}).$ (49)

In our simulations below we adopt an ohmic function,

$\displaystyle J_{ph}(\omega)=\frac{\pi K_{d}}{2}\omega
e^{-\omega/\omega_{c}},$ (50)

with the dimensionless Kondo parameter $K_{d}$, characterizing subsystem-bath
coupling, and the cutoff frequency $\omega_{c}$.

As an initial condition for the reservoirs, we assume canonical distributions
with the boson-phonon bath distribution following
$\rho_{B}=e^{-\beta_{ph}H_{ph}}/{\rm Tr}_{B}[e^{-\beta_{ph}H_{ph}}]$ and the
electronic-fermionic initial density matrix obeying
$\rho_{F}=\rho_{L}\otimes\rho_{R}$, with
$\rho_{\nu}=e^{-\beta_{\nu}(H_{\nu}-\mu_{\nu}N_{\nu})}/{\rm
Tr}_{F}[e^{-\beta_{\nu}(H_{\nu}-\mu_{\nu}N_{\nu})}]$, $\nu=L,R$. This results
in the expectation values of the exact eigenstates,

$\displaystyle\langle
a_{l}^{\dagger}a_{l^{\prime}}\rangle=\delta_{l,l^{\prime}}f_{L}(\epsilon_{l}),\,\,\,\,\,\langle
a_{r}^{\dagger}a_{r^{\prime}}\rangle=\delta_{r,r^{\prime}}f_{R}(\epsilon_{r}),$
(51)

where $f_{L}(\epsilon)=[\exp(\beta_{L}(\epsilon-\mu_{L}))+1]^{-1}$ denotes the
Fermi distribution function. An analogous expression holds for
$f_{R}(\epsilon)$. The reservoirs temperatures are denoted by $1/\beta_{\nu}$;
the chemical potentials are $\mu_{\nu}$.

Figure 3: Left panel: Energy of the donor (full line) and acceptor states
(dashed line). The dotted lines correspond to the chemical potentials at the
left and right sides. Right panel: Damping rate $K_{vib}$. The junction’s
parameters are $\Gamma_{\nu}=1$, $\beta_{\nu}=200$, $\kappa=0.1$,
$\omega_{0}=0.2$, and $\epsilon_{d}(\Delta\mu=0)=-0.2$,
$\epsilon_{a}(\Delta\mu=0)=0.4$. We used fermionic metals with a linear
dispersion relations for the original $H_{\nu}^{0}$ baths and sharp cutoffs at
$\pm 1$. All energy parameters are given in units of eV.

Figure 4: Scheme of the vibrational mode excitation and relaxation processes.
A full circle represents an electron transferred; a hollow circle depicts the
hole that has been left behind.

###  IV.3 Rectifying mechanism

We now explain the operation principles of the molecular rectifier. In our
construction the application of a bias voltage linearly shifts the energies of
the molecular electronic levels, D and A. In equilibrium, we set
$\epsilon_{a}<0$ and $\epsilon_{d}>0$. Under positive bias, defined as
$\mu_{L}-\mu_{R}>0$, the energy of the acceptor level increases, and the donor
level drops down, see Fig. 2. When both levels are buried within the bias
window, the junction can support large currents. At negative bias the
electronic level A is positioned above the bias window, resulting in small
currents. For a scheme of the energy organization of the system, see Fig. 3
panel, left panel.

A generic mechanism leading to vibrational instabilities (and eventually
junction rupture) in D-A molecular rectifiers has been discussed in Ref. Lu
_et al._ (2011): At large positive bias, when the D state is positioned above
the acceptor level, electron-hole pair excitations by the molecular vibration
(TLS) dominate the mode dynamics. This can be schematically seen in Fig. 4.
The second-order perturbation theory rate constant, to excite the vibrational
mode, while transferring an electron from $L$ to $R$, $k_{0\rightarrow
1}^{L\rightarrow R}$, overcomes other rates once the density of states at the
left end is positioned above the density of states at the right side. This is
the case at large positive bias, given our construction. The rate
$k_{0\rightarrow 1}^{L\rightarrow R}$ is defined next, in Sec. IV.4.

###  IV.4 Master equation ($v_{da}=0$)

In the limit of weak electron-vibration coupling, once the direct tunneling
term is neglected, $v_{da}$=0, it can be shown that the population of the
truncated vibrational mode satisfies a kinetic equation Simine and Segal
(2012),

$\displaystyle\dot{p}_{1}=-\left(k_{1\rightarrow 0}^{e}+k_{1\rightarrow
0}^{b}\right)p_{1}+\left(k_{0\rightarrow 1}^{e}+k_{0\rightarrow
1}^{b}\right)p_{0},$ $\displaystyle p_{0}+p_{1}=1.$ (52)

The excitation ($k_{0\rightarrow 1}$) and relaxation ($k_{1\rightarrow 0}$)
rate constants are given by a Fourier transform of bath correlation functions
of the operators $F_{e}$ and $F_{b}$, defined as

$\displaystyle F_{e}$ $\displaystyle=$
$\displaystyle\sum_{l,r}(\xi_{l,r}^{F}a_{l}^{\dagger}a_{r}+\xi_{r,l}^{F}a_{r}^{\dagger}a_{l}),$
$\displaystyle F_{b}$ $\displaystyle=$
$\displaystyle\sum_{p}\xi_{p}^{B}(b_{p}^{\dagger}+b_{p}),$ (53)

to yield

$\displaystyle k_{s\rightarrow s^{\prime}}^{e}$ $\displaystyle=$
$\displaystyle\int_{-\infty}^{\infty}e^{i(\epsilon_{s}-\epsilon_{s^{\prime}})\tau}{\rm
Tr}_{F}\left[\rho_{F}F_{e}(\tau)F_{e}(0)\right]d\tau$ $\displaystyle
k_{s\rightarrow s^{\prime}}^{b}$ $\displaystyle=$
$\displaystyle\int_{-\infty}^{\infty}e^{i(\epsilon_{s}-\epsilon_{s^{\prime}})\tau}{\rm
Tr}_{B}\left[\rho_{B}F_{b}(\tau)F_{b}(0)\right]d\tau.$ (54)

Here $s=0,1$ and $\epsilon_{1}-\epsilon_{0}=\omega_{0}$. The operators are
given in the interaction representation, e.g.,
$a_{l}^{\dagger}(t)=e^{iH_{L}t}a_{l}^{\dagger}e^{-iH_{L}t}$.

Phonon-bath induced rates. Expression (54) can be simplified, and the
contribution of the phonon bath to the vibrational rates reduces to

$\displaystyle k_{1\rightarrow 0}^{b}$ $\displaystyle=$
$\displaystyle\Gamma_{ph}(\omega_{0})[f_{B}(\omega_{0})+1],$ $\displaystyle
k_{0\rightarrow 1}^{b}$ $\displaystyle=$ $\displaystyle k_{1\rightarrow
0}^{b}e^{-\omega_{0}\beta_{ph}},$ (55)

where $f_{B}(\omega)=[e^{\beta_{ph}\omega}-1]^{-1}$ denotes the Bose-Einstein
distribution function. The damping rate is defined as
$\Gamma_{ph}(\omega)=2J_{ph}(\omega)$,

$\displaystyle\Gamma_{ph}(\omega)=2\pi\sum_{p}(\xi_{p}^{B})^{2}\delta(\omega_{p}-\omega).$
(56)

For brevity, we ignore below the direct reference to frequency.

Electronic-baths induced rates. The electronic rate constants (54) include the
following contributions Simine and Segal (2012),

$\displaystyle k_{1\rightarrow 0}^{e}=k_{1\rightarrow 0}^{L\rightarrow
R}+k_{1\rightarrow 0}^{R\rightarrow L};\,\,\,\ k_{0\rightarrow
1}^{e}=k_{0\rightarrow 1}^{L\rightarrow R}+k_{0\rightarrow 1}^{R\rightarrow
L},$ (57)

satisfying

$\displaystyle k_{1\rightarrow 0}^{L\rightarrow R}$ $\displaystyle=$
$\displaystyle
2\pi\kappa^{2}\sum_{l,r}|\lambda_{l}|^{2}|\lambda_{r}|^{2}f_{L}(\epsilon_{l})(1-f_{R}(\epsilon_{r}))\delta(\omega_{0}+\epsilon_{l}-\epsilon_{r})$
$\displaystyle k_{0\rightarrow 1}^{L\rightarrow R}$ $\displaystyle=$
$\displaystyle
2\pi\kappa^{2}\sum_{l,r}|\lambda_{l}|^{2}|\lambda_{r}|^{2}f_{L}(\epsilon_{l})(1-f_{R}(\epsilon_{r}))\delta(-\omega_{0}+\epsilon_{l}-\epsilon_{r}).$
(58)

Similar relations hold for the right-to-left going excitations. The energy in
the Fermi function $f_{\nu}(\epsilon)$ is measured with respect to the
(equilibrium) Fermi energy, placed at $(\mu_{L}+\mu_{R})$, and we assume that
the bias is applied symmetrically, $\mu_{L}=-\mu_{R}$. The rates can be
expressed in terms of the fermionic $\nu=L,R$ spectral density functions

$\displaystyle J_{\nu}(\epsilon)$ $\displaystyle=$ $\displaystyle
2\pi\kappa\sum_{j\in\nu}|\lambda_{j}|^{2}\delta(\epsilon_{j}-\epsilon).$ (59)

Using Eq. (42) we resolve this as a Lorentzian function, centered around
either the D or the A level,

$\displaystyle J_{L}(\epsilon)$ $\displaystyle=$
$\displaystyle\kappa\frac{\Gamma_{L}(\epsilon)}{(\epsilon-\epsilon_{d})^{2}+\Gamma_{L}(\epsilon)^{2}/4}$
$\displaystyle J_{R}(\epsilon)$ $\displaystyle=$
$\displaystyle\kappa\frac{\Gamma_{R}(\epsilon)}{(\epsilon-\epsilon_{a})^{2}+\Gamma_{R}(\epsilon)^{2}/4}.$
(60)

The electronic hybridization $\Gamma_{\nu}(\epsilon)$ is given in Eq. (44).
Using these definitions, we express the electronic rates [Eq. (58)] by
integrals ($s,s^{\prime}$=0,1)

$\displaystyle k_{s\rightarrow s^{\prime}}^{\nu\rightarrow\nu^{\prime}}$
$\displaystyle=$
$\displaystyle\frac{1}{2\pi}\int_{-\infty}^{\infty}f_{\nu}(\epsilon)\left[1-f_{\nu^{\prime}}(\epsilon+(s-s^{\prime})\omega_{0})\right]J_{\nu}(\epsilon)J_{\nu^{\prime}}(\epsilon+(s-s^{\prime})\omega_{0})d\epsilon.$
(61)

Observables. Within this simple kinetic approach, junction stability can be
recognized by watching the TLS population in the steady-state limit:
population inversion reflects on vibrational instability Simine and Segal
(2012). Solving Eq. (52) in the long time limit we find that

$\displaystyle p_{1}=\frac{k_{0\rightarrow 1}^{e}+k_{0\rightarrow
1}^{b}}{k_{0\rightarrow 1}^{e}+k_{0\rightarrow 1}^{b}+k_{1\rightarrow
0}^{e}+k_{1\rightarrow 0}^{b}},\,\,\,\,\,p_{0}=1-p_{1}.$ (62)

A related measure is the damping rate $K_{vib}$ Lu  _et al._ (2011), depicted
in Fig. 3 panel (b). It is defined as the difference between relaxation and
excitation rates,

$\displaystyle K_{vib}\equiv k_{1\rightarrow 0}^{e}+k_{1\rightarrow
0}^{b}-\left(k_{0\rightarrow 1}^{e}+k_{0\rightarrow 1}^{b}\right).$ (63)

Positive $K_{vib}$ indicates on a “normal” thermal-like behavior, when
relaxation processes overcome excitations. In this case, the junction remains
stable in the sense that the population of the ground state is larger than the
population of the excited level. A negative value for $K_{vib}$ evinces on the
process of an uncontrolled heating of the molecular mode, eventually leading
to vibrational instability and junction breakdown.

In the steady-state limit, the charge current $j$, flowing from $L$ to $R$, is
given by Simine and Segal (2012)

$\displaystyle j=p_{1}\left(k_{1\rightarrow 0}^{L\rightarrow
R}-k_{1\rightarrow 0}^{R\rightarrow L}\right)+p_{0}\left(k_{0\rightarrow
1}^{L\rightarrow R}-k_{0\rightarrow 1}^{R\rightarrow L}\right).$ (64)

This relation holds even when the TLS is coupled to an additional boson bath.
Note that in the long time limit the current that is evaluated at the left end
$j_{L}$ is equal to $j_{R}$. Therefore, we simple denote the current by $j$ in
that limit.

Master equation calculations proceed as follows. We set the hybridization
energy $\Gamma_{\nu}$ as an energy independent parameter, and evaluate the
fermionic spectral functions $J_{\nu}(\epsilon)$ of Eq. (60). With this at
hand, we integrate (numerically) Eq. (61), and gain the fermionic-bath induced
rates. The phonon bath-induced rates (55) are reached by setting the
parameters of the spectral function $J_{ph}$, to directly obtain
$\Gamma_{ph}$, see Eq. (56). Using this set of parameters, we evaluate the
levels occupation and the charge current directly in the steady-state limit.
We can also time evolve the set of differential equations (52), to obtain the
trajectory $p_{1,0}(t)$.

Figure 5: Absolute value of the quantity $\pi\rho\xi_{l,r}^{F}$. The figure
was generated by discretizing the reservoirs, using bands extending from $-D$
to $D$, $D=1$, with $N_{L}=200$ states per each band a linear dispersion
relation and a constant density of states for the $H_{L,R}^{0}$ reservoirs,
with a constant density of states $\rho=N_{L}/2D$. Electron-vibration coupling
is given by $\kappa=0.1$.

###  IV.5 Results

We simulate the dynamics of the subsystem in the spin-boson-fermion
Hamiltonian (47) using the path-integral approach of Sec. III. In order to
retrieve the vibrational mode occupation in the original basis in which Eq.
(45) is written, we rotate the reduced density matrix $\rho_{S}(t)$ back to
the original basis by applying the transformation
$U=\frac{1}{\sqrt{2}}(\sigma_{x}+\sigma_{z})$,

$\displaystyle\bar{\rho}_{S}(t)=U\rho_{S}(t)U.$ (65)

The diagonal elements of $\bar{\rho}_{S}(t)$, correspond to the vibrational
mode occupation, the ground state $|0\rangle$ and the excited state
$|1\rangle$,

$\displaystyle p_{0}(t)=\langle
0|\bar{\rho}_{S}(t)|0\rangle\,\,\,\,\,\,p_{1}(t)=\langle
1|\bar{\rho}_{S}(t)|1\rangle.$ (66)

As an initial condition we usually take
$\rho_{S}(0)=\frac{1}{2}(-\sigma_{x}+\hat{I}_{s})$, $\hat{I}_{s}$ is a
$2\times 2$ unit matrix. Under this choice, $\bar{\rho}_{S}(0)$ has only its
ground state populated.

Our simulations are performed with the following setup, displayed in the left
panel of Fig. 3: In the absence of a bias voltage we assign the donor the
energy $\epsilon_{d}=-0.2$ and the acceptor the value $\epsilon_{a}=0.4$.
These molecular electronic states are assumed to linearly follow the bias
voltage. The right panel in Fig. 3 depicts the damping rate $K_{vib}$ in the
absence of coupling to the phonon bath, as evaluated using the Master equation
method. This measure becomes negative beyond $\Delta\mu\sim 0.85$, which
corresponds to the situation where the (bias shifted) donor energy exceeds the
acceptor by $\omega_{0}$, $\epsilon_{d}-\epsilon_{a}\gtrsim\omega_{0}$;
$\omega_{0}=0.2$. This results in a significant exchange of electronic energy
to heat, affecting junction’s instability.

####  IV.5.1 Isolated mode

We study the time evolution of the vibrational mode occupation using
$v_{da}=0$ (unless otherwise stated), further decoupling it from a secondary
phonon bath, $K_{d}$=0.

Electron-vibration interaction energy. The interaction energy of the subsystem
(TLS) to the electronic degrees of freedom is encapsulated in the matrix
elements $\xi_{l,r}^{F}\equiv\kappa\lambda_{l}^{*}\lambda_{r}$, see Eq. (48).
The strength of this interaction is measured by the dimensionless parameter
$\pi\rho(\epsilon_{F})\xi_{l,r}^{F}$, which connects to the phase shift
experienced by Fermi sea electrons due to a scattering potential, introduced
here by the vibrational mode Roulet  _et al._ (1969). Here,
$\rho(\epsilon_{F})$ stands for the density of states at the Fermi energy.
Using the parameters of Fig. 3, taking $\kappa=0.1$, we show the absolute
value of these matrix elements in Fig. 5. The contour plot is mostly limited
to values smaller than 0.1, thus we conclude that this set of parameters
correspond to the weak coupling limit Roulet  _et al._ (1969). In this limit,
path-integral simulations should agree with Master equation calculations, as
we indeed confirm below. Deviations should be expected at larger values,
$\kappa\gtrsim 0.2$, and we study below these cases.

Units. We perform the simulations in arbitrary units with $\hbar\equiv 1$. One
can scale all energies with respect to the molecule-metal hybridization
$\Gamma_{\nu}$. With $\Gamma_{\nu}=1$, the weak coupling limit covers
$\kappa/\Gamma_{\nu}\lesssim 0.2$. To present results in physical units, we
assume that all energy parameters are given in eV, and scale correspondingly
the time unit and currents.

Figure 6: Population dynamics and convergence behavior of the truncated and
isolated vibrational mode (TLS) with increasing $N_{s}$. (a)-(b) Stable
behavior at $\mu_{L}=-\mu_{R}=0.2$. (c)-(d) Population inversion at
$\mu_{L}=-\mu_{R}=0.6$. Other parameters are the same as in Fig. 3. In all
figures $\delta t$=1, $N_{s}=3$ (heavy dotted), $N_{s}=4$ (heavy dashed),
$N_{s}=5$ (dashed-dotted), $N_{s}=6$ (dotted), $N_{s}=7$ (dashed) and
$N_{s}=8$ (full). We used $L_{s}=30$ electronic states at each fermionic bath
with sharp cutoffs at $\pm 1$.

Figure 7:  Independence of the population $p_{0}$ on the initial state for
different biases, $\Delta\mu=$ 0.4, 0.8, 1.2 top to bottom. Other parameters
are the same as in Fig. 3 and Fig. 6.

Figure 8: Population dynamics, $p_{0}(t)$. (a) Comparison between exact
simulations (dashed) and Master equation results (dashed-dotted) at
$\kappa=0.2$. (b) Deviations between exact results and Master equations for
$\kappa=0.1$ (dot and $\circ$) and for $\kappa=0.2$ ($+$ and $x$). Other
parameters are as determined in Fig. 3.

Dynamics. We first focus on two representative values for the bias voltage: In
the low-positive bias limit a stable operation is expected, reflected by a
normal population, $p_{0}>p_{1}$. At large positive bias population inversion
may take place, indicating on the onset of instability and potential junction
rupture Simine and Segal (2012).

Fig. 6 displays the TLS dynamics, and we present data for different memory
sizes $N_{s}\delta t$. At small positive bias,
$\epsilon_{d}-\epsilon_{a}<\omega_{0}$, the mode occupation is “normal”,
$p_{0}>p_{1}$. In particular, in panels (a)-(b) we discern the case
$\mu_{L}=-\mu_{R}=0.2$, resulting in the (shifted) electronic energies
$\epsilon_{d}=0$ and $\epsilon_{a}=0.2$. In this case the (converged)
asymptotic long-time population (representing steady-state values), are
$p_{0}^{ss}=0.76$ and $p_{1}^{ss}=0.24$. In contrast, when the bias is large,
$\mu_{L}=-\mu_{R}=0.6$, the electronic levels are shifted to
$\epsilon_{d}=0.4$ and $\epsilon_{a}=-0.2$, and electrons crossing the
junction discard their excess energy into the vibrational mode. Indeed, we see
in Fig. 6(c)-(d) the process of population inversion, $p_{0}^{ss}=0.43$ and
$p_{1}^{ss}=0.57$. The TLS approaches the steady-state value around
$t_{ss}\sim 0.1$ ps. Regarding convergence behavior, we note that at large
bias convergence is reached with a shorter memory size, compared to the small
bias case, as expected Segal  _et al._ (2010).

Fig. 7 exhibits the dynamics with different initial conditions, demonstrating
that the steady-state value is identical, yet the timescale to reach the
stationary limit may depend on the initial state.

We compare the exact dynamics to the Master equation time evolution behavior,
reached by solving Eq. (52). Panel (a) in Fig. 8 demonstrates excellent
agreement for $\kappa=0.2$, for both positive and negative biases. Below we
show that at this value Master equation’s predictions for the charge current
deviate from the exact result. Panel (b) in Fig. 8 focuses on the departure of
Master equation data from the exact values. These deviations are small, but
their dynamics indicate on the existence of high order excitation and
relaxation rates, beyond the second order rates of Sec. IV.4.

Steady-state characteristics. The full bias scan of the steady-state
population is displayed in Fig. 9, and we compare path-integral results with
Master equation calculations, revealing an excellent agreement in this weak
coupling limit ($\kappa=0.1$). The convergence behavior is presented in Fig.
10, and we plot the steady-state values as a function of memory size
($\tau_{c}$) for three different time steps, for representative biases. The
path-integral results well converge at intermediate-to-large positive biases,
$\Delta\mu\gtrsim 0.2$. We had difficulty converging our results in two
domains: (i) At small-positive potential bias, $\Delta\mu<0.2$. Here, large
memory size should be used for reaching full convergence; decorrelation time
approximately scales with $1/\Delta\mu$. (ii) At large negative biases,
$\Delta\mu<-0.4$ the current is very small as we show immediately. This
implies poor convergence at the range of $\tau_{c}$ employed. At these
negative biases the data oscillates with $\tau_{c}$, thus at negative bias it
is the averaged value for several-large $\tau_{c}$ which is plotted in Fig. 9.

Charge current. We show the current characteristics in Fig. 11, and confirm
that the junction acts as a charge rectifier. The insets display transient
data, affirming that at large bias steady-state is reached faster than in the
low bias case.

Strong coupling. Results at weak-to-strong couplings are shown in Fig. 12. The
value of the current, as reached from Master equation calculations, scale with
$\kappa^{2}$. In contrast, exact simulations indicate that the current grows
more slowly with $\kappa$, and it displays clear deviations (up to $50\%$)
from the perturbative Master equation result at $\kappa=0.3$. Interestingly,
the vibrational occupation (inset) shows little sensitivity to the coupling
strength, and even at $\kappa=0.3$ the Master equation technique provides an
excellent estimation for the levels occupation. This could be reasoned by the
fact that excited levels occupation is given by ratio of excitation rates to
the sum of excitation and relaxation rates. Such a ratio is (apparently) only
weakly sensitive to the value of $\kappa$ itself, even when high-order
processes do contribute to the current.

Figure 9: Converged data for the population of the isolated vibrational mode
in the steady-state limit with $\kappa=0.1$. Other parameters are the same as
in Fig. 3. We display path-integral data for $p_{0}$ ($\circ$) and $p_{1}$
($\square$). Master equation results appear as dashed line for $p_{0}$ and
dashed-dotted line for $p_{1}$.

Figure 10: (a) Convergence behavior of the population $p_{0}$ in the steady-
state limit for $\kappa=0.1$. Other parameters are the same as in Fig. 3.
Plotted are the steady-state values using different time steps, $\delta t=0.8$
($\circ$), $\delta t=1.0$ ($\square$), and $\delta t=1.2$ ($\diamond$) at
different biases, as indicated at the right end. (b) Population mean and its
standard deviation, utilizing the last six points from panel (a). (c) Current
mean and its standard deviation, similarly attained from the data in panel
(a).

Figure 11: Charge current in the steady-state limit for $\kappa=0.1$,
$K_{d}=0$. Other parameters are the same as in Fig. 3. Path-integral data
($\circ$), Master equation results (dashed). The insets display transient
results at $\Delta\mu=1.0$ eV (top) and $\Delta\mu=-0.5$ eV (bottom).

Figure 12: Charge current and vibrational occupation in the steady-state limit
at different electron-vibration coupling. Path-integral data is marked by
symbols, $\kappa=0.1$ ($\circ$), $\kappa=0.2$ ($\diamond$) and $\kappa=0.3$
($\square$). The corresponding Master equation results appear as dashed lines.
Inset: The population behavior in the steady-state limit for the three cases
$\kappa=0.1$ ($\circ$), $\kappa=0.2$ ($\diamond$) and $\kappa=0.3$
($\square$), with empty symbols for $p_{0}$ and filled ones for $p_{1}$. Other
parameters are the same as in Fig. 3.

Direct tunneling vs. vibrational assisted transport. Until this point (and
beyond this subsection) we have taken $v_{da}=0$. We now evaluate the
contribution of different transport mechanisms by adding a direct D-A
tunneling term, $v_{da}\neq 0$ to our model Hamiltonian. Electrons can now
either cross the junction in a coherent manner, or inelastically, by
exciting/de-exciting the vibrational mode. Fig. 13 demonstrates that when the
vibration assisted transport energy $\kappa$ is identical in strength to the
direct tunneling element $v_{da}$, the overall current is enhanced by about a
factor of two, compared to the case when only vibrational-assisted processes
are allowed. We also note that the occupation of the vibrational mode is
barely affected by the opening of the new electron transmission route
(deviations are within the convergence error). While we compare IF data to
Master equation results when $v_{da}=0$, in the general case of a nonzero D-A
tunneling term perturbative methods are more involved, and techniques similar
to those developed for the AH model should be used Galperin  _et al._ (2007a);
Mitra  _et al._ (2004); Galperin  _et al._ (2007b, 2009); Fransson and
Galperin (2010); Leijnse and Wegewijs (2008); Hartle and Thoss (2011a, b).

Figure 13: Study of the contribution of different transport mechanisms.
$v_{da}=0$ ($\circ$), with Master equation results noted by the dashed line,
and $v_{da}=0.1$ ($\square$). The main plot displays the charge current. The
inset presents the vibrational levels occupation, with empty symbols for
$p_{0}$ and filled symbols for $p_{1}$. Other parameters are the same as in
Fig. 3, particularly, the vibrational-electronic coupling is $\kappa=0.1$.

####  IV.5.2 Equilibration with a secondary phonon bath

We couple the isolated-truncated vibrational mode to a secondary phonon bath,
and follow the mode equilibration with this bath and the removal of the
vibrational instability effect, as we increase the vibrational mode-phonon
bath coupling. As an initial condition, the boson bath is assumed to be
thermal with an inverse temperature $\beta_{ph}$. This bath is characterized
by an ohmic spectral function (50) with the dimensionless Kondo parameter
$K_{d}$, characterizing subsystem-bath coupling, and the cutoff frequency
$\omega_{c}$.

Population behavior. We follow the mode dynamics to the steady-state limit
using the path-integral approach of Sec. III. The bosonic IF is given in the
appendix. We compare exact results with Master equation predictions, and Fig.
14 depicts our simulations. The following observations can be made: (i) The
vibrational instability effect is removed already for $K_{d}=0.01$, though
nonequilibrium effects are still largely visible in the mode occupation. (ii)
The vibrational mode is closed to be equilibrated with the phonon bath once
$K_{d}\sim 0.1$. (iii) For the present range of parameters (large
$\omega_{c}$, weak subsystem-bath couplings), Master equation tools reproduce
the behavior of the vibrational mode.

Charge Current. The role of the secondary phonon bath on the charge current
characteristics is displayed in Fig. 15. There are two main effects related to
the presence of the phonon bath: The step structure about zero bias is
flattened when $K_{d}\sim 0.1$, and the current-voltage characteristics as a
whole is slightly enhanced at finite $K_{d}$, at large bias. Both of these
effects are excellently reproduced with the Master equation, and we conclude
that in this weak-coupling regime the presence of the phonon bath does not
affect the rectifying behavior of the junction. We have also verified (not
shown) that at stronger coupling, $\kappa=0.2$ (where Master equation fails),
the thermal bath similarly affects the current-voltage behavior.

An important observation is that the current itself does not testify on the
state of the vibrational mode, whether it is in a stable or an unstable
nonequilibrium state, and whether it is thermalized. The study of the current
characteristics itself ($j$ vs. $\Delta\mu$) is therefore insufficient to
determine junction stability. More detailed information can be gained from the
structure of the first derivative, $dj/d(\Delta\mu)$, the local density of
states, and the second derivative, $d^{2}j/d\Delta\mu^{2}$, providing spectral
features Mii  _et al._ (2003); Galperin  _et al._ (2004a, b). In order to
examine these quantities, our simulations should be performed with many more
bath states, to eliminate possible spurious oscillations in the current (of
small amplitudes) that may result from the finite discretization of the fermi
baths.

Figure 14: Equilibration of the molecular vibrational mode with increasing
coupling to a secondary phonon bath. Path-integral results, (full symbols for
$p_{1}$, and empty symbols for $p_{0}$) with $K_{d}=0$ ($\circ$), $K_{d}=0.01$
($\diamond$), $K_{d}=0.1$ ($\square$), and, $K_{d}=0.1$, $\kappa=0$
($\triangleleft$). Unless otherwise specified, $\kappa=0.1$, $\beta_{ph}=5$
and the spectral function follows (50) with $\omega_{c}$=15. All other
electronic parameters are the same as in Fig. 3. Master equation results
appear in dotted lines.

Figure 15: Charge current for an isolated mode,$K_{d}=0$ ($\circ$), and an
equilibrated mode, $K_{d}=0.1$, $\beta_{ph}$=5, $\omega_{c}$=15 ($\square$).
Other electronic parameters are given in Fig. 3. Master equation results
appear in dashed-dotted lines.

###  IV.6 Convergence and Computational aspects

Convergence of the path-integral method should be verified with respect to
three numerical parameters: the number of states used to mimic a fermi sea,
$L_{s}$, the time step adopted, $\delta t$, and the memory time accounted for,
$\tau_{c}$. (i) Fermi sea discretization. We have found that excellent
convergence is achieved for relatively “small” fermi reservoirs, taking into
account $L_{s}>20$ states for each reservoir. In our simulations we
practically adopted $L_{s}=30$ for each Fermi bath. (ii) Time-step
discretization. The first criteria in selecting the value of the time step
$\delta t$ is that dynamical features of the isolated vibrational mode should
be observed. Using $\omega_{0}=0.2$, the period of the bath-free Rabbi
oscillation is $2\pi/(\omega_{0})\sim 30$, thus a time step of $\delta t\sim
1$ can capture the details of the TLS oscillation. This consideration serves
as an “upper bound” criteria. The second consideration connects to the time
discretization error, originates from the approximate splitting of the total
time evolution operator into a product of terms, see Eq. (5). For the
particular Trotter decomposition employed, the leading error grows with
$\delta t^{3}\times\left([H_{S},[V,H_{S}]]/12+[V,[V,H_{S}]]/24\right)$ Tannor
(2007) where $V=V_{SB}+V_{SF}+H_{B}+H_{F}$. The decomposition is exact when
the coupling of the subsystem to the reservoirs is weak and the time-step is
small, $\delta t\rightarrow 0$. For large coupling one should take a
sufficiently small time-step in order to avoid significant error buildup. In
the preset work, the dimensionless coupling to the fermi sea
$\pi\rho\kappa\lambda_{l}^{*}\lambda_{r}$ is typically maintained lower that
0.3; the dimensionless coupling to the boson bath is taken as $K_{d}=0.1$. The
value of $\delta t=0.6-1.2$ is thus sufficiently small for our simulations.
(iii) Memory error. Our approach assumes that bath correlations exponentially
decay resulting from the finite temperature and the nonequilibrium condition.
Based on this assumption, the total influence functional was truncated to
include only a finite number of time steps $N_{s}$, where
$\tau_{c}=N_{s}\delta t$. The total IF is retrieved by taking the limit
$N_{s}\rightarrow N$, ($N=t/\delta t$). Our simulations were performed for
$N_{s}=3...9$, covering memory time up to $\tau_{c}=N_{s}\delta t\sim 10$. The
results displayed converged for $N_{s}\sim 7-9$ for $\delta t=1$.

Computational efforts can be partitioned into two parts: In the initialization
step the (time invariant) IFs are computed. The size of the fermionic IF is
$d^{2N_{s}}$, where $d$ is the dimensionality of the subsystem (two in our
simulations). The power of two at the exponent results from the forward and
backward time evolution operators in the path-integral expression. This
initialization effort thus scales exponentially with the memory size accounted
for. The preparation of the bosonic IF is more efficient if the FV IF is used
Makri and Makarov (1995a); *QUAPI2. In the second, time evolution, stage, we
iteratively apply the linear map (27) or (34), a multiplication of two objects
of length $d^{2N_{s}}$. This operation linearly scales in time.

We now comment on the simulation time of a convergence analysis as presented
in Fig. 10, covering three different time steps and $N_{s}=3,...,9$. The
MATLAB implementation of the computational algorithm took advantage of the
MATLAB built-in multi-threaded parallel features and utilized 100$\%$ of all
available CPU cores on a node. When executed on one cluster node with two
quad-core 2.2GHz AMD Opteron cpus and 16GB memory, convergence analysis for of
the full voltage scan took about 4x24 hours and 250MB of memory. Computations
performed on the GPC supercomputer at the SciNet HPC Consortium Loken and et
al. (2010) were three times faster. Computational time scales linearly with
the simulated time $t$. For a fixed $N_{s}$ value, the computational effort
does not depend on the system temperature and other parameters employed.

##  V Summary

We have developed an iterative numerically-exact path-integral scheme that can
follow the dynamics, to the steady-state limit, of subsystems coupled to
multiple bosonic and fermionic reservoirs in an out-of-equilibrium initial
state. The method is based on the truncation of time correlations in the
influence functional, beyond the memory time dictated by temperature and
chemical biases. It combines two techniques: the QUAPI method Makri and
Makarov (1995a); *QUAPI2, for treating the dynamics of subsystems coupled to
harmonic baths, and the INFPI approach Segal  _et al._ (2010); *IF2, useful
for following the evolution of a subsystem when interacting with fermionic
baths.

The method is stable, efficient, and flexible, and it allows one to achieve
transient and steady-state data for both the reduced density matrix of the
subsystem and expectation values of operators, such as the charge current and
energy current. The method can be viewed as an extension of QUAPI, to
incorporate fermions in the dynamics. It could be further expanded to include
time-dependent Hamiltonians, e.g., pulsed fields.

To demonstrate the method usability in the field of molecular conduction, we
have applied the general scheme, and studied vibrational dynamics in a
molecular rectifier setup, where vibrational equilibration with an additional
phonon bath is allowed. Our main conclusions in this study are the following:
(i) The vibrational instability effect disappears once the vibrational mode is
weakly coupled ($K_{d}\sim 0.01$) to an additional phonon bath that can
dissipate the excess energy. (ii) When $K_{d}\sim 0.1$, the vibrational mode
is equilibrated with the secondary phonon bath. (iii) The charge current does
not testify on vibrational heating and instability. While we have performed
those simulations using a truncated vibrational mode, a TLS, representing an
anharmonic mode, we argue that the main characteristics of the vibrational
instability effect remain intact when the selected mode is made harmonic
Simine and Segal (2012).

Our simulations indicate that Master equation methods can excellently
reproduce exact results at weak coupling, in the markovian limit. More
significantly, Master equation tools can be used beyond the weak coupling
limit ($\kappa\sim 0.3$), if only a qualitative understanding of the junction
behavior is inquired. One should note that our Master equation technique
treats the D and A coupling to the metals exactly. It is perturbative only in
the interaction of the vibrational mode to the electrons, and to other phonon
degrees of freedom. In the case where tunneling transmission competes with
phonon-assisted transport, only path-integral simulations were provided, as
more involved Master equation methodologies should be developed in this case.

Our future objectives are twofold: (i) to improve the time-evolution
algorithm, and (ii) to employ the method for the study of other problems in
molecular electronics and phononics. By improving the methodology, we would
like to extend the usability of our method to difficult parameter regimes
(strong coupling), e.g., by developing an equation-of-motion for the memory
function Golosov  _et al._ (1999); Cohen and Rabani (2011). This will also
allow us to simulate more feasibly the dynamics of an $n-$level subsystem.
Another related objective is the study of heat current characteristics in the
spin-boson molecular junction Segal and Nitzan (2005a); *SegalR2. The single-
bath spin-boson model displays a rich dynamics with a complex phase diagram.
Similarly, we expect that the nonequilibrium version, with two harmonic baths
of different temperatures coupled to the TLS, will show complex behavior for
its heat current- temperature characteristics. Recent results, obtained using
an extension of the noninteracting blip approximation to the nonequilibrium
regime Nicolin and Segal (2011), demonstrate rich behavior. Other problems
that could be directed with our method include plexcitonics systems, as the
coupling between surface plasmons and molecular excitons should be treated
beyond the perturbative regime Manjavacas  _et al._ (2011). Finally, we have
discussed the calculation of reduced density matrix and currents in the path-
integral framework. It is of interest to generalize these expressions and gain
higher order cumulants, for the study of current, noise, and fluctuation
relations in many-body out-of-equilibrium systems.

###### Acknowledgements.

DS acknowledges support from an NSERC discovery grant. The work of LS was
supported by an Early Research Award of DS. Computations were performed on the
GPC supercomputer at the SciNet HPC Consortium Loken and et al. (2010). SciNet
is funded by: the Canada Foundation for Innovation under the auspices of
Compute Canada; the Government of Ontario; Ontario Research Fund - Research
Excellence; and the University of Toronto.

## Appendix A: Time-discrete Feynman-Vernon Influence functional

With the discretization of the path, the influence functional takes the form
(17). The coefficients $\eta_{k,k^{\prime}}$ were given in Makri and Makarov
(1995a); *QUAPI2 and we include them here for the completeness of our
presentation. The expressions are given here for the case of a single boson
bath with the initial temperature $1/\beta_{ph}$ and the spectral function
$J_{ph}(\omega)=\pi\sum_{p}(\xi^{B}_{p})^{2}\delta(\omega-\omega_{p})$,
$J_{ph}(\omega)=J_{ph}(-\omega)$,

$\displaystyle\eta_{k,k^{\prime}}$ $\displaystyle=$
$\displaystyle\frac{2}{\pi}\int_{-\infty}^{\infty}d\omega\frac{J_{ph}(\omega)}{\omega^{2}}\frac{\exp(\beta_{ph}\omega/2)}{\sinh(\beta_{ph}\omega/2)}\sin^{2}(\omega\delta
t/2)e^{-i\omega\delta t(k-k^{\prime})},\,\,\,\,\,0<k^{\prime}<k<N$
$\displaystyle\eta_{k,k}$ $\displaystyle=$
$\displaystyle\frac{1}{2\pi}\int_{-\infty}^{\infty}d\omega\frac{J_{ph}(\omega)}{\omega^{2}}\frac{\exp(\beta_{ph}\omega/2)}{\sinh(\beta_{ph}\omega/2)}\left(1-e^{-i\omega\delta
t}\right),\,\,\,\,0<k<N$ $\displaystyle\eta_{k,0}$ $\displaystyle=$
$\displaystyle\frac{2}{\pi}\int_{-\infty}^{\infty}d\omega\frac{J_{ph}(\omega)}{\omega^{2}}\frac{\exp(\beta_{ph}\omega/2)}{\sinh(\beta_{ph}\omega/2)}\sin(\omega\delta
t/4)\sin(\omega\delta t/2)e^{-i\omega(k\delta t-\delta
t/4)},\,\,\,\,\,\,\,0<k<N$ $\displaystyle\eta_{N,k^{\prime}}$ $\displaystyle=$
$\displaystyle\frac{2}{\pi}\int_{-\infty}^{\infty}d\omega\frac{J_{ph}(\omega)}{\omega^{2}}\frac{\exp(\beta_{ph}\omega/2)}{\sinh(\beta_{ph}\omega/2)}\sin(\omega\delta
t/4)\sin(\omega\delta t/2)e^{-i\omega(N\delta t-k^{\prime}\delta t-\delta
t/4)},\,\,\,\,\,\,\,0<k^{\prime}<N$ $\displaystyle\eta_{N,0}$ $\displaystyle=$
$\displaystyle\frac{2}{\pi}\int_{-\infty}^{\infty}d\omega\frac{J_{ph}(\omega)}{\omega^{2}}\frac{\exp(\beta_{ph}\omega/2)}{\sinh(\beta_{ph}\omega/2)}\sin^{2}(\omega\delta
t/4)e^{-i\omega(N\delta t-\delta t/2)}$ $\displaystyle\eta_{0,0}$
$\displaystyle=$
$\displaystyle\eta_{N,N}=\frac{1}{2\pi}\int_{-\infty}^{\infty}d\omega\frac{J_{ph}(\omega)}{\omega^{2}}\frac{\exp(\beta_{ph}\omega/2)}{\sinh(\beta_{ph}\omega/2)}\left(1-e^{-i\omega\delta
t/2}\right)$ (A1)

## References

  * Bonca and Trugmann (1995) J. Bonca and S. Trugmann, Phys. Rev. Lett. 75, 2566 (1995).
  * Ness  _et al._ (2001) H. Ness, S. Shevlin, and A. Fisher, Phys. Rev. B 63, 125422 (2001).
  * Cízek  _et al._ (2004) M. Cízek, M. Thoss, and W. Domcke, Phys. Rev. B 70, 125406 (2004).
  * Galperin  _et al._ (2007a) M. Galperin, M. A. Ratner, and A. Nitzan, J. Phys.: Condens. Matter 19, 103201 (2007a).
  * Mitra  _et al._ (2004) A. Mitra, I. Aleiner, and A. J. Millis, Phys. Rev. B 69, 245302 (2004).
  * Galperin  _et al._ (2007b) M. Galperin, A. Nitzan, and M. A. Ratner, Phys. Rev. B 75, 155312 (2007b).
  * Galperin  _et al._ (2009) M. Galperin, M. A. Ratner, and A. Nitzan, J. Chem. Phys 130, 144109 (2009).
  * Fransson and Galperin (2010) J. Fransson and M. Galperin, Phys. Rev. B 81, 075311 (2010).
  * Segal and Nitzan (2002) D. Segal and A. Nitzan, J. Chem. Phys. 117, 3915 (2002).
  * Leijnse and Wegewijs (2008) M. Leijnse and M. R. Wegewijs, Phys. Rev. B 78, 235424 (2008).
  * Hartle and Thoss (2011a) R. Hartle and M. Thoss, Phys. Rev. B 83, 125419 (2011a).
  * Hartle and Thoss (2011b) R. Hartle and M. Thoss, Phys. Rev. B 83, 115414 (2011b).
  * Volkovich  _et al._ (2010) R. Volkovich, R. Hartle, and M. Thoss, Phys. Chem. Chem. Phys. 32, 14333 (2010).
  * Swenson  _et al._ (2011) D. W. H. Swenson, T. Levy, G. Cohen, E. Rabani, and W. H. Miller, J. Chem. Phys. 134, 164103 (2011).
  * Swenson  _et al._ (2012) D. W. H. Swenson, G. Cohen, and E. Rabani, Molecular Physics 110, 743 (2012).
  * Makri and Makarov (1995a) N. Makri and D. E. Makarov, J. Chem. Phys. 102, 4600 (1995a).
  * Makri and Makarov (1995b) N. Makri and D. E. Makarov, J. Chem. Phys. 102, 4611 (1995b).
  * Segal  _et al._ (2010) D. Segal, A. J. Millis, and D. R. Reichman, Phys. Rev. B 82, 205323 (2010).
  * Segal  _et al._ (2011) D. Segal, A. J. Millis, and D. R. Reichman, Phys. Chem. Chem. Phys. 13, 14378 (2011).
  * Bedkihal and Segal (2012) S. Bedkihal and D. Segal, Phys. Rev. B 85, 155324 (2012).
  * Kulkarni  _et al._ (2012) M. Kulkarni, K. L. Tiwari, and D. Segal, Phys. Rev. B 86, 155424 (2012).
  * Kulkarni  _et al._ (2013) M. Kulkarni, K. L. Tiwari, and D. Segal, New J. Phys. 15, 013014 (2013).
  * Muhlbacher and Rabani (2008) L. Muhlbacher and E. Rabani, Phys. Rev. Lett. 100, 176403 (2008).
  * Wang  _et al._ (2011) H. Wang, I. Pshenichnyuk, R. Hartle, and M. Thoss, J. Chem. Phys. 135, 244506 (2011).
  * Weiss  _et al._ (2008) S. Weiss, J. Eckel, M. Thorwart, and R. Egger, Phys. Rev. B 77, 195316 (2008).
  * Eckel  _et al._ (2010) J. Eckel, F. Heidrich-Meisner, S. G. Jakobs, M. Thorwart, M. Pletyukhov, and R. Egger, New J. Phys. 12, 043042 (2010).
  * Hutzen  _et al._ (2012) R. Hutzen, S. Weiss, M. Thorwart, and R. Egger, Phys. Rev. B 85, 121408 (2012).
  * Aviram and Ratner (1974) A. Aviram and M. A. Ratner, Chem. Phys. Lett. 29, 277 (1974).
  * Lu  _et al._ (2011) J.-T. Lu, P. Hedegard, and M. Brandbyge, Phys. Rev. Lett. 107, 046801 (2011).
  * Simine and Segal (2012) L. Simine and D. Segal, Phys. Chem. Chem. Phys. 14, 13820 (2012).
  * Kondo (1964) J. Kondo, Prog. Theor. Phys. 32, 37 (1964).
  * Mitra and Millis (2005) A. Mitra and A. J. Millis, Phys. Rev. B 72, 121102(R) (2005).
  * Mitra and Millis (2007) A. Mitra and A. J. Millis, Phys. Rev. B 76, 085342 (2007).
  * Segal  _et al._ (2007) D. Segal, D. R. Reichman, and A. J. Millis, Phys. Rev. B 76, 195316 (2007).
  * Lutchyn  _et al._ (2008) R. M. Lutchyn, L. Cywinski, C. P. Nave, and S. D. Sarma, Phys. Rev. B 78, 024508 (2008).
  * Feynman and Hibbs (1965) R. P. Feynman and A. R. Hibbs,  _Quantum Mechanics and Path Integrals_ (McGraw-Hill, New-York, 1965).
  * Klich (2003) I. Klich, “in quantum noise in mesoscopic systems,” (Kluwer, 2003).
  * Yu.  _et al._ (2004) L. H. Yu., Z. K. Keane, J. W. Ciszek, L. Cheng, M. P. Stewart, J. M. Tour, and D. Natelson, Phys. Rev. Lett. 93, 266802 (2004).
  * Djukic  _et al._ (2005) D. Djukic, K. S. Thygesen, C. Untiedt, R. H. M. Smit, K. Jacobsen, and J. M. van Ruitenbeek, Phys. Rev. B 71, 161402 (2005).
  * Kumar  _et al._ (2012) M. Kumar, R. Avriller, A. L. Yeyati, and J. M. van Ruitenbeek, Phys. Rev. Lett. 108, 146602 (2012).
  * Pasupathy  _et al._ (2005) A. N. Pasupathy, J. Park, C. Chang, A. V. Soldatov, R. C. B. S. Lebedkin, J. E. Grose, L. A. K. Donev, J. P. Sethna, D. C. Ralph, and P. L. McEuen, Nano Lett. 5, 203 (2005).
  * Ioffe  _et al._ (2008) Z. Ioffe, T. Shamai, A. Ophir, G. Noy, I. Yutsis, K. Kfir, O. Cheshnovsky, and Y. Selzer, Nature Nanotech. 3, 727 (2008).
  * Ward  _et al._ (2011) D. R. Ward, D. A. Corley, J. M. Tour, and D. Natelson, Nature Nanotech. 6, 33 (2011).
  * Huang  _et al._ (2007) Z. Huang, F. Chen, R. Dagosta, P. A. Bennett, M. D. Ventra, and N. Tao, Nature Nanotech. 2, 698 (2007).
  * Entin-Wohlman  _et al._ (2010) O. Entin-Wohlman, Y. Imry, and A. Aharony, Phys. Rev. B 82, 115314 (2010).
  * Jiang  _et al._ (2012) J.-H. Jiang, O. Entin-Wohlman, and Y. Imry, Phys. Rev. B 85, 075412 (2012).
  * Entin-Wohlman and Aharony (2012) O. Entin-Wohlman and A. Aharony, Phys. Rev. B 85, 085401 (2012).
  * Mahan (2000) G. D. Mahan,  _Many-particle physics_ (Plenum press, New York, 2000).
  * Roulet  _et al._ (1969) B. Roulet, J. Gavoret, and P. Nozieres, Phys. Rev. 178, 1072 (1969).
  * Mii  _et al._ (2003) T. Mii, S. Tikhodeev, and H. Ueba, Phys. Rev. B 68, 205406 (2003).
  * Galperin  _et al._ (2004a) M. Galperin, M. Ratner, and A. Nitzan, J. Chem. Phys. 121, 11965 (2004a).
  * Galperin  _et al._ (2004b) M. Galperin, M. Ratner, and A. Nitzan, Nano. Lett. 4, 1605 (2004b).
  * Tannor (2007) D. J. Tannor,  _Intriduction to Quantum Mechanics: A time dependent perspective_ (University Science Books, 2007).
  * Loken and et al. (2010) C. Loken and et al., J. Phys.: Conf. Ser. 256, 012026 (2010).
  * Golosov  _et al._ (1999) A. A. Golosov, R. A. Friesner, and P. Pechukas, J. Chem. Phys 110, 138 (1999).
  * Cohen and Rabani (2011) G. Cohen and E. Rabani, Phys. Rev. B 84, 075150 (2011).
  * Segal and Nitzan (2005a) D. Segal and A. Nitzan, Phys. Rev. Lett. 94, 034301 (2005a).
  * Segal and Nitzan (2005b) D. Segal and A. Nitzan, J. Chem. Phys. 122, 194704 (2005b).
  * Nicolin and Segal (2011) L. Nicolin and D. Segal, J. Chem. Phys. 135, 164106 (2011).
  * Manjavacas  _et al._ (2011) A. Manjavacas, F. J. G. de Abajo, and P. Nordlander, Nano Lett. 11, 2318 (2011).