The Molecular Photo-Cell: Quantum Transport and Energy Conversion at Strong Non-Equilibrium

The molecular photo-cell is a single molecular donor-acceptor complex attached to electrodes and subject to external illumination. Besides the obvious relevance to molecular photo-voltaics, the molecular photo-cell is of interest being a paradigmatic example for a system that inherently operates in out-of-equilibrium conditions and typically far from the linear response regime. Moreover, this system includes electrons, phonons and photons, and environments which induce coherent and incoherent processes, making it a challenging system to address theoretically. Here, using an open quantum systems approach, we analyze the non-equilibrium transport properties and energy conversion performance of a molecular photo-cell, including both coherent and incoherent processes and treating electrons, photons, and phonons on an equal footing. We find that both the non-equilibrium conditions and decoherence play a crucial role in determining the performance of the photovoltaic conversion and the optimal energy configuration of the molecular system.

complex, a phonon bath (at ambient temperature) and a photon bath (at the solar temperature). Non-equilibrium is induced by two sources, namely the temperature difference between the incoming photons (originating from the sun) and the phonons with ambient temperature, and the bias voltage between the electrodes.
Using this formalism, we calculate the non-equilibrium densities of electrons, photons and phonons, the electric current and power output, and the thermodynamic efficiency at maximal power. We find that under certain conditions, the distribution functions for the phonons and photons can be very different from the equilibrium distributions, and therefore approximating the system as close to equilibrium is not a valid approximation. We then study the signature of non-equilibrium on the energy conversion efficiency of the system.
In addition to be able to include strong non-equilibrium effects and to account for all constituents on equal footing, our formalism allows one to introduce effects of dephasing and decoherence in a simple and physically transparent way. We study how the loss of quantum coherence affects the efficiency, and show that classical electron transfer (as opposed to coherent electron propagation) enhances the efficiency.
Results: efficiency for the coherent system far from equilibrium The system under consideration is composed of a molecular junction, in which a donor-acceptor (D-A) complex is placed between two metallic electrodes ( Fig. 1) 4 . This is an idealization of an envisioned future single-layer molecular photovoltaic cells, where a selfassembled layer of D-A pairs is placed on a conducting substrate, and is covered by a top transparent electrode. For the donor (D), we consider the highest occupied molecular orbital (HOMO) with energy E D,1 and lowest unoccupied molecular orbital (LUMO) with energy E D,2 . For the acceptor (A), we consider only the LUMO with energy E A (the energy of the A's HOMO is much lower and is unaffected by any dynamics in the photocell 4 ). The reader is referred to the Methods section and supplementary materials for a detailed description of the model and calculation. Recent advances in the experimental ability to measure photo-conductivity and PV conversion in single-molecule junctions [11][12][13][14][15] make our theoretical model experimentally relevant.
We begin by examining the electron, photon and phonon densities at zero bias voltage. Note that the system is still out of equilibrium due to the temperature difference between the photon and phonon baths. In Fig. 2(a) steady state averages of the photon density n pht~npht , the phonon density n phn~nphn , and the electronic occupation of the D-LUMO, n D,LUMO~nD,LUMO h i , are plotted as a function of e-pht coupling l e-pht , wheren pht ,n phn andn D,LUMO are particle number operators of the photon, the phonon, and the D,LUMO, respectively. Here, AE?ae represents the steady state average, i.e., AE?ae ; tr(?r SS ) with the steady state density operator r SS (see supplementary material). We set the orbital energies to be E D,1~{ 0:1 eV, E D,2~1 :4eV, and E A~1 :25eV 4 . The e -phn interaction is fixed at l e-phn 5 0.1 eV (dashed vertical line in Fig. 2(a)). Three distinct regimes are observed: (i) At small e -pht coupling, l e{pht =l e{phn , the system is close to equilibrium, and the photon occupation is defined by the solar temperature T s (dotted line). The phonons are excited according to the ambient temperature (which is very small compared to v 0 , and consequently the phonon occupation is very small), and the occupation of the D-LUMO level is also very small. (ii) As the e -pht and e -phn interactions become comparable, energy is transferred from the photons to the phonons, mediated by excitation of electrons from the D-HOMO to the D-LUMO level. As a result, the D-LUMO and phonon occupations increase, while the photon occupation decreases. Alternatively, this situation can be described in terms of heating (although the notion of   temperature is not applicable out of equilibrium, it is still useful to think in terms of an effective temperature): the junction is locally (and efficiently) heated by the photons. This heat is transferred to the phonons, resulting in an elevated effective phonon temperature and, consequently, enhanced phonon occupation. (iii) When epht coupling is large, l epht ?l ephn , there is no longer efficient transfer of heat to the phonons. However, the D-LUMO occupation continues to rise, due to energy pumping from the photons to the D-LUMO.
An important measure of the operational efficiency of the molecular PV cell is the efficiency at maximal power, g mx , defined as the ratio between the cell's maximal output power P out and the corresponding input power P in supplied by the photons 16,17 . The output power is given by P out 5 J 3 V, where J is the particle current through the system, and the input power is calculated in a similar way 4 from , which is related to the photon-induced part of the particle current. The inset of Fig. 2(b) shows a typical J 2 V characteristics (blue) and the bias voltage V dependence of the output power (purple). As seen in the figure, there exists a bias voltage V mx which gives the maximum efficiency. In Fig. 2(b), the efficiency at V mx , i.e., the efficiency at maximum power, is plotted as a function of the epht coupling (we use the same parameters as in Fig. 2(a)). For very small epht coupling (l epht , 0.002 eV), the efficiency is very small, and grows linearly with l epht . In this regime, the time it takes for an electron to absorb a photon is larger than the time the photons spent in the cell (defined by c pht n B (T s ) which translates to ,0.002 eV), and so the photon absorption is very small leading to poor efficiency. This regime is followed by a plateau regime, where any energy transferred from the photons to the electrons is quickly dissipated by phonons and is not converted into electrical power. There is thus little change in the efficiency, as long as the rate of photon absorption is smaller than the electron-photon-relaxation time. Only when the e -pht coupling reaches the e -phn coupling l ephn 5 0.1 eV the efficiency begins to increase: in this regime, energy is transferred to the electrons by the photons faster than can be dissipated by the phonons, and as a result an increasing amount of this energy is transferred into electronic power, resulting in a rise of efficiency.
It is important to note that the results described in Fig. 2, especially in the region where the electron-phonon and electronphoton couplings are of the same order, cannot be obtained by assuming equilibrium distributions for the phonons and photons, and this situation is described here for the first time. Since the strength of the electron-photon and electron-phonon interactions of future realistic devices are unknown, a situation where they are of similar magnitude may occur, in which case the system dynamics cannot be described as close to equilibrium, and the full non-equilibrium dynamics need to be taken into account.
To further demonstrate the power of this method, we next discuss the effect of Coulomb interactions on the efficiency. In the Hamiltonian of Eq. (1), the A-LUMO energy, E A , already includes the Coulomb repulsion energy on the acceptor 4 . In excitonic systems, the Coulomb interaction is typically considered through the ''exciton binding energy'', defined by the Coulomb interaction term between an electron at the D-LUMO and a hole in the D-HOMO. While in methods such as non-equilibrium Green's function adding Coulomb interaction requires substantial effort, the present method does not require either additional technical complexity or additional computational power to account for any Coulomb interaction effects. In Fig. 3 we show the efficiency at maximum power g mx as a function of the exciton Coulomb energy U (which can be estimated from, e.g. density-functional calculations). We set l epht 5 0.1 eV and l ephn 5 0.2 eV. We find an almost linear decrease in the efficiency, with a reduction of ,15% for U 5 0.2 eV.
Results: the role of decoherence The next question we wish to address is the extent to which the quantum nature of the system affects the PV conversion efficiency, a question which is beyond the reach of the formalism presented in Refs. 4, 5. The formalism we present here allows us to access, in addition to fully quantum-coherent processes described above, also incoherent processes. To do so, we define a variable j such that 0 # j # 1, and define t mx 5 0.05 eV and c mx 5 10 12 s 21 (as in Ref. 4). We now parameterize t D-A and c D-A with j, t D-A 5 t mx (1 2 2(j 2 0.5)H(j 2 0.5)), c D-A 5 c mx (1 2 2(0.5 2 j)H(0.5 2 j)) (H(j) is the Heaviside unit stepfunction). This parameterization is shown on the right inset of Fig. 4, and is constructed such that for j 5 0 the system is fully coherent, for j 5 0.5 the system is mixed (both quantum and classical processes), and for j 5 1 the system is fully incoherent, so the range 0 , j , 1 interpolates between all three cases.
In Fig. 4 we plot the efficiency at maximum power g mx as a function of the position of the A-LUMO, E A and the parameter j. We set c D-A 5 10 12 s 21 as in Ref. 4. We find that the quantum coherence or classical decoherence (parameterized by j) has a profound effect on the efficiency of the molecular PV-cell in two important aspects.
First, the optimal position of the A-LUMO energy differs according to the nature of the transition under consideration: quantum (coherent), both quantum and classical, or classical D-A transitions (solid lines in Fig. 4). For the last case (t D-A 5 0, c D-A 5 10 12 s 21 ), we find that E A is optimal at E A *1:25 eV, verifying the result of Ref. 4. When quantum correlations are added (t D-A 5 0.05 eV, c D-A 5 10 12 s 21 ), two peaks emerge at E A *0:9, 1.3 eV, and a lower peak emerges at E A *1:6 eV. For a system with only quantum transitions (t D-A 5 0.05 eV, c D-A 5 0, enlarged in the back inset in Fig. 4), the lower peak vanishes, and the optimal LUMO positions are at 1.2 eV and 1.6 eV. Thus, in the design of optimal molecular PV cells, it is important to take into account the quantum nature. Second, as can be clearly seen in Fig. 4, the addition of classical D-A transitions increases the efficiency substantially by more than an order of magnitude. This finding is surprising, since one would expect that incoherent (and dissipative) transitions would lead to a decrease in the efficiency. To understand the origin of this effect, we performed time-dependent calculations (not shown) for a system composed of D-LUMO, acceptor, and the coupling with right electrode, and found that for the coherent case, an electron that is excited to the D-LUMO coherently oscillates between the D-and A-LUMO, while for the incoherent case the electron decays from the D-to the A-LUMO exponentially (and its return rate is exponentially small). This implies that in the coherent case the electron spends much more time in the D-LUMO than in the incoherent case, before transferring to the right electrode. Since the electron can decay back to the D-HOMO (emitting a phonon) only directly from the D-LUMO, the longer it spends on the D-LUMO, the higher the probability for nonradiative decay back to the D-HOMO, leading to a decrease in efficiency. This phenomena is similar to dephasing-assisted transport conjectured to occur in biological systems [18][19][20][21] , but here is the first time it is discussed and demonstrated in the context of molecular PV cells. Since in realistic single-molecule junctions both coherent and incoherent effects may be important, they must be included in a theoretical description of the system.
To illustrate this connection between dynamics and efficiency, we address the relaxation dynamics of the donor-acceptor system. Considering only the D-A LUMOs and the right electrode (without the photons and phonons), we construct the LindbladianM-matrix of Eq. (2) (at zero bias), by constructing a vector form for the Lindblad equation of Eq. (2), _ r~{Mr. TheM-matrix has a zero eigenvalue, which defines the steady state. The (real part of the) rest of the eigenvalues define the relaxation rates towards the steady state. The minimal rate C min (i.e., eigenvalues ofM with smallest real part which, we numerically check, is non-zero) represents the longest relaxation time for the system to reach the steady state from any general state (see supplementary material).
In Fig. 5, we plot the decay rate C min of the model which only contains D-A LUMOs and the right electrode (without the photons and phonons) as a function of the position of the A-LUMO E A for the three cases of fully coherent (t D- acceptor electron transfer processes. As can be seen, the decay rate is much smaller for the fully coherent case, indicating a longer relaxation time. This is in line with the observation above that slower relaxation dynamics lead to lower efficiency. In addition, we point out that the relaxation times of the D-A LUMOs and the right electrode (which are numerically much easier to calculate than the efficiency of the full system including D-HOMO, photons, phonons, and related Lindblad dissipators) serve as an indicator for the efficiency, even though they do not capture the fine details required for an optimal design of the system (see Fig. 4(a) and (b)).
To further examine possible effects of coherence on the efficiency of the molecular PV cell, we study a system where the donor has two degenerate D-LUMO levels which have been introduced experimentally 22,23 . Here, we study a simplified system (schematically depicted on the right side of Fig. 6), in which the photons and the phonons excite electrons with equal amplitudes from the D-HOMO to the two D-LUMO levels. Each of the levels is coupled to the A-LUMO with the same hopping amplitude t D-A , and they are coupled to each other with a complex hopping amplitude he 2ipw . In Fig. 6, the efficiency at maximum power g mx is plotted as a function of the inter-LUMO coupling h (solid line) and the phase w (dashed line) for c D-A 5 0, i.e., no incoherent D-A transfer (l ephn 5 l epht 5 0.1 eV). We find that while h has little effect on the efficiency, the phase w has a significant effect (increasing the efficiency by up to ,15%).  Surprisingly, we also find that this quantum interference effect persists even when incoherent D-A transfer is included (c D-A 5 10 12 s 21 as in Fig. 4), and an substantial increase of g mx , 30% is observed by varying w (dotted line).

Conclusions
In Summary, we have proposed a novel formalism to study nonequilibrium quantum transport in molecular junctions, and applied it to investigate a minimal model of PV energy conversion in ideal, single-molecule PV cells. The results shown above indicate that quantum coherence effects are important in determining the nonequilibrium energy conversion performance of molecular PV cells. The formalism presented here sets the stage for a fully coherent quantum mechanical calculation of energy conversion in more realistic models for molecular PV cells and can be directly linked to quantum chemistry methods (such as density-functional theory). The progress in the experimental ability to measure photo-conductivity and PV conversion in single-molecule junctions [11][12][13][14][15] allows one to envision real PV devices composed of a single molecular junction or a molecular monolayer, making our theoretical model experimentally relevant. Furthermore, our method can include both coherent and incoherent effects, making it a useful tool in the study of other energy conversion processes such as photosynthesis, where both classical and quantum processes take place [18][19][20][21][24][25][26][27][28][29][30][31][32][33][34][35][36][37][38] , or other chemical and photo-chemical processes 39,40 .
The essence of the Lindblad approach is that instead of describing the environment by encoding it into a self-energy (as is done in the non-equilibrium Green's function approach 42 ) the environment is characterized by its action on the system. This action is mapped onto so-called LindbladV-operators, which describe incoherent transitions of the system elements due to the presence of an environment. The Lindblad equation was recently employed to address various aspects of electron transport [43][44][45][46][47][48][49] , yet in these studies the interaction with an environment was limited to electrons only, and the non-equilibrium dynamics of other constituents (i.e. phonons or photons) was not considered.
We assume that the left electrode is coupled only to the D-HOMO and that the right electrode is coupled only to the A-LUMO 4 , as in Fig. 1. The correspondingVoperators are then 44,48,49 : where c L,R are electron transfer rates to the left and right electrodes, T is the ambient temperature (we take T 5 300 K), m L 5 0 is the left-electrode chemical potential, m R 5 V is the right electrode chemical potential, V is the bias voltage, where c pht , c phn are photon and phonon relaxation rates (set to c pht 5 c phn 5 10 12 s 21 ), T s , 5700 K is the solar temperature, is the Bose-Einstein distribution, andñ B~1 zn B . Once theV-operators are set, the dynamics are determined by the propagation of the density matrix via the Lindblad equation. We numerically study the steady-state density matrix r SS with truncated bosonic space containing n excitations and find that all the physical properties converge at n 5 6. Therefore, we demonstrate our results with n 5 6 (for comparison, in Ref. 4 the photons and phonons were treated in a nonself-consistent mean-field approximation). The expression for the current is obtained