††thanks<EMAIL_ADDRESS>
# Top-Ranked Cycle Flux Network Analysis of Molecular Photocells

Nikhil Gupt  Indian Institute of Technology Kanpur, Kanpur, Uttar Pradesh
208016, India  Shuvadip Ghosh  Indian Institute of Technology Kanpur, Kanpur,
Uttar Pradesh 208016, India  Arnab Ghosh  Indian Institute of Technology
Kanpur, Kanpur, Uttar Pradesh 208016, India

###### Abstract

We introduce a top-ranked cycle flux ranking scheme of network analysis to
assess the performance of molecular junction solar cells. By mapping the
Lindblad master equation to the quantum-transition network, we propose a
microscopic Hamiltonian description underpinning the rate equations commonly
used to characterize molecular photocells. Our approach elucidates the
paramount significance of edge flux and unveils two pertinent electron
transfer pathways that play equally important roles in robust photocurrent
generation. Furthermore, we demonstrate that non-radiative loss processes
impede the maximum power efficiency of photocells, which may otherwise be
above the Curzon-Ahlborn limit. These findings shed light on the intricate
functionalities that govern molecular photovoltaics and offer a comprehensive
approach to address them in a systematic way.

##  I Introduction

Molecular junctions consisting of a single donor and acceptor, placed between
two external leads is an active field of research Nitzan (2006); Bittner
(2010); Aradhya and Venkataraman (2013); Su  _et al._ (2016) that combines the
fundamental aspects of quantum transport Scully (2010); Scully  _et al._
(2011); Thoss and Evers (2018) and their possible practical implementations
Chen (2005); Pop (2010); Killoran  _et al._ (2015); Arp  _et al._ (2016). One
such application of molecular junctions, among many others, is in photovoltaic
(PV) cells, where solar energy of incident photons is converted into electric
power Deibel and Dyakonov (2010); Nicholson and Castro (2010); Kamat (2013).
The interplay between the heat current mediated by the temperature difference
of the solar radiation and the PV cell at ambient temperature, together with
the charge current arising due to the bias voltage across the electrodes,
drive the system towards strong non-equilibrium steady states Altaner  _et
al._ (2012); Simine and Segal (2012); Agarwalla  _et al._ (2015); Ajisaka  _et
al._ (2015). A great deal is therefore focused on theoretical modeling to
optimize the performance of PV cells Creatore  _et al._ (2013); Fruchtman  _et
al._ (2016) and to explore the underlying transport mechanisms that could
facilitate the realization of sophisticated on-chip complex quantum thermal
devices Giazotto and Martínez-Pérez (2012); Joulain  _et al._ (2016); Whitney
_et al._ (2018); Fong  _et al._ (2019); Pekola and Karimi (2021).

Though, the master equation in Lindblad form is a conventional theoretical
technique for the description of non-equilibrium open quantum systems Breuer
and Petruccione (2007), particularly in quantum optics Agarwal (2012) and
quantum thermodynamics communities Sinha  _et al._ (2013); Gupt  _et al._
(2021); Gelbwaser-Klimovsky  _et al._ (2015); Kosloff (2013); Deffner and
Campbell (2019); Dong  _et al._ (2021); Gupt  _et al._ (2022); Ghosh  _et al._
(2022); Kumar  _et al._ (2023) where spins and atomic degrees are involved,
its extension towards molecular systems is relatively a novel area of
exploration Aviram and Ratner (1974); Harbola  _et al._ (2006); Galperin  _et
al._ (2007); Wang  _et al._ (2020); Sarkar and Dubi (2020); Oz  _et al._
(2023). As an alternative, in recent times, a state space representation in
the framework of network theory by Nitzan et. al. Einax  _et al._ (2011,
2013); Einax and Nitzan (2014, 2016); Craven  _et al._ (2018) has become a
popular method to study the non-equilibrium charge transport behavior of
molecular photovoltaics Einax  _et al._ (2011, 2013); Einax and Nitzan (2014,
2016) and thermoelectric devices Whitney  _et al._ (2018); Craven  _et al._
(2018); Mazza  _et al._ (2014). In this context, the latest finding by Wang et
al. Wang  _et al._ (2022), based on the works of Zurek Zurek (1981), and Cao
Wu and Cao (2013), is worth mentioning. They have shown that the quantum
Lindblad master equation can be cast into a Pauli master equation without loss
of any generality. This facilitates the representation of dissipative quantum
dynamics as a weighted network with nodes and edges, where nodes (vertices)
denote quantum states and edges denote the non-equilibrium transition from one
quantum state to another with nonzero flux rates Chernyak and Sinitsyn (2008);
Rahav  _et al._ (2008). Thus, optimization of the performance of multi-
component quantum thermal devices reduces to identifying the major working
cycles amongst various possible pathways as is the case with photovoltaics
with a multitude of electron transfer channels. Although the concept of cycle
flux was well developed in algebraic graph theory Hill and Chen (1975); Kohler
and Vollmerhaus (1980); Schnakenberg (1976); Tutte (2001); Balakrishnan and
Ranganathan (2012); Polettini and Esposito (2014); Owen  _et al._ (2020) from
the early works of Hill Hill and Chen (1975), Kohler and Vollmerhaus Kohler
and Vollmerhaus (1980), as well as Schnakenberg Schnakenberg (1976), the
recent developments by Wang et. al. have provided an efficient cycle flux
ranking scheme to fully comprehend the intricate functionality of complex
quantum systems, with particular emphasis on spin-Seebeck effect within the
linear response regime Wang  _et al._ (2022) .

In this paper, we extend this idea to molecular systems and establish the
equivalence between the dynamical formulation of the Lindblad master equation
and the state space representation of molecular photocells pioneered by Nitzan
et. al. Einax  _et al._ (2011, 2013); Einax and Nitzan (2014, 2016). We point
out however the state space method is quite effective in computing the steady-
state currents in PV cells, it falls short to unravel the underlying working
mechanism of the photovoltaic devices. On the contrary, top-ranked cycle flux
analysis provides a natural and alternative gateway to capture the underlying
features, which is otherwise challenging due to multiple electron transfer
pathways in molecular junctions. Thus, our present findings demonstrate that
the cycle flux ranking scheme could go beyond its standard applications of
near-equilibrium situations and could serve as a potential candidate for
decoding the fundamental working principle of complex molecular systems even
far from equilibrium scenarios.

The work is organized as follows: In section II, we introduce the basic model
of the PV cell and derive the open quantum dynamics of Lindbladian form which
is shown to be equivalent to the Pauli master equations employed by the Nitzan
et. al. Einax  _et al._ (2011) within a state space formulation. Next, we
elaborate on the fundamental principles of the cycle flux ranking scheme in
the context of our present model in Sec. III, and summarize important outcomes
and findings of our analysis in Sec. IV. Finally, we conclude in Sec. V.

##  II Model and Dynamics

The basic model of a PV cell is comprised of two “effective” sites,
representing the donor (D) and acceptor (A) molecules placed between two
metallic leads (L and R), as depicted in Fig. 1. We consider each site as a
two-state system with the ground ($|n_{\rm s1}\rangle$) and excited ($|n_{\rm
s2}\rangle$) states, corresponding to the highest occupied molecular orbital
(HOMO) and lowest unoccupied molecular orbital (LUMO) of two molecules having
energies $\varepsilon_{\rm s1}$ and $\varepsilon_{\rm s2}$ (${\rm s=D,A}$)
respectively. The general form of the electronic part of the Hamiltonian for
the donor-acceptor system reads as

$H_{el}=\sum_{\rm s=D,A}\sum_{\rm i=1,2}\varepsilon_{\rm si}n_{\rm
si}+\sum_{\rm s=D,A}U_{\rm s}n_{\rm s1}n_{\rm s2}.$ (1)

Here, $U_{\rm D}$ $(U_{\rm A})$ stands for the positive Coulomb repulsion
energy if two electrons are present in the donor (acceptor) sites and $n_{\rm
si}=c^{\dagger}_{\rm si}c_{\rm si}$, ${\rm i}=\\{1,2\\}$ represents the number
operator corresponding to the state $|n_{\rm si}\rangle$, satisfying fermionic
anti-commutation relation Ghosh  _et al._ (2012) $\\{c_{\rm
si},c^{\dagger}_{\rm
s^{\prime}j}\\}=\delta_{\rm{ss^{\prime}}}\delta_{\rm{ij}}$. The most important
energy scales of the problem are notably the energy gap between the donor
levels ($\Delta E=\varepsilon_{\rm D2}-\varepsilon_{\rm D1}$) and the donor-
acceptor levels ($\Delta\varepsilon=\varepsilon_{\rm D2}-\varepsilon_{\rm
A2}$), while the acceptor’s ground state energy has no influence on the
overall cell operation as one expects Einax  _et al._ (2011); Sarkar and Dubi
(2020). This can be attributed to the following two facts (i) solar radiation
triggers the radiative transition of electrons from the donor ground to the
excited state; and (ii) electron transfer between the excited states of the
donor to the acceptor could in principle be governed by the cell’s vibrational
degrees of freedom, such as thermal phonons modes. To make a further realistic
assumption about our model in view of the prior literature Einax  _et al._
(2011, 2013); Einax and Nitzan (2014); Sarkar and Dubi (2020), we restrict the
donor to be either in the ground or excited state (i.e., $n_{\rm D1}n_{\rm
D2}=0$), whereas the acceptor can be either in the singly or doubly occupied
state. As the ground state of the acceptor is always occupied (i.e. $n_{\rm
A1}=1$), it can be either in state $|n_{\rm A1}n_{\rm A2}\rangle$=$|10\rangle$
or in state $|n_{\rm A1}n_{\rm A2}\rangle$=$|11\rangle$, while, the donor
state $|n_{\rm D1}n_{\rm D2}\rangle$ can exist in any one of the three
possible configurations $\\{|00\rangle,|10\rangle,|01\rangle\\}$. This allows
us to directly compare our results with the state space model of Nitzan et.
al. Einax  _et al._ (2011), characterized by a single Coulomb repulsion
parameter arising out of $U_{\rm{A}}$ in Eq. (1).

Figure 1: Schematic diagram of a molecular photovoltaic solar cell. The system
consists of the donor (D) and the acceptor (A) molecules placed between two
free-electron reservoirs. The left reservoir ($\rm L$) is exclusively coupled
to the ground state of the donor and the right reservoir ($\rm R$) is coupled
only to the excited state of the acceptor. Photon-induced transition is
indicated by a wiggly line and phonon-induced processes are depicted by broken
lines.

The Hamiltonian of the free-electron reservoirs Ghosh  _et al._ (2022) to
which the donor and acceptor molecules are coupled are respectively given by

$H_{\rm L}=\sum_{l}(\epsilon_{l}-\mu_{\rm
L})d^{\dagger}_{l}d_{l},\quad\text{and}\quad H_{\rm
R}=\sum_{r}(\epsilon_{r}-\mu_{\rm R})d^{\dagger}_{r}d_{r},$ (2)

where $d_{l}(d_{r}$) and $d^{\dagger}_{l}(d^{\dagger}_{r}$) are the electron
annihilation and creation operators for left (L) and right (R) electrodes and
$\mu_{\alpha}$ is the corresponding chemical potential of $\alpha$-th
($\alpha={\rm L,R}$) reservoirs. In other words, the bias voltage between two
metallic electrodes is given by $U=(\mu_{\rm L}-\mu_{\rm R})/|e|$, where $e$
is the charge of the electron. Tunneling Hamiltonian between the electrodes
and the molecule is chosen to be of the form of Simine and Segal (2012); Ghosh
_et al._ (2022)

$H_{I}=\sum_{l}\hbar g_{l}(d^{\dagger}_{l}c_{\rm D1}+c^{\dagger}_{\rm
D1}d_{l})+\sum_{r}\hbar g_{r}(d^{\dagger}_{r}c_{\rm A2}+c^{\dagger}_{\rm
A2}d_{r}),$ (3)

where $g_{l(r)}$ is the respective coupling strength. The first term indicates
that the left electrode swaps electrons solely with the ground state of the
donor which gets triggered by the solar radiation to its excited state. The
second term denotes that the right electrode can only exchange electrons with
the excited state of the acceptor, as the ground state of the acceptor
molecule is already occupied. Finally, electron transfer processes within the
system are governed by both photons ($H_{\rm pht}$) and phonon ($H_{\rm phn}$)
baths Sarkar and Dubi (2020) with $H_{\rm
pht}=\sum_{k}\epsilon_{k}a^{\dagger}_{k}a_{k}$, and $H_{\rm
phn}=\sum_{q}\epsilon^{\prime}_{q}b^{\dagger}_{q}b_{q}$, where
$a^{\dagger}_{k}$ ($b^{\dagger}_{q}$) and $a_{k}$ ($b_{q}$) are bosonic
creation and annihilation operators for the $k$-th ($q$-th) bath mode with
energy $\epsilon_{k}$ ($\epsilon^{\prime}_{q}$) respectively. The interactions
between the photon and phonon baths with the molecule are taken as Simine and
Segal (2012); Ajisaka  _et al._ (2015); Sarkar and Dubi (2020)

$H^{\rm pht}_{I}=\sum_{k}\hbar g^{\rm D}_{k}(a^{\dagger}_{k}c^{\dagger}_{\rm
D1}c_{\rm D2}+c^{\dagger}_{\rm D2}c_{\rm D1}a_{k}),$ (4)

and

$\displaystyle H^{\rm phn}_{I}$ $\displaystyle=$ $\displaystyle\sum_{q}\hbar
g^{\rm D}_{q}(b^{\dagger}_{q}c^{\dagger}_{\rm D1}c_{\rm D2}+c^{\dagger}_{\rm
D2}c_{\rm D1}b_{q})$ (5) $\displaystyle+$ $\displaystyle\sum_{q}\hbar g^{\rm
DA}_{q}(b^{\dagger}_{q}c^{\dagger}_{\rm A2}c_{\rm D2}+c^{\dagger}_{\rm
D2}c_{\rm A2}b_{q}),$

respectively, where $g^{\rm D}_{k}$ and $g^{\rm D}_{q}(g^{\rm DA}_{q})$ are
the corresponding coupling constants. From Eqs. (4) and (5), it is clear that
the transition between the ground to the excited state at the donor site is
governed by both photon (radiative process) and phonon (nonradiative process)
modes, while the electron transfer between the donor and the acceptor, is
solely driven by vibrational thermal phonon modes at the ambient temperature.

Consequently, the time evolution of the system dynamics in the interaction
picture under the Born-Markov approximation is described by the quantum master
equation Breuer and Petruccione (2007); Kosloff (2013); Gupt  _et al._ (2022);
Ghosh  _et al._ (2022) (Appendix A)

$\displaystyle\frac{d\rho}{dt}=\mathcal{L}_{\rm L}[\rho]+\mathcal{L}_{\rm
R}[\rho]+\mathcal{L}^{\rm D}_{\rm pht}[\rho]+\mathcal{L}^{\rm D}_{\rm
phn}[\rho]+\mathcal{L}^{\rm DA}_{\rm phn}[\rho],$ (6)

where $\rho$ is the reduced density matrix of the system and
$\mathcal{L}_{\nu}$ is the Lindblad superoperator describing the effect of
dissipation induced by the $\nu$-th thermal bath. Since the Hamiltonian in Eq.
(1) is diagonal in the number state basis for the donor and acceptor
molecules, the reduced density matrix $\rho$ of the above Lindblad master
equation effectively decouples the diagonal and off-diagonal matrix elements
in the eigenbasis of $H_{el}$ Ghosh  _et al._ (2022). This, in turn, allows
for a closed-form equation of motion for the occupation probabilities or the
population of the $\mathbb{i}$-th site as
$P_{\mathbb{i}}=\langle\mathbb{i}|\rho|\mathbb{i}\rangle$, w.r.t the various
system eigenstates $\\{|\mathbb{i}\rangle\\}$. As we mentioned before, the
acceptor can only be in two states $\\{|10\rangle,|11\rangle\\}$, while the
donor can be any of the three possible configurations
$\\{|00\rangle,|10\rangle,|01\rangle\\}$, so there are six possible
eigenstates for the overall system. We label them as follows:
$|\mathbb{0}\rangle=|0010\rangle$, $|\mathbbm{1}\rangle=|1010\rangle$,
$|\mathbb{2}\rangle=|0110\rangle$, $|\mathbb{3}\rangle=|0011\rangle$,
$|\mathbb{4}\rangle=|1011\rangle$ and $|\mathbb{5}\rangle=|0111\rangle$. The
corresponding eigen-energies $\varepsilon_{\mathbb{j}}$
($\mathbb{j}=\mathbb{0,1,2,..,5}$) for the states are given by:
$\varepsilon_{\mathbb{0}}=\varepsilon_{\rm A1}$,
$\varepsilon_{\mathbb{1}}=\varepsilon_{\rm D1}+\varepsilon_{\rm A1}$,
$\varepsilon_{\mathbb{2}}=\varepsilon_{\rm D2}+\varepsilon_{\rm A1}$,
$\varepsilon_{\mathbb{3}}=\varepsilon_{\rm A1}+\varepsilon_{\rm A2}+U_{\rm
A}$, $\varepsilon_{\mathbb{4}}=\varepsilon_{\rm D1}+\varepsilon_{\rm
A1}+\varepsilon_{\rm A2}+U_{\rm A}$ and
$\varepsilon_{\mathbb{5}}=\varepsilon_{\rm D2}+\varepsilon_{\rm
A1}+\varepsilon_{\rm A2}+U_{\rm A}$ respectively. As a result, the time
evolution equations for the population are governed by the kinetic equations
(detailed derivation in Appendix A)

$\frac{dP_{\mathbb{0}}}{dt}=(k_{\mathbb{01}}P_{\mathbb{1}}-k_{\mathbb{10}}P_{\mathbb{0}})+(k_{\mathbb{03}}P_{\mathbb{3}}-k_{\mathbb{30}}P_{\mathbb{0}}),$
(7)
$\displaystyle\frac{dP_{\mathbb{1}}}{dt}=(k_{\mathbb{10}}P_{\mathbb{0}}-k_{\mathbb{01}}P_{\mathbb{1}})+(k_{\mathbb{12}}P_{\mathbb{2}}-k_{\mathbb{21}}P_{\mathbb{1}})$
$\displaystyle+(k_{\mathbb{14}}P_{\mathbb{4}}-k_{\mathbb{41}}P_{\mathbb{1}}),$
(8)
$\displaystyle\frac{dP_{\mathbb{2}}}{dt}=(k_{\mathbb{21}}P_{\mathbb{1}}-k_{\mathbb{12}}P_{\mathbb{2}})+(k_{\mathbb{23}}P_{\mathbb{3}}-k_{\mathbb{32}}P_{\mathbb{2}})$
$\displaystyle+(k_{\mathbb{25}}P_{\mathbb{5}}-k_{\mathbb{52}}P_{\mathbb{2}}),$
(9)
$\displaystyle\frac{dP_{\mathbb{3}}}{dt}=(k_{\mathbb{30}}P_{\mathbb{0}}-k_{\mathbb{03}}P_{\mathbb{3}})+(k_{\mathbb{32}}P_{\mathbb{2}}-k_{\mathbb{23}}P_{\mathbb{3}})$
$\displaystyle+(k_{\mathbb{34}}P_{\mathbb{4}}-k_{\mathbb{43}}P_{\mathbb{3}}),$
(10)
$\displaystyle\frac{dP_{\mathbb{4}}}{dt}=(k_{\mathbb{41}}P_{\mathbb{1}}-k_{\mathbb{14}}P_{\mathbb{4}})+(k_{\mathbb{43}}P_{\mathbb{3}}-k_{\mathbb{34}}P_{\mathbb{4}})$
$\displaystyle+(k_{\mathbb{45}}P_{\mathbb{5}}-k_{\mathbb{54}}P_{\mathbb{4}}),$
(11)
$\displaystyle\frac{dP_{\mathbb{5}}}{dt}=(k_{\mathbb{52}}P_{\mathbb{2}}-k_{\mathbb{25}}P_{\mathbb{5}})+(k_{\mathbb{54}}P_{\mathbb{4}}-k_{\mathbb{45}}P_{\mathbb{5}}).$
(12)

This is a “classical” looking Pauli master equation where the transition rates
involved are quantum mechanical in nature. For instance, $k_{\mathbb{ji}}$
($k_{\mathbb{j}\leftarrow\mathbb{i}}$) denotes the rate of transition from
quantum state $|\mathbb{i}\rangle$ to $|\mathbb{j}\rangle$
($\mathbb{i},\mathbb{j}=\mathbb{0,1,2,3,4,5}$ but $\mathbb{i}\neq\mathbb{j}$),
which are given by:

$\displaystyle k_{\mathbb{10}}$ $\displaystyle=$ $\displaystyle
k_{\mathbb{43}}=\gamma_{\rm L}f(\epsilon_{\rm L}),$ (13) $\displaystyle
k_{\mathbb{01}}$ $\displaystyle=$ $\displaystyle k_{\mathbb{34}}=\gamma_{\rm
L}[1-f(\epsilon_{\rm L})],$ (14) $\displaystyle k_{\mathbb{30}}$
$\displaystyle=$ $\displaystyle k_{\mathbb{41}}=k_{\mathbb{52}}=\gamma_{\rm
R}f(\epsilon_{\rm R}),$ (15) $\displaystyle k_{\mathbb{03}}$ $\displaystyle=$
$\displaystyle k_{\mathbb{14}}=k_{\mathbb{25}}=\gamma_{\rm
R}[1-f(\epsilon_{\rm R})],$ (16) $\displaystyle k_{\mathbb{21}}$
$\displaystyle=$ $\displaystyle k_{\mathbb{54}}\equiv k_{\rm r}+k_{\rm nr},$
(17) $\displaystyle=$ $\displaystyle\gamma^{\rm D}_{\rm pht}n(\epsilon_{\rm
r})+\gamma^{\rm D}_{\rm phn}n(\epsilon_{\rm nr}),$ $\displaystyle
k_{\mathbb{12}}$ $\displaystyle=$ $\displaystyle
k_{\mathbb{45}}\equiv\tilde{k}_{\rm r}+\tilde{k}_{\rm nr},$ (18)
$\displaystyle=$ $\displaystyle\gamma^{\rm D}_{\rm pht}[n(\epsilon_{\rm
r})+1]+\gamma^{\rm D}_{\rm phn}[n(\epsilon_{\rm nr})+1],$ $\displaystyle
k_{\mathbb{23}}$ $\displaystyle=$ $\displaystyle\gamma^{\rm DA}_{\rm
phn}n(\epsilon_{\rm DA}),$ (19) $\displaystyle k_{\mathbb{32}}$
$\displaystyle=$ $\displaystyle\gamma^{\rm DA}_{\rm phn}[n(\epsilon_{\rm
DA})+1].$ (20)

The transition rates consist of two terms: (a) the rate coefficient $\gamma$
depends on the coupling strength through the bath spectral function that
characterizes the inverse time scale associated with the corresponding
processes; (b) the second term contains the information about the statistical
properties of the quantum bath (fermionic and bosonic) through its
temperature-dependent auto-correlation functions Ghosh  _et al._ (2012). It
can be classified into two categories, namely, the absorption and relaxation
processes. For a fermionic bath, the absorption or excitation process is
governed by $f(\epsilon_{\alpha})$, while the de-excitation process is
controlled by $1-f(\epsilon_{\alpha})$ factor Gupt  _et al._ (2021); Ghosh
_et al._ (2022). The $f(\epsilon_{\alpha})$ is the Fermi-Dirac distribution
which is given by $f(\epsilon_{\alpha})={1}/(e^{\epsilon_{\alpha}}+1)\equiv
f(\varepsilon_{\mathbb{j}\mathbb{i}},\mu_{\alpha},T_{\alpha})$, where
$\alpha={\rm L},{\rm R}$,
$\epsilon_{\alpha}=(\varepsilon_{\mathbb{j}\mathbb{i}}-\mu_{\alpha})/k_{B}T_{\alpha}$
and $T_{\alpha}$ is the temperature of the $\alpha$-th reservoir. In the case
of bosonic (photon and phonon) bath, the same factors are given by
$n(\epsilon_{\beta})$ and $1+n(\epsilon_{\beta})$ Gupt  _et al._ (2021), with
$n(\epsilon_{\beta})$ as the Bose-Einstein distribution
$n(\epsilon_{\beta})={1}/(e^{\epsilon_{\beta}}-1)$, for $\beta={\rm r},{\rm
nr},{\rm DA}$, where $\epsilon_{\rm r}=(\varepsilon_{\rm D2}-\varepsilon_{\rm
D1})/k_{B}T_{\rm S}$, $\epsilon_{\rm nr}=(\varepsilon_{\rm
D2}-\varepsilon_{\rm D1})/k_{B}T_{\rm ph}$, $\epsilon_{\rm
DA}=(\Delta\varepsilon-U_{\rm A})/k_{B}T_{\rm ph}$ and $T_{\rm S}$ ($T_{\rm
ph}$) stands for the temperature of the photon (phonon) bath respectively.

We emphasize that Eqs. (7)-(12), along with their rate coefficients calculated
from the above microscopic picture, are analogous to the phenomenological rate
equations considered by Nitzan et. al. Einax  _et al._ (2011, 2013); Einax and
Nitzan (2014) in their study of molecular photocells. In terms of the
occupation probabilities of the individual states, empirical definitions for
the electron currents used by Nitzan et. al. Einax  _et al._ (2011), within
the state-space model, can be summarized as

$\displaystyle J_{\rm R}$ $\displaystyle=$
$\displaystyle[k_{\mathbb{30}}P_{\mathbb{0}}-k_{\mathbb{03}}P_{\mathbb{3}}]+[k_{\mathbb{41}}P_{\mathbb{1}}-k_{\mathbb{14}}P_{\mathbb{4}}]$
(21) $\displaystyle+$
$\displaystyle[k_{\mathbb{52}}P_{\mathbb{2}}-k_{\mathbb{25}}P_{\mathbb{5}}],$
$\displaystyle J_{\rm L}$ $\displaystyle=$
$\displaystyle[k_{\mathbb{10}}P_{0}-k_{\mathbb{01}}P_{\mathbb{1}}]+[k_{\mathbb{43}}P_{\mathbb{3}}-k_{\mathbb{34}}P_{\mathbb{4}}],$
(22) $\displaystyle J_{\rm S}$ $\displaystyle=$ $\displaystyle k_{\rm
r}[P_{\mathbb{1}}+P_{\mathbb{4}}]-\tilde{k}_{\rm
r}[P_{\mathbb{2}}+P_{\mathbb{5}}],$ (23) $\displaystyle J_{\rm nr}$
$\displaystyle=$ $\displaystyle k_{\rm
nr}[P_{\mathbb{1}}+P_{\mathbb{4}}]-\tilde{k}_{\rm
nr}[P_{\mathbb{2}}+P_{\mathbb{5}}],$ (24) $\displaystyle J_{\rm DA}$
$\displaystyle=$ $\displaystyle
k_{23}P_{\mathbb{3}}-k_{\mathbb{\mathbb{32}}}P_{\mathbb{2}}.$ (25)

Here $J_{\rm R}$ and $J_{\rm L}$ are the electron currents leaving and
entering the molecular system to and from the electrodes, $J_{\rm S}$ and
$J_{\rm nr}$ are, respectively, the radiative (photon-induced) and
nonradiative (phonon induced losses) electron currents between the ground and
excited states of the donor, and $J_{\rm DA}$ is the average current due to
transfer of electrons (phonon induced) between donor and acceptor molecules.
While the above rate equations can be solved numerically to obtain the steady-
state flux-current, $J_{\rm L}=-J_{\rm R}=-J_{\rm{DA}}=J_{\rm
S}+J_{\rm{nr}}=J$, it does not reveal much about the actual mechanisms of the
transport processes involved. In what follows, we take advantage of the cycle
flux analysis scheme of the algebraic graph theory Hill and Chen (1975);
Kohler and Vollmerhaus (1980); Schnakenberg (1976); Tutte (2001); Balakrishnan
and Ranganathan (2012); Polettini and Esposito (2014); Owen  _et al._ (2020)
to gain a deeper insight into the underlying transport channels of the
molecular solar cells and explore the effect of non-radiative losses on the
thermodynamic efficiency of the photocell.

##  III The Cycle flux Analysis

From Eqs. (7)-(12), we may rewrite the time evolution equation of the
population distribution $P_{\mathbb{i}}$ over the quantum mechanical system
states in the following compact form

$\displaystyle\frac{dP_{\mathbb{i}}}{dt}=\sum^{\mathbb{5}}_{\mathbb{j=0}}[k_{\mathbb{ij}}P_{\mathbb{j}}-k_{\mathbb{ji}}P_{\mathbb{i}}];\quad\mathbb{j}\neq{\mathbb{i}},$
(26)

subject to the condition $\sum_{\mathbb{i=1}}^{\mathbb{5}}P_{\mathbb{i}}=1$.
Thus it completely characterizes the overall system where $k_{\mathbb{ji}}$
($k_{\mathbb{ij}}$) depicts the open quantum system under non-equilibrium
conditions with forward (backward) transition rates listed in Eqs. (13)-(20).
With the help of higher-order kinetic expansion, dissipative dynamics in
presence of quantum coherence can be even mapped onto kinetic network Wu and
Cao (2013). Therefore, without loss of any generality, and following the Ref.
Wang  _et al._ (2022), the above dissipative quantum dynamics can be
effectively represented by a network or graph, where nodes or vertices are
assigned to each quantum state, and edges or lines relate the allowed
transitions between them. In the present case, the basic graph ($\mathbb{G}$)
depicted in Fig. 2 provides a useful visualization of the system, with each
vertex denoting a quantum state and its associated probability
$P_{\mathbb{i}}$, and each edge representing a pair of transition
probabilities $k_{\mathbb{ji}}$ and $k_{\mathbb{ij}}$, one for each possible
direction. Since the existence of the forward transition rate
$k_{\mathbb{ji}}>0$ implies the reverse transition $k_{\mathbb{ij}}>0$ in all
practical circumstances, we may assume without loss of any generality that our
basic graph $\mathbb{G}$ is connected in the sense that for each pair of
states or vertices $(\mathbb{i},\mathbb{j}),\mathbb{i}\neq\mathbb{j}$, there
exists at least one sequence of transitions or edges that connects them in
both directions. If $\mathbb{G}$ is not connected, the physical system
represented by $\mathbb{G}$ can be decomposed into non-interacting subsystems,
which can be analyzed independently Schnakenberg (1976).

Figure 2: The fundamental basic graph ($\mathbb{G}$) as well as its undirected
subcycles, encompasses all possible quantum transport channels (networks) of
the photovoltaic solar cell under non-equilibrium conditions.

Despite the fact that the underlying master equation (26) is linear and the
uniqueness of the steady-state solution is guaranteed by the properties of the
master equation, finding out the complete analytical solutions of the steady-
state populations is a nontrivial task even for the simplest physical
conditions Schnakenberg (1976). However, the diagrammatic representation of
the system in terms of its basic graph offers a highly versatile and effective
tool to tackle such complicated problems Ren (2017); Dutta  _et al._ (2020);
Wang  _et al._ (2022). For instance, the steady state solution of
$P_{\mathbb{i}}$ of our basic graph, is defined as Hill and Chen (1975);
Kohler and Vollmerhaus (1980); Schnakenberg (1976)

$\bar{P}_{\mathbb{i}}=\Lambda_{\mathbb{i}}/\Lambda,$ (27)

where $\Lambda_{\mathbb{i}}$ is the sum of the weight of the spanning trees
rooted on $\mathbb{i}^{th}$ state and $\Lambda$ is the sum of weights of
spanning trees rooted on every individual state, i.e.,
$\Lambda=\sum^{\mathbb{5}}_{\mathbb{i=0}}\Lambda_{\mathbb{i}}$. The above
method was first invented by Kirchhoff in 1847 Kirchhoff (1847) in the context
of network theory and later on, rigorously formulated by King and Altmann King
and Altman (1956) in the context of biochemical reactions. As a result, it is
commonly referred to as Kirchhoff’s theorem in the theory of network analysis
Kirchhoff (1847); King and Altman (1956); Schnakenberg (1976); Dutta  _et al._
(2020). Here by spanning tree, we refer to a covering subgraph of the basic
graph that contains all the vertices with the minimum number of edges which is
always connected but contains no circuits (cyclic sequence of edges or cycle
trajectories). From Fig. 2, it is easy to understand that a basic graph
generally contains a large number of undirected subcycles, and each subcycle
is a combination of a pair of two one-directional circuits or cycle
trajectories, $C^{+}$ (counterclockwise) and $C^{-}$ (clockwise). The net
cycle flux is therefore given by the difference between the two circuit fluxes
such as $J_{C}=J^{+}_{C}-J^{-}_{C}$. The notion of the “circuit flux” was
introduced by Kohler and Vollmerhaus Kohler and Vollmerhaus (1980), and has
since been widely employed to characterize a large variety of biological
systems Qian (2007); Ren (2017); Barato and Seifert (2015); Dutta  _et al._
(2020); King and Altman (1956); Giebink  _et al._ (2011); Rutten  _et al._
(2009). In essence, it captures the frequency of circuit completion along a
specific cycle trajectory ($C^{+}$ or $C^{-}$) and can be used to quantify the
edge flux
$J_{\mathbb{i}\rightarrow\mathbb{j}}=\sum_{\mathcal{C}}J^{+}_{\mathcal{C}}-J^{-}_{\mathcal{C}}$,
as the summation of the differences between the circuit fluxes along all cycle
trajectories that include the edge $\mathbb{i}\rightarrow\mathbb{j}$. For
example, the net edge flux $J_{\mathbb{2}\rightarrow\mathbb{3}}$ in Fig. 2 can
be obtained by
$J_{\mathbb{2}\rightarrow\mathbb{3}}=(J^{+}_{C_{1}}+J^{+}_{C_{3}}+J^{+}_{C_{4}})-(J^{-}_{C_{1}}+J^{-}_{C_{3}}+J^{-}_{C_{4}})$.

Finally, the flux associated with each one-directional “circuit” is determined
by Hill and Chen (1975); Kohler and Vollmerhaus (1980); Schnakenberg (1976)

$J^{\pm}_{C}=\Pi^{\pm}_{C}\frac{\Lambda_{C}}{\Lambda}.$ (28)

Here, $\Pi^{\pm}_{C}$ represents the weight factor which is defined by the
product of transition rates along the cycle trajectory $C^{\pm}$, whereas the
sum of the weights of the spanning trees rooted on cycle $C$ is given by
$\Lambda_{C}$ and $\Lambda$ measures the total weight of the spanning trees
rooted on each individual state. As an example, Fig 3, displays three such
cycles of our basic graph $\mathbb{G}$ and the corresponding spanning trees
that are rooted on them, where the cycle $C^{+}_{1}$ possesses a weight factor
$\Pi^{+}_{C_{1}}=k_{\mathbb{03}}k_{\mathbb{32}}k_{\mathbb{21}}k_{\mathbb{10}}$.
In reality, enumerating the vast number of spanning trees that are rooted in
each individual vertices requires an inconceivable amount of work,
particularly as the graph size escalates. Fortunately, we can navigate this
problem by leveraging the generalized matrix-tree theorem from the algebraic
graph Wang  _et al._ (2022); Ren (2017). Upon rewriting the master equation in
the following form: $\frac{d\textbf{P}}{dt}=-\textbf{M}\textbf{P}$, where
$\textbf{P}=\\{P_{\mathbb{0}},P_{\mathbb{1}},P_{\mathbb{2}},P_{\mathbb{3}},P_{\mathbb{4}},P_{\mathbb{5}}\\}$
is a column matrix and

$\displaystyle\textbf{M}=\left[{\begin{array}[]{cccccc}k_{\mathbb{10}}+k_{\mathbb{30}}&-k_{\mathbb{01}}&0&-k_{\mathbb{03}}&0&0\\\
-k_{\mathbb{10}}&k_{\mathbb{01}}+k_{\mathbb{41}}+k_{\mathbb{21}}&-k_{\mathbb{12}}&0&-k_{\mathbb{14}}&0\\\
0&-k_{\mathbb{21}}&k_{\mathbb{12}}+k_{\mathbb{32}}+k_{\mathbb{52}}&-k_{\mathbb{23}}&0&-k_{\mathbb{25}}\\\
-k_{\mathbb{30}}&0&-k_{\mathbb{32}}&k_{\mathbb{03}}+k_{\mathbb{23}}+k_{\mathbb{43}}&-k_{\mathbb{34}}&0\\\
0&-k_{\mathbb{41}}&0&-k_{\mathbb{43}}&k_{\mathbb{34}}+k_{\mathbb{14}}+k_{\mathbb{54}}&-k_{\mathbb{45}}\\\
0&0&-k_{\mathbb{52}}&0&-k_{\mathbb{54}}&k_{\mathbb{25}}+k_{\mathbb{45}}\\\
\end{array}}\right],$ (35)

is the Laplacian (transition) matrix of the weighted graph, the matrix-tree
theorem provides a powerful recipe for calculating the number of directed
spanning trees rooted on a particular cycle. Specifically, the matrix-tree
theorem asserts that an effective expression for $\Lambda_{C}$ can be
evaluated as the determinant of the principal minor of the Laplacian matrix M
of the basic graph, namely $\det(\textbf{M}[C;C])$. In other words, the
determinant of the reduced matrix $\textbf{M}[C;C]$ obtained by removing the
rows and columns indexed by $\mathbb{i}\in C$ of the original matrix M, is
equal to the sum of the weights of directed spanning trees rooted on C, i.e.,
$\Lambda_{C}=\det(\textbf{M}[C;C])$. As an example, for the cycle $C_{1}$, the
reduced Laplacian matrix $\textbf{M}[C_{1};C_{1}]$ or
$\textbf{M}[\mathbb{0,1,2,3;0,1,2,3}]$ and its determinant are given by:

$\displaystyle\textbf{M}[C_{1};C_{1}]=\left[{\begin{array}[]{cc}k_{\mathbb{34}}+k_{\mathbb{14}}+k_{\mathbb{54}}&-k_{\mathbb{45}}\\\
-k_{\mathbb{54}}&k_{\mathbb{25}}+k_{\mathbb{45}}\\\ \end{array}}\right],$ (39)
$\displaystyle\det(\textbf{M}[C_{1};C_{1}])=k_{\mathbb{25}}(k_{\mathbb{34}}+k_{\mathbb{14}}+k_{\mathbb{54}})+k_{\mathbb{45}}(k_{\mathbb{34}}+k_{\mathbb{14}}).$

Likewise, in the case of $\Lambda_{\mathbb{i}}$, the determinant related to
the principal minor of the M matrix can be derived by excluding the relevant
row and column that correspond to state $\mathbb{i}$. Consequently,
$\Lambda_{\mathbb{i}}$ can be represented by the determinant of
$\textbf{M}[\mathbb{i};\mathbb{i}]$ and the steady-state population,
$\bar{P}_{\mathbb{i}}$, can be precisely written as the ratio of the
determinants Wang  _et al._ (2022); Ren (2017)

$\bar{P}_{\mathbb{i}}=\frac{\Lambda_{\mathbb{i}}}{\Lambda}\equiv\frac{\det(\textbf{M}[\mathbb{i};\mathbb{i}])}{\sum_{\mathbb{i}}\det(\textbf{M}[\mathbb{i};\mathbb{i}])},$
(41)

where we identify $\Lambda$ as the sum of determinants of the principal minors
of M, i.e., $\sum_{\mathbb{i}}\textbf{M}[\mathbb{i};\mathbb{i}]$. Similarly,
the one-directional circuit flux can directly be calculated as a product of
two factors Wang  _et al._ (2022); Ren (2017)

$J^{\pm}_{C}=\Pi^{\pm}_{C}\frac{\det(\textbf{M}[C;C])}{\sum_{\mathbb{i}}\det(\textbf{M}[\mathbb{i};\mathbb{i}])}.$
(42)

The first coefficient captures the weight of the cycle in the particular
direction, for instance,
$C_{1}(\mathbb{0}\rightarrow\mathbb{1}\rightarrow\mathbb{2}\rightarrow\mathbb{3}\rightarrow\mathbb{0})$
in Fig 3(a), has a counterclockwise weight factor
$\Pi^{+}_{C_{1}}=k_{\mathbb{03}}k_{\mathbb{32}}k_{\mathbb{21}}k_{\mathbb{10}}$.
The subsequent term involves a ratio of two determinants: the numerator
$\det(\textbf{M}[C_{1}\\{\mathbb{0},\mathbb{1},\mathbb{2},\mathbb{3}\\};C_{1}\\{\mathbb{0},\mathbb{1},\mathbb{2},\mathbb{3}\\}])$
tallies with the weighted summation of all five spanning trees rooted on
$C_{1}$ (Fig 3b), while the denominator
$\sum_{\mathbb{i=0}}^{\mathbb{5}}\det(\textbf{M}[\mathbb{i};\mathbb{i}])$
serves as a constant of normalization factor via a common term, representing
the total weight of the spanning trees rooted on each individual state. In
summary, the graph-theoretic representation of the cycle flux can be
intuitively understood as the flow of weighted edges on spanning trees
directed towards a cycle, which is intricately linked to the frequencies at
which the cycle trajectory occurs. As a result, the determinant
$\det(\textbf{M}[C_{1};C_{1}])$ exemplified above, can be alternatively
computed by means of the determinant $\det(\textbf{N}[\mathbb{k};\mathbb{k}])$
of a “new graph” N which is obtained by merging the set of vertices
$\\{\mathbb{0},\mathbb{1},\mathbb{2},\mathbb{3}\\}$ into a new vertex
“$\mathbb{k}$” within the initial graph M Wang  _et al._ (2022), represented
in Fig 3(b) by the shaded region.

Figure 3:  (a) First three top-ranked cycle trajectories of our basic graph
$\mathbb{G}$ shown in Fig. 2. (b) Spanning trees rooted on the top-ranked
cycles $C_{1}$, $C_{3}$, and $C_{4}$.

Following the same procedure, we can evaluate the cycle fluxes of all the
cycle trajectories of our basic graph and efficiently ranks out the top-ranked
cycle fluxes. The generalized matrix-tree theorem thus provides valuable
machinery by accomplishing the behavior of weighted graphs and their
corresponding cycle fluxes. Notably, the proficiency of the cycle flux ranking
scheme improves with the size of the graph Wang  _et al._ (2022).

##  IV Results and Discussion

In what follows, we utilize the aforementioned ranking scheme to examine the
non-equilibrium transport characteristic of PV solar cells. To unfold the
working principle of the donor-acceptor molecular junctions, we begin by
decomposing the PV cell network into subcycles or a complete set of paired
cycle trajectories. In the present case, the basic graph [Fig. 2] consists of
5 subcycles or 10 paired one-directional circuits (cycle trajectories).
Subsequently, the efficient cycle flux ranking scheme [Cf. Eq. (42)] is
implemented to rank out the top-ranked cycle fluxes and identify the major
working cycles of the molecular solar cell. In Fig. 4, we have plotted all the
circuit fluxes against the bias voltage $U$ and notice that the cycle flux
trajectories $C^{+}_{1}$
($|\mathbb{1}\rangle\rightarrow|\mathbb{2}\rangle\rightarrow|\mathbb{3}\rangle\rightarrow|\mathbb{0}\rangle\rightarrow|\mathbb{1}\rangle$)
or
($|1010\rangle\rightarrow|0110\rangle\rightarrow|0011\rangle\rightarrow|0010\rangle\rightarrow|1010\rangle$)
and $C^{+}_{4}$
($|\mathbb{1}\rangle\rightarrow|\mathbb{2}\rangle\rightarrow|\mathbb{3}\rangle\rightarrow|\mathbb{4}\rangle\rightarrow|\mathbb{1}\rangle$)
or
($|1010\rangle\rightarrow|0110\rangle\rightarrow|0011\rangle\rightarrow|1011\rangle\rightarrow|1010\rangle$),
are comparable in magnitudes and can be classified as the first-ranked
circuits. The second-ranked circuit is found to be $C^{+}_{3}$
($|\mathbb{3}\rangle\rightarrow|\mathbb{4}\rangle\rightarrow|\mathbb{5}\rangle\rightarrow|\mathbb{2}\rangle\rightarrow|\mathbb{3}\rangle$)
or
($|0011\rangle\rightarrow|1011\rangle\rightarrow|0111\rangle\rightarrow|0110\rangle\rightarrow|0011\rangle$).
Both the first and second-ranked circuits are plotted in Fig. 4a, while
$C^{-}_{1}$
($|\mathbb{1}\rangle\rightarrow|\mathbb{0}\rangle\rightarrow|\mathbb{3}\rangle\rightarrow|\mathbb{2}\rangle\rightarrow|\mathbb{1}\rangle$)
or
($|1010\rangle\rightarrow|0010\rangle\rightarrow|0011\rangle\rightarrow|0110\rangle\rightarrow|1010\rangle$),
$C^{-}_{3}$
($|\mathbb{3}\rangle\rightarrow|\mathbb{2}\rangle\rightarrow|\mathbb{5}\rangle\rightarrow|\mathbb{4}\rangle\rightarrow|\mathbb{3}\rangle$)
or
($|0011\rangle\rightarrow|0110\rangle\rightarrow|0111\rangle\rightarrow|1011\rangle\rightarrow|0011\rangle$),
and $C^{-}_{4}$
($|\mathbb{1}\rangle\rightarrow|\mathbb{4}\rangle\rightarrow|\mathbb{3}\rangle\rightarrow|\mathbb{2}\rangle\rightarrow|\mathbb{1}\rangle$)
or
($|1010\rangle\rightarrow|1011\rangle\rightarrow|0011\rangle\rightarrow|0110\rangle\rightarrow|1010\rangle$),
are the lowest ranked cycle trajectories shown in Fig. 4b. On the contrary,
the third-ranked paired circuits $C^{\pm}_{2}$ correspond to sequence of
states
($|\mathbb{1}\rangle\leftrightarrow|\mathbb{4}\rangle\leftrightarrow|\mathbb{3}\rangle\leftrightarrow|\mathbb{0}\rangle\leftrightarrow|\mathbb{1}\rangle$)
or
($|1010\rangle\leftrightarrow|1011\rangle\leftrightarrow|0011\rangle\leftrightarrow|0010\rangle\leftrightarrow|1010\rangle$),
and the fourth-ranked paired circuits, $C^{\pm}_{5}$ correspond to sequence of
states
($|\mathbb{1}\rangle\leftrightarrow|\mathbb{2}\rangle\leftrightarrow|\mathbb{5}\rangle\leftrightarrow|\mathbb{4}\rangle\leftrightarrow|\mathbb{1}\rangle$)
or
($|1010\rangle\leftrightarrow|0110\rangle\leftrightarrow|0111\rangle\leftrightarrow|1011\rangle\leftrightarrow|1010\rangle$),
are plotted in Fig. 4c. From Fig. 4d, we can understand that circuit fluxes
corresponding to $C^{\pm}_{2}$ and $C^{\pm}_{5}$ are not only several orders
of magnitude smaller than the first ($C^{+}_{1}$ and $C^{+}_{4}$) and second-
ranked ($C^{+}_{3}$) cycle trajectories, they also have an equal amount of
flux current in both directions. Therefore, the cycle affinity Einax and
Nitzan (2016), $\mathcal{A}$ becomes zero upon the addition of two counter-
pair cycle trajectories, where the cycle affinity is defined as
$\mathcal{A}=-\ln{\mathcal{K}}$, $\mathcal{K}={\Pi^{+}_{C}}/{\Pi^{-}_{C}}$
being the ratio of the forward and the backward rates for a specific cycle.
Since the magnitude of the cycle flux is the difference between the circuit
fluxes $J_{C}=J^{+}_{C}-J^{-}_{C}$, it follows immediately that cycles
$C_{1}$, $C_{3}$ and $C_{4}$ possess nonzero flux current fulfilling
$\mathcal{K}\neq 1$, while cycle flux associated with $C_{2}$ and $C_{5}$ are
identically zero, satisfying $\mathcal{K}=1$.

Figure 4: The circuit fluxes $J^{\pm}_{C}$ (in Ampere) as a function of bias
voltage $U$ (in Volt) are plotted as follows: (a) top-ranked circuit fluxes
$J^{+}_{C_{1}}$, $J^{+}_{C_{3}}$, $J^{+}_{C_{4}}$ vs. $U$ (b) lowest ranked
circuit fluxes $J^{-}_{C_{1}}$, $J^{-}_{C_{3}}$, $J^{-}_{C_{4}}$ vs. $U$ (c)
third and fourth-ranked paired circuits $J^{\pm}_{C_{2}}$ and
$J^{\pm}_{C_{5}}$ vs. $U$, and (d) all circuit fluxes (in log scale) vs. $U$.
The parameters are $\mu_{\rm L}=0$ ${\rm eV}$, $\mu_{\rm R}=\mu_{\rm L}+|e|U$,
$\varepsilon_{\rm D1}=-0.1$ ${\rm eV}$, $\varepsilon_{\rm D2}=1.4$ ${\rm eV}$,
$\varepsilon_{\rm A2}=0.9$ ${\rm eV}$, $U_{\rm A}=0.25$ ${\rm eV}$, $T_{\rm
ph}=T_{\rm L}=T_{\rm R}=T$, $T=300$ K, $T_{\rm S}=6000$ K, $\gamma_{\rm
L}=\gamma_{\rm R}=\gamma^{\rm D}_{\rm pht}=\gamma^{\rm D}_{\rm
phn}=0.01\gamma^{\rm DA}_{\rm phn}$ and $\gamma^{\rm DA}_{\rm phn}=10^{12}$
${\rm s}^{-1}$ Einax  _et al._ (2011).

As a result, we conclude that the $C_{2}$ and $C_{5}$ cycles do not contribute
to the overall photocurrent, and the entire contribution to the electron
current generated within the photocell solely comes from the first ($C_{1}$
and $C_{4}$) and second-ranked ($C_{3}$) cycles. We have made a dual-axis plot
to illustrate the total current ($J(U)$) and the corresponding power
($P=UJ(U)$), where $P$ passes through a maximum in Fig. 5, as expected. Our
analysis reveals that the total current is really the sum of the top three
cycle fluxes, i.e, $J=J_{C_{1}}+J_{C_{4}}+J_{C_{3}}$, regardless of the
presence or the absence of the nonradiative loss processes due to $\gamma^{\rm
D}_{\rm phn}$. To make a proper comparison, we have taken the parameter set
used by Nitzan’s group Einax  _et al._ (2011) and found excellent agreement
with their results. Qualitatively, one can understand that among the five
possible cycles, those cycles ($C_{1}$, $C_{3}$ and $C_{4}$) which involve the
$|\mathbb{2}\rangle\leftrightarrow|\mathbb{3}\rangle$ transition, contribute a
finite amount to the overall electron current, whereas sub-cycles $C_{2}$ and
$C_{5}$ which do not possess this edge, have zero contribution to the total
current. This can be attributed to the fact that among various possible
electron transfer channels, those pathways are only relevant that do involve
the electron transfer between donor and acceptor molecules. In graph-
theoretical language, this corresponds to subcycles with nonzero edge flux
along $|\mathbb{2}\rangle\leftrightarrow|\mathbb{3}\rangle$. Furthermore, it
is evident from Fig. 5 that compared to first-ranked cycles, the second-ranked
cycle $C_{3}$ makes a negligible contribution to the overall electron current
for the majority of the parameter range. Consequently, the total flux current
$J$ is simply twice that of the two individual first-ranked cycles.
Intuitively, one expects the $C_{1}$ electron transfer channel to be the one
and only exclusive route for the electron transfer pathways in photovoltaic
devices Einax and Nitzan (2014, 2016). Surprisingly, our analysis reveals the
existence of an equally important $C_{4}$ pathway that contributes to a
similar extent to the overall photocurrent, along with the natural $C_{1}$
electron transfer channel. These findings represent the first set of important
results for our efficient ranking scheme, which goes beyond the conventional
understanding of the electron transfer pathways in molecular photovoltaics.

Figure 5: Dual axis plot of top-ranked cycle fluxes $J_{C_{1}}$ (azure blue),
$J_{C_{3}}$ (black), $J_{C_{4}}$ (green), and electron current $J$ (blue) on
the left vertical axis (in Ampere) as a function of voltage bias $U$ and $P$
(red) on the right vertical axis (in Watt). The parameters are the same as
mentioned in Fig. 4. The solid lines represents $J(P)$ in the presence of
nonradiative loss process ($\gamma^{\rm D}_{\rm phn}=10^{10}$ ${\rm s}^{-1}$)
whereas the dash-dot lines represents $J(P)$ in the absence of nonradiative
losses ($\gamma^{\rm D}_{\rm phn}=0$).

Secondly, it can be inferred from Fig.5 that $C_{1}$ and $C_{4}$ operate as
first-ranked cycles, while $C_{3}$ acts as a second-ranked cycle for bias
voltage $U$ around $\sim 1$ eV. After that, both the net current as well as
the power experience a simultaneous drop. The explanation for this behavior
can be acquired from the population plot shown in Fig.6, as a function of $U$.
For bias voltages $U\sim 1$ eV, $|1010\rangle$ or $|\mathbb{1}\rangle$ is the
maximally populated state, and both the highest-ranked cycles start with the
initial state $|1010\rangle$, represented by the curve $P_{\mathbb{1}}$ in
Fig. 6. In contrast, the second-ranked cycle $C_{3}$ starts with a less
populated initial state $|0011\rangle$, denoted by $P_{\mathbb{0}}$. Notably,
the first two steps of the highest-ranked cycles are identical: it starts with
the electron transfer between the levels ${\rm D1}$ to ${\rm D2}$ at the donor
site that involve radiative (photon-induced) and nonradiative (phonon-induced)
processes, while the second step designates the electron transfer between
donor and acceptor molecules. It is worth pointing out that the above two
steps of the highest-ranked cycles represent their counter-clockwise (forward)
cycle trajectories, namely, $C^{+}_{1}$ and $C^{+}_{4}$, respectively. This is
reasonable since the backward (clockwise) cycle flux trajectories are
negligible compared to their forward counterpart and because of this reason it
is the counter-clockwise or forward cycle trajectories that predominately
determines the underlying directions of the cycle flux of the individual
subcycles.

Figure 6: Populations of the different states as a function of $U$ (in Volts).
The parameters are $\mu_{\rm L}=0$ ${\rm eV}$, $\mu_{\rm R}=\mu_{\rm L}+|e|U$,
$\varepsilon_{\rm D1}=-0.1$ ${\rm eV}$, $\varepsilon_{\rm D2}=1.4$ ${\rm eV}$,
$\varepsilon_{\rm A2}=0.9$ ${\rm eV}$, $U_{\rm A}=0.25$ ${\rm eV}$, $T_{\rm
ph}=T_{\rm L}=T_{\rm R}=T$, $T=300$ K, $T_{\rm S}=6000$ K, $\gamma_{\rm
L}=\gamma_{\rm R}=\gamma^{\rm D}_{\rm pht}=\gamma^{\rm D}_{\rm
phn}=0.01\gamma^{\rm DA}_{\rm phn}$ and $\gamma^{\rm DA}_{\rm phn}=10^{12}$
${\rm s}^{-1}$ Einax  _et al._ (2011).

In particular, the dynamical steps of the cycle $C^{+}_{1}$ are as follows:
starting from the neutral state $|1010\rangle$ ($|\mathbb{1}\rangle$), the
system transits sequentially into $|0110\rangle$ ($|\mathbb{2}\rangle$) by
absorbing the photon (electron transfer between ${\rm D1}$ and ${\rm D2}$ at
the donor site), followed by $|\mathbb{3}\rangle$ or $|0011\rangle$ (via
phonon relaxation), $|\mathbb{0}\rangle$ or $|0010\rangle$ (one electron
tunnels from ${\rm A2}$ of the acceptor into the right electron reservoir) and
finally returns to its initial state $|\mathbb{1}\rangle$ or $|1010\rangle$
(where one electron tunnels from the left reservoir into ${\rm D1}$ of the
donor). Similarly, for cycle, $C^{+}_{4}$, the first two processes are common
with the cycle $C^{+}_{1}$ (i.e., $|1010\rangle\rightarrow|0110\rangle$
($|\mathbb{1}\rangle\rightarrow|\mathbb{2}\rangle$) and
$|0110\rangle\rightarrow|0011\rangle$
($|\mathbb{2}\rangle\rightarrow|\mathbb{3}\rangle$)). The system then changes
its state from $|0011\rangle$ to $|1011\rangle$ or
$|\mathbb{3}\rangle\rightarrow|\mathbb{4}\rangle$, i.e., one electron tunnel
from the left reservoir into ${\rm D1}$ of the donor, and finally, returns
backs to its initial state $|\mathbb{1}\rangle$ or $|1010\rangle$ (where one
electron tunnels from the ${\rm A2}$ of the acceptor into the right
reservoir). So, these are the two pertinent electron transfer pathways that
play equivalent roles in a typical donor-acceptor molecular junction
photocell. For all the numerical plots presented in Figs. 4\- 6, we have used
the same set the parameters as Nitzan et. al. Einax  _et al._ (2011), which
are also reasonable from the experimental perspective as reported in Ref. Rice
and Gartstein (1996); Soci  _et al._ (2007); Pensack and Asbury (2010).

To understand the significant drop in the total electron current and power,
after a certain bias voltage in Fig. 5, we need to examine the steady-state
populations as shown in Fig. 6. We observe that as long as the chemical
potential of the right reservoir, $\mu_{\rm R}$, reaches the energy level of
$\varepsilon_{\rm A2}$, the trend of the populations obeys
$P_{\mathbb{1}}>>P_{\mathbb{0}}>P_{\mathbb{3}}\simeq
P_{\mathbb{4}}>P_{\mathbb{2}}\simeq P_{\mathbb{5}}$. However, once $\mu_{\rm
R}$ exceeds the energy of the acceptor level ${\rm A2}$, the probability of
finding the system in states $|1010\rangle$ ($P_{\mathbb{1}}$), $|0010\rangle$
($P_{\mathbb{0}}$), and $|0110\rangle$ ($P_{\mathbb{2}}$) decreases sharply.
Meanwhile, the probability of finding the system in states $|1011\rangle$
($P_{\mathbb{4}}$) and $|0111\rangle$ ($P_{\mathbb{3}}$) increases rapidly.
This indicates that electron tunneling from the excited state of the acceptor
into the right reservoir becomes less favorable. As a result, the circuit
fluxes $J^{+}_{C_{1}}$ and $J^{+}_{C_{4}}$ remain constant until $\mu_{\rm R}$
reaches the energy level of $\varepsilon_{\rm A2}$, and then decrease steadily
after surpassing this energy level. We can also analyze the second-ranked
circuit flux, $J^{+}_{C_{3}}$, which remains constant until $\mu_{\rm R}$
reaches the energy level $\varepsilon_{\rm A2}$, and becomes maximum precisely
at the midpoint of $\varepsilon_{\rm D2}$ and $\varepsilon_{\rm A2}$, i.e.,
$1.15$ ${\rm eV}$. By carefully analyzing Fig. 5, we can see that we have
varied the bias voltages by changing $\mu_{\rm R}$ while keeping all other
parameters fixed. Now, Fig.2 shows that $\mu_{\rm R}$ governs transitions
between pairs of states where the acceptor’s excited state is occupied in at
least one of them. For the counter-clock $C_{1}$ cycle, the second last edge
($|\mathbb{3}\rangle\leftrightarrow|\mathbb{0}\rangle$) is controlled by
$\mu_{\rm R}$, whereas, for $C_{4}$, it is the last edge
($|\mathbb{4}\rangle\leftrightarrow|\mathbb{1}\rangle$) that is controlled by
the bias voltage. Thus, we get a single drop in the current vs. voltage
diagram as $\mu_{\rm R}$ crosses the value of $\varepsilon_{\rm A2}$. In the
case of the $C_{3}$ cycle, the last two edges
(($|\mathbb{4}\rangle\leftrightarrow|\mathbb{5}\rangle$) and
($|\mathbb{5}\rangle\leftrightarrow|\mathbb{2}\rangle$)) are controlled by
$\mu_{\rm R}$. Since $\mu_{\rm R}$ controls both transitions, the cycle
current will be maximum when both transitions are favored optimally Ghosh  _et
al._ (2022), i.e., at the average value of both energy levels,
$\varepsilon_{\rm D2}$ and $\varepsilon_{\rm A2}$. For the present choice of
parameters, this turns out to be precisely $1.15$ eV, as obtained in Fig. 5.

Figure 7: Plot of $\eta^{*}$ ($U_{max}J_{C}(U_{max})/\dot{Q}_{\rm S}$) vs
$\Delta\varepsilon$ (eV) in absence of any nonradiative loss processes and the
parameters are $\mu_{\rm L}=0$ ${\rm eV}$, $\mu_{\rm R}=\mu_{\rm L}+|e|U$,
$\varepsilon_{\rm D1}=-0.1$ ${\rm eV}$, $\varepsilon_{\rm D2}=1.4$ ${\rm eV}$,
$\varepsilon_{\rm A2}=\varepsilon_{\rm D2}-\Delta\varepsilon$, $U_{\rm
A}=0.25$ ${\rm eV}$, $T_{\rm ph}=T_{\rm L}=T_{\rm R}=T$, $T=300$ K, $T_{\rm
S}=6000$ K, $\gamma_{\rm L}=\gamma_{\rm R}=\gamma^{\rm D}_{\rm
pht}=0.01\gamma^{\rm DA}_{\rm phn}$, and $\gamma^{\rm DA}_{\rm phn}=10^{12}$
${\rm s}^{-1}$ Einax  _et al._ (2011). Red and black dashed lines represent
the Carnot and Curzon-Ahlborn bounds respectively. Figure 8: Plot of
$\eta^{*}$ ($U_{max}J_{C}(U_{max})/\dot{Q}_{\rm S}$) vs $\Delta\varepsilon$
(eV) in presence of nonradtive losses arising from $\gamma^{\rm D}_{\rm{phn}}$
and the parameters are $\mu_{\rm L}=0$ ${\rm eV}$, $\mu_{\rm R}=\mu_{\rm
L}+|e|U$, $\varepsilon_{\rm D1}=-0.1$ ${\rm eV}$, $\varepsilon_{\rm D2}=1.4$
${\rm eV}$, $\varepsilon_{\rm A2}=\varepsilon_{\rm D2}-\Delta\varepsilon$,
$U_{\rm A}=0.25$ ${\rm eV}$, $T_{\rm ph}=T_{\rm L}=T_{\rm R}=T$, $T=300$ K,
$T_{\rm S}=6000$ K, $\gamma_{\rm L}=\gamma_{\rm R}=\gamma^{\rm D}_{\rm
pht}=0.01\gamma^{\rm DA}_{\rm phn}$ and $\gamma^{\rm DA}_{\rm phn}=10^{12}$
${\rm s}^{-1}$ Einax  _et al._ (2011). Red and black dashed lines represent
the Carnot and Curzon-Ahlborn bounds respectively.

From the above analysis, it becomes clear that the donor-acceptor energy gap
plays a crucial role in determining the performance of the photocell. To this
end, in Figs. 7 and 8, we consider for our photocell, a real practical
quantity of interest, i.e., the efficiency at the maximum power as a function
of the donor-acceptor energy gap $\Delta\varepsilon=\varepsilon_{\rm
D2}-\varepsilon_{\rm A2}$. To be precise, we define the thermodynamic
efficiency Einax  _et al._ (2011); Rutten  _et al._ (2009); Giebink  _et al._
(2011) at maximum power as follows:

$\displaystyle\eta^{*}$ $\displaystyle=$
$\displaystyle\frac{P_{max}}{\dot{Q}_{\rm
S}}=\frac{U_{max}J(U_{max})}{\dot{Q}_{\rm S}},$ (43) $\displaystyle\approx$
$\displaystyle\frac{U_{max}J_{C_{1}}(U_{max})}{\dot{Q}_{\rm
S}}+\frac{U_{max}J_{C_{4}}(U_{max})}{\dot{Q}_{\rm S}},$ (44)

where $U_{max}$ denotes the bias voltage at which power passes through a
maximum, i.e., at $U=U_{max}$, power becomes $P_{max}=U_{max}J(U_{max})$. In
Eq. (43), the denominator $\dot{Q}_{\rm S}=\Delta EJ_{\rm S}(U_{max})$,
represents the total energy absorbed per unit time from the radiation field,
with $\Delta E$ being the energy gap between the donor excited and ground
levels. We rewrite the expression of efficiency at maximum power in terms of
the two first-ranked cycle fluxes in Eq. (44), where each term represents the
contribution of its respective cycle to the overall efficiency. From Figs. 7
and 8, we can conclude that the cycles $C_{1}$ and $C_{4}$ are the primary
electron transfer channels across the entire parameter space, regardless of
whether radiative process or non-radiative losses are at play. Interestingly,
the plots exhibit a distinctive peak in the behavior of $\eta^{*}$, in the
presence of non-radiative loss processes. In such situations, $\eta^{*}$
attains a maximum value as a function of $\Delta\varepsilon$, where the energy
gap $\Delta\varepsilon$ regulates the electron transfer process between donor
and acceptor [Fig. 8]. Moreover, it is worth noting that although $\eta^{*}$
is bounded by the Carnot ($\eta_{C}=1-T/T_{\rm S}$), it can exceed the
standard Curzon-Ahlborn bound Curzon and Ahlborn (1975)
($\eta_{CA}=1-\sqrt{T/T_{\rm S}}$) in the absence of any nonradiative losses
[Fig. 7]. In other words, the strength of nonradiative loss processes is found
to have a deleterious effect on the overall performance of the solar cell. It
may even reduce the efficiency at maximum power at a value well below the
$\eta_{CA}$ as shown in Fig. 8. Hence, in order to improve the device
performance to its ultimate level, it becomes crucial to diligently mitigate
all forms of nonradiative losses to a considerable extent.

##  V Conclusion

We introduce the top-ranked cycle flux ranking scheme of network analysis as a
tool to elucidate the complex working principles of molecular junction solar
cells. Our approach takes advantage of the mapping between the dissipative
Lindblad master equation for molecular systems and the quantum-transition
network that characterizes the nonequilibrium transport behavior of molecular
photocells. We now summarize our key insights as follows:

(i) We have provided a microscopic Hamiltonian description of the
phenomenological rate equations that are commonly used to characterize
molecular junction solar cells. Based on a minimal model Hamiltonian,
classical looking rate equations are derived from detailed quantum Lindblad
master equations. The resulting rate equations, though appear “classical”, the
underlying transition rates are shown to be quantum mechanical in nature.

(ii) With the implementation of the effective ranking scheme, we have
predicted the existence of a counterintuitive electron transfer pathway, which
provides valuable insights into the detailed working principles of molecular
photocells. Through rigorous analysis, we have clearly justified that cycles
that contribute to the overall photocurrent must involve nonzero edge flux
between donor and acceptor molecules embedded in a complex graph of quantum
transition networks. We provide explanations for the drop in the current and
power when measured against the bias voltage, as well as the various cycle
fluxes. Our result advances the conventional understanding of nonequilibrium
electron transfer pathways in donor-acceptor molecular junction solar cells
which might be a significant step toward making efficient photovoltaic devices
in the near future.

(iii) Finally, we obtain a crucial insight into the efficiency of photocells,
revealing that their maximum power efficiency in the absence of non-radiative
losses, can exceed the conventional Curzon-Ahlborn bound, yet abide by the
Carnot limit. Nonetheless, the incorporation of a non-radiative recombination
process at the donor site has been found to have a detrimental effect on the
photocell’s performance, reducing its efficiency at maximum power, below the
Curzon-Ahlborn limit. These findings underscore the intricacies of various
factors that govern the overall performance of molecular photovoltaics. While
in our present work, top-rankled cycle fluxes capture the essential physics of
molecular junction solar cells, future research can take into account other
additional effects, such as environment-assisted electron transfer and
recombination rates Sarkar and Dubi (2020), Marcus homogeneous and
heterogeneous electron transfer rates Einax  _et al._ (2013), from a similar
perspective. Thus the current approach offers important insights for further
research in the field of photovoltaics.

## Acknowledgement

AG acknowledges financial support from the Initiation grant of IITK (Grant No.
IITK/CHM/2018513). N.G. is thankful to CSIR for the fellowship, and it is
supported in part by the International Centre for Theoretical Sciences (ICTS)
during a visit for participating in the online program - Bangalore School on
Statistical Physics XII (code: ICTS/bssp2021/6). S.G. is grateful to the
Ministry of Education, Government of India, for the Prime Minister Research
Fellowship (PMRF).

## References

  * Nitzan (2006) A. Nitzan,  _Chemical Dynamics in Condensed Phases: Relaxation, Transfer and Reactions in Condensed Molecular Systems_ (Oxford University Press, 2006).
  * Bittner (2010) E. R. Bittner,  _Quantum Dynamics: Applications in Biological and Materials Systems_ (CRC Press, 2010).
  * Aradhya and Venkataraman (2013) S. V. Aradhya and L. Venkataraman, Nature Nanotechnology 8, 399 (2013).
  * Su  _et al._ (2016) T. A. Su, M. Neupane, M. L. Steigerwald, L. Venkataraman, and C. Nuckolls, Nature Reviews Materials 1, 16002 (2016).
  * Scully (2010) M. O. Scully, Phys. Rev. Lett. 104, 207701 (2010).
  * Scully  _et al._ (2011) M. O. Scully, K. R. Chapin, K. E. Dorfman, M. B. Kim, and A. Svidzinsky, Proc. Natl. Acad. Sci. USA 108, 15097 (2011).
  * Thoss and Evers (2018) M. Thoss and F. Evers, The Journal of Chemical Physics 148, 030901 (2018).
  * Chen (2005) G. Chen,  _Nanoscale Energy Transport and Conversion: A Parallel Treatment of Electrons, Molecules, Phonons, and Photons_ (Oxford University Press, Cambridge, England, 2005).
  * Pop (2010) E. Pop, Nano Research 3, 147 (2010).
  * Killoran  _et al._ (2015) N. Killoran, S. F. Huelga, and M. B. Plenio, The Journal of Chemical Physics 143 (2015), 10.1063/1.4932307, 155102.
  * Arp  _et al._ (2016) T. B. Arp, Y. Barlas, V. Aji, and N. M. Gabor, Nano Letters 16, 7461 (2016).
  * Deibel and Dyakonov (2010) C. Deibel and V. Dyakonov, Reports on Progress in Physics 73, 096401 (2010).
  * Nicholson and Castro (2010) P. G. Nicholson and F. A. Castro, Nanotechnology 21, 492001 (2010).
  * Kamat (2013) P. V. Kamat, The Journal of Physical Chemistry Letters 4, 908 (2013).
  * Altaner  _et al._ (2012) B. Altaner, S. Grosskinsky, S. Herminghaus, L. Katthän, M. Timme, and J. Vollmer, Phys. Rev. E 85, 041133 (2012).
  * Simine and Segal (2012) L. Simine and D. Segal, Phys. Chem. Chem. Phys. 14, 13820 (2012).
  * Agarwalla  _et al._ (2015) B. K. Agarwalla, J.-H. Jiang, and D. Segal, Phys. Rev. B 92, 245418 (2015).
  * Ajisaka  _et al._ (2015) S. Ajisaka, B. Žunkovič, and Y. Dubi, Scientific Reports 5, 8312 (2015).
  * Creatore  _et al._ (2013) C. Creatore, M. A. Parker, S. Emmott, and A. W. Chin, Phys. Rev. Lett. 111, 253601 (2013).
  * Fruchtman  _et al._ (2016) A. Fruchtman, R. Gómez-Bombarelli, B. W. Lovett, and E. M. Gauger, Phys. Rev. Lett. 117, 203603 (2016).
  * Giazotto and Martínez-Pérez (2012) F. Giazotto and M. J. Martínez-Pérez, Nature 492, 401 (2012).
  * Joulain  _et al._ (2016) K. Joulain, J. Drevillon, Y. Ezzahri, and J. Ordonez-Miranda, Phys. Rev. Lett. 116, 200601 (2016).
  * Whitney  _et al._ (2018) R. S. Whitney, R. Sánchez, and J. Splettstoesser, “Quantum thermodynamics of nanoscale thermoelectrics and electronic devices,” in  _Thermodynamics in the Quantum Regime: Fundamental Aspects and New Directions_, edited by F. Binder, L. A. Correa, C. Gogolin, J. Anders, and G. Adesso (Springer International Publishing, Cham, 2018) pp. 175–206.
  * Fong  _et al._ (2019) K. Y. Fong, H.-K. Li, R. Zhao, S. Yang, Y. Wang, and X. Zhang, Nature 576, 243 (2019).
  * Pekola and Karimi (2021) J. P. Pekola and B. Karimi, Rev. Mod. Phys. 93, 041001 (2021).
  * Breuer and Petruccione (2007) H.-P. Breuer and F. Petruccione,  _The Theory of Open Quantum Systems_ (Oxford University Press, 2007).
  * Agarwal (2012) G. S. Agarwal,  _Quantum Optics_ (Cambridge University Press, 2012).
  * Sinha  _et al._ (2013) S. S. Sinha, A. Ghosh, and D. S. Ray, Phys. Rev. E 87, 042112 (2013).
  * Gupt  _et al._ (2021) N. Gupt, S. Bhattacharyya, and A. Ghosh, Phys. Rev. E 104, 054130 (2021).
  * Gelbwaser-Klimovsky  _et al._ (2015) D. Gelbwaser-Klimovsky, W. Niedenzu, and G. Kurizki, Adv. At. Mol. Opt. Phys. 64, 329 (2015).
  * Kosloff (2013) R. Kosloff, Entropy 15, 2100 (2013).
  * Deffner and Campbell (2019) S. Deffner and S. Campbell,  _Quantum Thermodynamics_, 2053-2571 (Morgan & Claypool Publishers, 2019).
  * Dong  _et al._ (2021) H. Dong, A. Ghosh, M. B. Kim, S.-W. Li, A. A. Svidzinsky, Z. Zhang, G. Kurizki, and M. O. Scully, The European Physical Journal Special Topics 230, 867 (2021).
  * Gupt  _et al._ (2022) N. Gupt, S. Bhattacharyya, B. Das, S. Datta, V. Mukherjee, and A. Ghosh, Phys. Rev. E 106, 024110 (2022).
  * Ghosh  _et al._ (2022) S. Ghosh, N. Gupt, and A. Ghosh, Entropy 24 (2022), 10.3390/e24121810.
  * Kumar  _et al._ (2023) S. Kumar, N. Gupt, S. Ghosh, and A. Ghosh, Journal of Chemical Sciences 135, 36 (2023).
  * Aviram and Ratner (1974) A. Aviram and M. A. Ratner, Chemical Physics Letters 29, 277 (1974).
  * Harbola  _et al._ (2006) U. Harbola, M. Esposito, and S. Mukamel, Phys. Rev. B 74, 235309 (2006).
  * Galperin  _et al._ (2007) M. Galperin, M. A. Ratner, and A. Nitzan, Journal of Physics: Condensed Matter 19, 103201 (2007).
  * Wang  _et al._ (2020) Y.-S. Wang, P. Nijjar, X. Zhou, D. I. Bondar, and O. V. Prezhdo, The Journal of Physical Chemistry B 124, 4326 (2020).
  * Sarkar and Dubi (2020) S. Sarkar and Y. Dubi, The Journal of Physical Chemistry C 124, 15115 (2020).
  * Oz  _et al._ (2023) A. Oz, A. Nitzan, O. Hod, and J. E. Peralta, (2023), arXiv:2305.04343 .
  * Einax  _et al._ (2011) M. Einax, M. Dierl, and A. Nitzan, The Journal of Physical Chemistry C 115, 21396 (2011).
  * Einax  _et al._ (2013) M. Einax, M. Dierl, P. R. Schiff, and A. Nitzan, Europhysics Letters 104, 40002 (2013).
  * Einax and Nitzan (2014) M. Einax and A. Nitzan, The Journal of Physical Chemistry C 118, 27226 (2014).
  * Einax and Nitzan (2016) M. Einax and A. Nitzan, The Journal of Chemical Physics 145, 014108 (2016).
  * Craven  _et al._ (2018) G. T. Craven, D. He, and A. Nitzan, Phys. Rev. Lett. 121, 247704 (2018).
  * Mazza  _et al._ (2014) F. Mazza, R. Bosisio, G. Benenti, V. Giovannetti, R. Fazio, and F. Taddei, New Journal of Physics 16, 085001 (2014).
  * Wang  _et al._ (2022) L. Wang, Z. Wang, C. Wang, and J. Ren, Phys. Rev. Lett. 128, 067701 (2022).
  * Zurek (1981) W. H. Zurek, Phys. Rev. D 24, 1516 (1981).
  * Wu and Cao (2013) J. Wu and J. Cao, The Journal of Chemical Physics 139, 044102 (2013).
  * Chernyak and Sinitsyn (2008) V. Y. Chernyak and N. A. Sinitsyn, Phys. Rev. Lett. 101, 160601 (2008).
  * Rahav  _et al._ (2008) S. Rahav, J. Horowitz, and C. Jarzynski, Phys. Rev. Lett. 101, 140602 (2008).
  * Hill and Chen (1975) T. L. Hill and Y. D. Chen, Proceedings of the National Academy of Sciences 72, 1291 (1975).
  * Kohler and Vollmerhaus (1980) H.-H. Kohler and E. Vollmerhaus, Journal of Mathematical Biology 9, 275 (1980).
  * Schnakenberg (1976) J. Schnakenberg, Rev. Mod. Phys. 48, 571 (1976).
  * Tutte (2001) W. T. Tutte,  _Graph Theory_ , Vol. 21 (Cambridge University Press, 2001).
  * Balakrishnan and Ranganathan (2012) R. Balakrishnan and K. Ranganathan,  _A Textbook of Graph Theory_ (Springer Science & Business Media, 2012).
  * Polettini and Esposito (2014) M. Polettini and M. Esposito, Journal of Statistical Mechanics: Theory and Experiment 2014, P10033 (2014).
  * Owen  _et al._ (2020) J. A. Owen, T. R. Gingrich, and J. M. Horowitz, Phys. Rev. X 10, 011066 (2020).
  * Ghosh  _et al._ (2012) A. Ghosh, S. S. Sinha, and D. S. Ray, Phys. Rev. E 86, 011138 (2012).
  * Ren (2017) J. Ren, Frontiers of Physics 12, 120505 (2017).
  * Dutta  _et al._ (2020) A. Dutta, G. M. Schütz, and D. Chowdhury, Phys. Rev. E 101, 032402 (2020).
  * Kirchhoff (1847) G. Kirchhoff, Poggendorffs Ann. Phys. Chem. 148, 497 (1847).
  * King and Altman (1956) E. L. King and C. Altman, The Journal of Physical Chemistry 60, 1375 (1956).
  * Qian (2007) H. Qian, Annual Review of Physical Chemistry 58, 113 (2007).
  * Barato and Seifert (2015) A. C. Barato and U. Seifert, Phys. Rev. Lett. 114, 158101 (2015).
  * Giebink  _et al._ (2011) N. C. Giebink, G. P. Wiederrecht, M. R. Wasielewski, and S. R. Forrest, Phys. Rev. B 83, 195326 (2011).
  * Rutten  _et al._ (2009) B. Rutten, M. Esposito, and B. Cleuren, Phys. Rev. B 80, 235122 (2009).
  * Rice and Gartstein (1996) M. J. Rice and Y. N. Gartstein, Phys. Rev. B 53, 10764 (1996).
  * Soci  _et al._ (2007) C. Soci, I.-W. Hwang, D. Moses, Z. Zhu, D. Waller, R. Gaudiana, C. Brabec, and A. Heeger, Advanced Functional Materials 17, 632 (2007).
  * Pensack and Asbury (2010) R. D. Pensack and J. B. Asbury, The Journal of Physical Chemistry Letters 1, 2255 (2010).
  * Curzon and Ahlborn (1975) F. L. Curzon and B. Ahlborn, Am. J. Phys. 43, 22 (1975).

##  Appendix A Derivation of the Lindblad Master Equation and the electron
currents

The total Hamiltonian of the system interacting with the environments (baths)
can be written as

$H_{T}=H_{el}+H_{B}+H_{SB},$ (45)

where $H_{el}$ is the Hamiltonian of the molecular system and $H_{B}$ is the
total bath Hamiltonian which is given by $H_{B}=H_{\rm L}+H_{\rm R}+H_{\rm
pht}+H_{\rm phn}$. The interaction Hamiltonian between different baths and
molecular sites is given by

$H_{SB}=H_{I}+H^{\rm pht}_{I}+H^{\rm phn}_{I}.$ (46)

Here, $H_{I}=H^{\rm L}_{I}+H^{\rm R}_{I}$ represents the interaction of the
molecular sites (donor HOMO and acceptor LUMO) with the left and right
electrodes. $H^{\rm pht}_{I}$ characterizes the interaction of the photon bath
with the donor site only, and $H^{\rm phn}_{I}$ denotes the interaction of the
phonon bath with the respective molecular sites Ajisaka  _et al._ (2015);
Sarkar and Dubi (2020); Agarwalla  _et al._ (2015). For convenience, we
rewrite the interaction Hamiltonians [Eqs. (3)-(5)] of the main text as

$\displaystyle H^{\rm L}_{I}=\sum_{l}\hbar g_{l}(d^{\dagger}_{l}c_{\rm
D1}+c^{\dagger}_{\rm D1}d_{l});\quad H^{\rm R}_{I}=\sum_{r}\hbar
g_{r}(d^{\dagger}_{r}c_{\rm A2}+c^{\dagger}_{\rm A2}d_{r}),\quad H^{\rm
pht}_{I}=\sum_{k}\hbar g^{\rm D}_{k}(a^{\dagger}_{k}c^{\dagger}_{\rm D1}c_{\rm
D2}+c^{\dagger}_{\rm D2}c_{\rm D1}a_{k}),$ $\displaystyle H^{\rm
phn}_{I}=\sum_{q}\hbar g^{\rm D}_{q}(b^{\dagger}_{q}c^{\dagger}_{\rm D1}c_{\rm
D2}+c^{\dagger}_{\rm D2}c_{\rm D1}b_{q})+\sum_{q}\hbar g^{\rm
DA}_{q}(b^{\dagger}_{q}c^{\dagger}_{\rm A2}c_{\rm D2}+c^{\dagger}_{\rm
D2}c_{\rm A2}b_{q}).$ (47)

In order to derive the master equation, one can start with the von Neumann
equation for the total density matrix $\rho_{T}$ in the interaction picture
which is given by

$\frac{d\rho_{T}}{dt}=-\frac{i}{\hbar}[H_{SB}(t),\rho_{T}(t)].$ (48)

Under the Born-Markov approximation, the master equation in terms of the
reduced density matrix $\rho$ of the system can be written as Breuer and
Petruccione (2007); Gupt  _et al._ (2022); Ghosh  _et al._ (2022)

$\frac{d\rho(t)}{dt}=-\frac{1}{\hbar^{2}}\Tr_{B}\int^{\infty}_{0}ds[H_{SB}(t),[H_{SB}(t-s),\rho_{T}(t)]],$
(49)

where $\rho(t)=\Tr_{B}\\{\rho_{T}\\}\equiv\Tr_{{\rm L},{\rm R},{\rm pht},{\rm
phn}}\\{\rho(t)\otimes\rho_{\rm L}\otimes\rho_{\rm R}\otimes\rho_{\rm
pht}\otimes\rho_{\rm phn}\\}$ and $\Tr_{{\rm L},{\rm R},{\rm pht},{\rm phn}}$
stands for the trace over each bath degrees of freedom. As a result, we can
rewrite the Eq. (49) as

$\frac{d\rho(t)}{dt}=-\frac{1}{\hbar^{2}}\Tr_{{\rm L},{\rm R},{\rm pht},{\rm
phn}}\int^{\infty}_{0}ds[H_{SB}(t),[H_{SB}(t-s),\rho(t)\otimes\rho_{\rm
L}\otimes\rho_{\rm R}\otimes\rho_{\rm pht}\otimes\rho_{\rm phn}]].$ (50)

Since the bath operators obey the following relations

$\displaystyle\Tr_{\rm L}(d_{l}(t)\rho_{\rm L})$ $\displaystyle=$
$\displaystyle 0=\Tr_{\rm L}(d^{\dagger}_{l}(t)\rho_{\rm
L}),\quad\quad\quad\Tr_{\rm R}(d_{r}(t)\rho_{\rm R})=0=\Tr_{\rm
R}(d^{\dagger}_{r}(t)\rho_{\rm R}),$ $\displaystyle\Tr_{\rm
pht}(a_{k}(t)\rho_{\rm pht})$ $\displaystyle=$ $\displaystyle 0=\Tr_{\rm
pht}(a^{\dagger}_{k}(t)\rho_{\rm pht}),\quad\Tr_{\rm phn}(b_{q}(t)\rho_{\rm
phn})=0=\Tr_{\rm phn}(b^{\dagger}_{q}(t)\rho_{\rm phn}),$ (51)

one can verify that Gupt  _et al._ (2022); Ghosh  _et al._ (2022)

$\Tr_{{\rm L},{\rm R},{\rm pht},{\rm
phn}}\\{[H^{\alpha}_{I}(t),[H^{\beta}_{I}(t-s),\rho(t)\otimes\rho_{\rm
L}\otimes\rho_{\rm R}\otimes\rho_{\rm pht}\otimes\rho_{\rm
phn}]]\\}=0;\quad\alpha\neq\beta,\quad\alpha,\beta={\rm L},{\rm R},{\rm
pht},{\rm phn}.$ (52)

This simplifies Eq. (50) to

$\frac{d\rho(t)}{dt}=-\frac{1}{\hbar^{2}}\sum_{\alpha}\Tr_{{\rm L},{\rm
R},{\rm pht},{\rm
phn}}\int^{\infty}_{0}ds[H^{\alpha}_{I}(t),[H^{\alpha}_{I}(t-s),\rho(t)\otimes\rho_{\rm
L}\otimes\rho_{\rm R}\otimes\rho_{\rm pht}\otimes\rho_{\rm phn}]].$ (53)

In the above equation, we use the interaction picture system operators as
Breuer and Petruccione (2007)

$\displaystyle c_{\rm si}(t)$ $\displaystyle=$ $\displaystyle
e^{{iH_{el}t}/{\hbar}}c_{\rm
si}e^{{-iH_{el}t}/{\hbar}}=\sum_{\\{\varepsilon_{\mathbb{ji}}\\}}e^{{-i{\varepsilon_{\mathbb{ji}}}t}/{\hbar}}c_{\rm
si};$ $\displaystyle c^{\dagger}_{\rm si}(t)$ $\displaystyle=$ $\displaystyle
e^{{iH_{el}t}/{\hbar}}c^{\dagger}_{\rm
si}e^{{-iH_{el}t}/{\hbar}}=\sum_{\\{\varepsilon_{\mathbb{ji}}\\}}e^{{i{\varepsilon_{\mathbb{ji}}}t}/{\hbar}}c^{\dagger}_{\rm
si};\quad\quad{\rm s=D,A}\quad\text{and}\quad{\rm i}=1,2;$ (54)

where
$\varepsilon_{\mathbb{ji}}=\varepsilon_{\mathbb{j}}-\varepsilon_{\mathbb{i}}>0$
is the transition energy associated with the $\mathbb{j}$-th and
$\mathbb{i}$-th states of the molecular system. Similarly, one can also
evaluate the expressions for the interaction picture bath operators Breuer and
Petruccione (2007). With the above prescriptions, we can derive after a little
bit of algebra, the compact form of the following Lindblad master equation

$\frac{d\rho}{dt}=\mathcal{L}_{\rm L}[\rho]+\mathcal{L}_{\rm
R}[\rho]+\mathcal{L}^{\rm D}_{\rm pht}[\rho]+\mathcal{L}^{\rm D}_{\rm
phn}[\rho]+\mathcal{L}^{\rm DA}_{\rm phn}[\rho].$ (55)

The explicit forms of the Lindblad superoperator $\mathcal{L}$ in the above
equation are given by

$\displaystyle\mathcal{L}_{\rm L}[\rho]=\sum_{\\{\varepsilon_{\rm
L}\\}}\gamma_{\rm L}\Big{[}f(\varepsilon_{\rm L},\mu_{\rm L},T_{\rm
L})\Big{(}c^{\dagger}_{\rm D1}(\varepsilon_{\rm L})\rho c_{\rm
D1}(\varepsilon_{\rm L})-\frac{1}{2}\\{c_{\rm D1}(\varepsilon_{\rm
L})c^{\dagger}_{\rm D1}(\varepsilon_{\rm L}),\rho\\}\Big{)}$
$\displaystyle+(1-f(\varepsilon_{\rm L},\mu_{\rm L},T_{\rm L}))\Big{(}c_{\rm
D1}(\varepsilon_{\rm L})\rho c^{\dagger}_{\rm D1}(\varepsilon_{\rm
L})-\frac{1}{2}\\{c^{\dagger}_{\rm D1}(\varepsilon_{\rm L})c_{\rm
D1}(\varepsilon_{\rm L}),\rho\\}\Big{)}\Big{]},$ (56)
$\displaystyle\mathcal{L}_{\rm R}[\rho]=\sum_{\\{\varepsilon_{\rm
R}\\}}\gamma_{\rm R}\Big{[}f(\varepsilon_{\rm R},\mu_{\rm R},T_{\rm
R})\Big{(}c^{\dagger}_{\rm A2}(\varepsilon_{\rm R})\rho c_{\rm
A2}(\varepsilon_{\rm R})-\frac{1}{2}\\{c_{\rm A2}(\varepsilon_{\rm
R})c^{\dagger}_{\rm A2}(\varepsilon_{\rm R}),\rho\\}\Big{)}$
$\displaystyle+(1-f(\varepsilon_{\rm R},\mu_{\rm R},T_{\rm R}))\Big{(}c_{\rm
A2}(\varepsilon_{\rm R})\rho c^{\dagger}_{\rm A2}(\varepsilon_{\rm
R})-\frac{1}{2}\\{c^{\dagger}_{\rm A2}(\varepsilon_{\rm R})c_{\rm
A2}(\varepsilon_{\rm R}),\rho\\}\Big{)}\Big{]},$ (57)
$\displaystyle\mathcal{L}^{\rm D}_{\rm pht}[\rho]=\sum_{\\{\varepsilon_{\rm
D}\\}}\gamma^{\rm D}_{\rm pht}\Big{[}n(\varepsilon_{\rm D},T_{\rm
S})\Big{(}V^{\dagger}_{\rm D}(\varepsilon_{\rm D})\rho V_{\rm
D}(\varepsilon_{\rm D})-\frac{1}{2}\\{V_{\rm D}(\varepsilon_{\rm
D})V^{\dagger}_{\rm D}(\varepsilon_{\rm D}),\rho\\}\Big{)}$
$\displaystyle+(n(\varepsilon_{\rm D},T_{\rm S})+1)\Big{(}V_{\rm
D}(\varepsilon_{\rm D})\rho V^{\dagger}_{\rm D}(\varepsilon_{\rm
D})-\frac{1}{2}\\{V^{\dagger}_{\rm D}(\varepsilon_{\rm D})V_{\rm
D}(\varepsilon_{\rm D}),\rho\\}\Big{)}\Big{]},$ (58)
$\displaystyle\mathcal{L}^{\rm D}_{\rm phn}[\rho]=\sum_{\\{\varepsilon_{\rm
D}\\}}\gamma^{\rm D}_{\rm phn}\Big{[}n(\varepsilon_{\rm D},T_{\rm
ph})\Big{(}V^{\dagger}_{\rm D}(\varepsilon_{\rm D})\rho V_{\rm
D}(\varepsilon_{\rm D})-\frac{1}{2}\\{V_{\rm D}(\varepsilon_{\rm
D})V^{\dagger}_{\rm D}(\varepsilon_{\rm D}),\rho\\}\Big{)}$
$\displaystyle+(n(\varepsilon_{\rm D},T_{\rm ph})+1)\Big{(}V_{\rm
D}(\varepsilon_{\rm D})\rho V^{\dagger}_{\rm D}(\varepsilon_{\rm
D})-\frac{1}{2}\\{V^{\dagger}_{\rm D}(\varepsilon_{\rm D})V_{\rm
D}(\varepsilon_{\rm D}),\rho\\}\Big{)}\Big{]},$ (59)
$\displaystyle\mathcal{L}^{\rm DA}_{\rm phn}[\rho]=\sum_{\\{\varepsilon_{\rm
DA}\\}}\gamma^{\rm DA}_{\rm phn}\Big{[}n(\varepsilon_{\rm DA},T_{\rm
ph})\Big{(}V^{\dagger}_{\rm DA}(\varepsilon_{\rm DA})\rho V_{\rm
DA}(\varepsilon_{\rm DA})-\frac{1}{2}\\{V_{\rm DA}(\varepsilon_{\rm
DA})V^{\dagger}_{\rm DA}(\varepsilon_{\rm DA}),\rho\\}\Big{)}$
$\displaystyle+(n(\varepsilon_{\rm DA},T_{\rm ph})+1)\Big{(}V_{\rm
DA}(\varepsilon_{\rm DA})\rho V^{\dagger}_{\rm DA}(\varepsilon_{\rm
DA})-\frac{1}{2}\\{V^{\dagger}_{\rm DA}(\varepsilon_{\rm DA})V_{\rm
DA}(\varepsilon_{\rm DA}),\rho\\}\Big{)}\Big{]},$ (60)

where we define the operators, $\\{V^{\dagger}_{\rm D}=c^{\dagger}_{\rm
D2}c_{\rm D1}$, $V_{\rm D}=c^{\dagger}_{\rm D1}c_{\rm D2}\\}$ and
$\\{V^{\dagger}_{\rm DA}=c^{\dagger}_{\rm D2}c_{\rm A2}$, $V_{\rm
DA}=c^{\dagger}_{\rm A2}c_{\rm D2}\\}$, as the combination of system operators
that are responsible for the transition between donor ground and excited
states, as well as the transition between the excited states of the donor and
acceptor molecules, respectively. The electron transfer rates corresponding to
their respective bath are characterized by the various $\gamma$’s. Their
explicit forms in terms of the system-reservoir coupling constants can be
calculated by Fermi’s golden rule, as
$\gamma_{L(R)}=\sum_{l(r)}2\pi\hbar|g_{l(r)}|^{2}\delta\big{(}\varepsilon-\epsilon_{l(r)}\big{)}$,
$\gamma^{D}_{\rm
pht}=\sum_{k}2\pi\hbar|g^{D}_{k}|^{2}\delta\big{(}\varepsilon-\epsilon_{k}\big{)}$,
and $\gamma^{r}_{\rm
phn}=\sum_{q}2\pi\hbar|g^{r}_{q}|^{2}\delta\big{(}\varepsilon-\epsilon^{\prime}_{q}\big{)}$
where $r={\rm D}$, ${\rm DA}$. The function
$f(\varepsilon,\mu,T)=[e^{(\varepsilon-\mu)/k_{B}T}+1]^{-1}$ is the Fermi-
Dirac distribution for the left ($\rm L$) and right ($\rm R$) bath with energy
$\varepsilon$, chemical potential $\mu$ and temperature $T$ respectively.
Analogously, $n(\varepsilon,T)=[e^{\varepsilon/k_{B}T}-1]^{-1}$ is the Bose-
Einstein distribution corresponding to the photon and phonon baths with energy
$\varepsilon$ and temperature $T$ respectively. The distribution functions are
obtained by tracing over the respective bath density operator. For examples,
$\Tr_{\rm L(R)}\big{(}d^{\dagger}_{l(r)}d_{l(r)}\rho_{\rm
L(R)}\big{)}=f\big{(}\epsilon_{l(r)},\mu_{\rm L(R)},T_{\rm L(R)}\big{)}$, and
$\Tr_{\rm L(R)}\big{(}d_{l(r)}d^{\dagger}_{l(r)}\rho_{\rm
L(R)}\big{)}=1-f\big{(}\epsilon_{l(r)},\mu_{\rm L(R)},T_{\rm L(R)}\big{)}$,
where the bath operators $d^{\dagger}_{l(r)}$ and $d_{l(r)}$ obey anti-
commutation relation, whereas, $\Tr_{\rm
pht}\big{(}a^{\dagger}_{k}a_{k}\rho_{\rm
pht}\big{)}=n\big{(}\epsilon_{k},T_{\rm S}\big{)}$, $\Tr_{\rm
pht}\big{(}a_{k}a^{\dagger}_{k}\rho_{\rm
pht}\big{)}=1+n\big{(}\epsilon_{k},T_{\rm S}\big{)}$, and $\Tr_{\rm
phn}\big{(}b^{\dagger}_{q}b_{q}\rho_{\rm
phn}\big{)}=n\big{(}\epsilon^{\prime}_{q},T_{\rm ph}\big{)}$, $\Tr_{\rm
phn}\big{(}b_{q}b^{\dagger}_{q}\rho_{\rm
phn}\big{)}=1+n\big{(}\epsilon^{\prime}_{q},T_{\rm ph}\big{)}$ where the
operators $a^{\dagger}_{k}$ ($b^{\dagger}_{q}$) and $a_{k}$ ($b_{q}$) follow
commutation relations. In both cases, $k_{B}$ is the Boltzmann constant and
the transition energies driven by the $\rm L(R)$ baths are $\varepsilon_{\rm
L}=\varepsilon_{\mathbb{10}},\varepsilon_{\mathbb{43}}$, and $\varepsilon_{\rm
R}=\varepsilon_{\mathbb{30}},\varepsilon_{\mathbb{41}},\varepsilon_{\mathbb{52}}$
respectively, while $\varepsilon_{\rm
D}=\varepsilon_{\mathbb{21}},\varepsilon_{\mathbb{54}}$ (driven by both photon
and phonon bath), and $\varepsilon_{\rm DA}=\varepsilon_{\mathbb{23}}$ (driven
by phonon bath). Finally, we note that the creation and annihilation operators
can be expressed in terms of the system eigenstates in the following form

$\displaystyle c^{\dagger}_{\rm D1}$ $\displaystyle=$
$\displaystyle|\mathbb{1}\rangle\langle\mathbb{0}|+|\mathbb{4}\rangle\langle\mathbb{3}|;\quad
c_{\rm
D1}=|\mathbb{0}\rangle\langle\mathbb{1}|+|\mathbb{3}\rangle\langle\mathbb{4}|,$
$\displaystyle c^{\dagger}_{\rm D2}$ $\displaystyle=$
$\displaystyle|\mathbb{2}\rangle\langle\mathbb{0}|+|\mathbb{5}\rangle\langle\mathbb{3}|;\quad
c_{\rm
D2}=|\mathbb{0}\rangle\langle\mathbb{2}|+|\mathbb{3}\rangle\langle\mathbb{5}|,$
$\displaystyle c^{\dagger}_{\rm A2}$ $\displaystyle=$
$\displaystyle|\mathbb{3}\rangle\langle\mathbb{0}|+|\mathbb{4}\rangle\langle\mathbb{1}|+|\mathbb{5}\rangle\langle\mathbb{2}|,$
$\displaystyle c_{\rm A2}$ $\displaystyle=$
$\displaystyle|\mathbb{0}\rangle\langle\mathbb{3}|+|\mathbb{1}\rangle\langle\mathbb{4}|+|\mathbb{2}\rangle\langle\mathbb{5}|.$
(61)

So, the time evolution of the occupation probabilities which are the diagonal
elements of the reduced density matrix
$P_{\mathbb{i}}=\langle\mathbb{i}|\rho|\mathbb{i}\rangle$, can be obtained
using the Lindblad master equation [Eq. (55)] in the following way. For
example,

$\displaystyle\frac{dP_{\mathbb{0}}}{dt}=\langle\mathbb{0}|\frac{d\rho}{dt}|\mathbb{0}\rangle=\langle\mathbb{0}|\mathcal{L}_{\rm
L}[\rho]|\mathbb{0}\rangle+\langle\mathbb{0}|\mathcal{L}_{\rm
R}[\rho]|\mathbb{0}\rangle+\langle\mathbb{0}|\mathcal{L}^{\rm D}_{\rm
pht}[\rho]|\mathbb{0}\rangle+\langle\mathbb{0}|\mathcal{L}^{\rm D}_{\rm
phn}[\rho]|\mathbb{0}\rangle+\langle\mathbb{0}|\mathcal{L}^{\rm DA}_{\rm
phn}[\rho]|\mathbb{0}\rangle,$ (62)

and the first term of the R.H.S. can be calculated by using Eqs. (A) and (A)
as:

$\displaystyle\langle\mathbb{0}|\mathcal{L}_{\rm L}[\rho]|\mathbb{0}\rangle$
$\displaystyle=\langle\mathbb{0}|\sum_{\\{\varepsilon_{\rm L}\\}}\gamma_{\rm
L}\Big{[}f(\varepsilon_{\rm L},\mu_{\rm L},T_{\rm L})\Big{(}c^{\dagger}_{\rm
D1}\rho c_{\rm D1}-\frac{1}{2}\\{c_{\rm D1}c^{\dagger}_{\rm
D1},\rho\\}\Big{)}+(1-f(\varepsilon_{\rm L},\mu_{\rm L},T_{\rm
L}))\Big{(}c_{\rm D1}\rho c^{\dagger}_{\rm D1}-\frac{1}{2}\\{c^{\dagger}_{\rm
D1}c_{\rm D1},\rho\\}\Big{)}\Big{]}|\mathbb{0}\rangle$ (63)
$\displaystyle={\gamma_{\rm L}f(\varepsilon_{\mathbb{10}},\mu_{\rm L},T_{\rm
L})}\Big{(}\langle\mathbb{0}|\mathbb{1}\rangle\langle\mathbb{0}|\rho|\mathbb{0}\rangle\langle\mathbb{1}|\mathbb{0}\rangle-\frac{1}{2}\langle\mathbb{0}|\mathbb{0}\rangle\langle\mathbb{0}|\rho|\mathbb{0}\rangle-\frac{1}{2}\langle\mathbb{0}|\rho|\mathbb{0}\rangle\langle\mathbb{0}|\mathbb{0}\rangle\Big{)}$
$\displaystyle+\gamma_{\rm L}(1-f(\varepsilon_{\mathbb{10}},\mu_{\rm L},T_{\rm
L}))\Big{(}\langle\mathbb{0}|\mathbb{0}\rangle\langle\mathbb{1}|\rho|\mathbb{1}\rangle\langle\mathbb{0}|\mathbb{0}\rangle-\frac{1}{2}\langle\mathbb{0}|\mathbb{1}\rangle\langle\mathbb{1}|\rho|\mathbb{0}\rangle-\frac{1}{2}\langle\mathbb{0}|\rho|\mathbb{1}\rangle\langle\mathbb{1}|\mathbb{0}\rangle\Big{)}$
$\displaystyle=\gamma_{\rm L}(1-f(\varepsilon_{\mathbb{10}},\mu_{\rm L},T_{\rm
L}))P_{\mathbb{1}}-\gamma_{L}f(\varepsilon_{\mathbb{10}},\mu_{\rm L},T_{\rm
L})P_{\mathbb{0}}\equiv
k_{\mathbb{01}}P_{\mathbb{1}}-k_{\mathbb{10}}P_{\mathbb{0}}.$

Here we identify
$\varepsilon_{\mathbb{10}}=\varepsilon_{\mathbb{1}}-\varepsilon_{\mathbb{0}}=\varepsilon_{\rm
D1}$. Similarly, one can check that the second term reduces to

$\displaystyle\langle\mathbb{0}|\mathcal{L}_{\rm R}[\rho]|\mathbb{0}\rangle$
$\displaystyle=\gamma_{\rm R}(1-f(\varepsilon_{\mathbb{30}},\mu_{\rm R},T_{\rm
R}))\langle\mathbb{3}|\rho|\mathbb{3}\rangle-\gamma_{\rm
R}f(\varepsilon_{\mathbb{30}},\mu_{\rm R},T_{\rm
R})\langle\mathbb{0}|\rho|\mathbb{0}\rangle$ (64) $\displaystyle=\gamma_{\rm
R}(1-f(\varepsilon_{\mathbb{30}},\mu_{\rm R},T_{\rm
R}))P_{\mathbb{3}}-\gamma_{\rm R}f(\varepsilon_{\mathbb{30}},\mu_{\rm
R},T_{\rm R})P_{\mathbb{0}}\equiv
k_{\mathbb{03}}P_{\mathbb{3}}-k_{\mathbb{30}}P_{\mathbb{0}},$

and last three terms $\langle\mathbb{0}|\mathcal{L}^{\rm D}_{\rm
pht}[\rho]|\mathbb{0}\rangle$, $\langle\mathbb{0}|\mathcal{L}^{\rm D}_{\rm
phn}[\rho]|\mathbb{0}\rangle$ and $\langle\mathbb{0}|\mathcal{L}^{\rm DA}_{\rm
phn}[\rho]|\mathbb{0}\rangle$ are equal to zero. Combing these results, we
obtain from Eq. (62)

$\frac{dP_{\mathbb{0}}}{dt}=(k_{\mathbb{01}}P_{\mathbb{1}}-k_{\mathbb{10}}P_{\mathbb{0}})+(k_{\mathbb{03}}P_{\mathbb{3}}-k_{\mathbb{30}}P_{\mathbb{0}}),$
(65)

which corresponds to our main text Eq. (7) and the corresponding transition
rates are summarized in Eqs. (13)-(20) of the main text. Following a similar
procedure, we can derive all the main text Eqs. (II)-(12) for the population
dynamics. In terms of the occupation probabilities of the individual states,
phenomenological expressions for the electron currents can be written
following Nitzan et.al Einax  _et al._ (2011) as

$\displaystyle J_{\rm L}$ $\displaystyle=$
$\displaystyle[k_{\mathbb{10}}P_{0}-k_{\mathbb{01}}P_{\mathbb{1}}]+[k_{\mathbb{43}}P_{\mathbb{3}}-k_{\mathbb{34}}P_{\mathbb{4}}],$
(66) $\displaystyle J_{\rm R}$ $\displaystyle=$
$\displaystyle[k_{\mathbb{30}}P_{\mathbb{0}}-k_{\mathbb{03}}P_{\mathbb{3}}]+[k_{\mathbb{41}}P_{\mathbb{1}}-k_{\mathbb{14}}P_{\mathbb{4}}]+[k_{\mathbb{52}}P_{\mathbb{2}}-k_{\mathbb{25}}P_{\mathbb{5}}],$
(67) $\displaystyle J_{\rm S}$ $\displaystyle=$ $\displaystyle k_{\rm
r}[P_{\mathbb{1}}+P_{\mathbb{4}}]-\tilde{k}_{\rm
r}[P_{\mathbb{2}}+P_{\mathbb{5}}],$ (68) $\displaystyle J_{\rm nr}$
$\displaystyle=$ $\displaystyle k_{\rm
nr}[P_{\mathbb{1}}+P_{\mathbb{4}}]-\tilde{k}_{\rm
nr}[P_{\mathbb{2}}+P_{\mathbb{5}}],$ (69) $\displaystyle J_{\rm DA}$
$\displaystyle=$ $\displaystyle
k_{23}P_{\mathbb{3}}-k_{\mathbb{\mathbb{32}}}P_{\mathbb{2}}.$ (70)

Here $J_{\rm L}(J_{\rm R})$ is the electron currents entering (leaving) the
molecular system from (to) the electrodes, $J_{\rm S}$ and $J_{\rm nr}$ are,
respectively, the radiative (photon-induced) and nonradiative (phonon induced
losses) electron currents between ground and excited states of the donor, and
$J_{\rm DA}$ is the average current due to transfer of electrons between donor
and acceptor species. Following Ref. Sarkar and Dubi (2020), one can derive
the same from the microscopic picture, starting from the Lindblad master
equation Eq. (55). In particular, steady-state electron currents can be
obtained by using the following definitions

$\displaystyle J_{\rm L}$ $\displaystyle=$ $\displaystyle\Tr{N_{\rm
L}\mathcal{L}_{\rm L}[\rho]},$ (71) $\displaystyle J_{\rm R}$ $\displaystyle=$
$\displaystyle\Tr{N_{\rm R}\mathcal{L}_{\rm R}[\rho]},$ (72) $\displaystyle
J_{\rm S}$ $\displaystyle=$ $\displaystyle\Tr{N_{\rm D}\mathcal{L}^{\rm
D}_{\rm pht}[\rho]},$ (73) $\displaystyle J_{\rm nr}$ $\displaystyle=$
$\displaystyle\Tr{N_{\rm D}\mathcal{L}^{\rm D}_{\rm phn}[\rho]},$ (74)
$\displaystyle J_{\rm DA}$ $\displaystyle=$ $\displaystyle\Tr{N_{\rm
DA}\mathcal{L}^{\rm DA}_{\rm phn}[\rho]}.$ (75)

of current where $N_{\rm K}$ (${\rm K=L}$, ${\rm R}$, ${\rm D}$, ${\rm DA}$)
is the number operator which is defined as

$\displaystyle N_{\rm L}$ $\displaystyle=$ $\displaystyle c^{\dagger}_{\rm
D1}c_{\rm
D1}=|\mathbb{1}\rangle\langle\mathbb{1}|+|\mathbb{4}\rangle\langle\mathbb{4}|,$
$\displaystyle N_{\rm R}$ $\displaystyle=$ $\displaystyle c^{\dagger}_{\rm
A2}c_{\rm
A2}=|\mathbb{3}\rangle\langle\mathbb{3}|+|\mathbb{4}\rangle\langle\mathbb{4}|+|\mathbb{5}\rangle\langle\mathbb{5}|,$
$\displaystyle N_{\rm D}$ $\displaystyle=$ $\displaystyle V^{\dagger}_{\rm
D}V_{\rm D}=c^{\dagger}_{\rm D2}c_{\rm D1}c^{\dagger}_{\rm D1}c_{\rm
D2}=|\mathbb{2}\rangle\langle\mathbb{2}|+|\mathbb{5}\rangle\langle\mathbb{5}|,$
$\displaystyle N_{\rm DA}$ $\displaystyle=$ $\displaystyle V^{\dagger}_{\rm
DA}V_{\rm DA}=c^{\dagger}_{\rm D2}c_{\rm A2}c^{\dagger}_{\rm A2}c_{\rm
D2}=|\mathbb{2}\rangle\langle\mathbb{2}|.$ (76)