The existence of light induced coherences in biological systems and their effects on efficiency has been a topic of considerable interest. [1] [2] [3] Most of these studies have focused on excitation with ultra-fast coherent light sources, 4,5 which is not typical of natural conditions. Rather, incident sunlight is incoherent and provides a constant intensity over the excitation period. 6, 7 Studies of energy transfer in donor-acceptor(DA) systems under incoherent radiation [8] [9] [10] [11] [12] have provided evidence that it is possible for long time coherences to occur under the influence of a natural incoherent pumping source when the system is coupled to a secondary bath.
There have been a variety of approaches to describe the effects of natural radiation on quantum subsystems. Many of the early studies focused on the incoherent light-matter interactions using Bloch-Redfield theory, although a white noise description 13 has also been developed. These results were then extended to include effects of trapping and bath dephasing through the Lindblad formalism, which treats the incoherent light in the secular approximation. 14 However, describing light-matter interactions at the secular level treats the populations and coherences on different footings, and important effects such as Fano coherences cannot be observed. 15 A recent study 9 remedied this situation by adopting a hybrid approach that combines the non-secular description of the Bloch-Redfield treatment with a Lindblad description for the rest of the environment. With the use of this formalism analytic results were obtained for a symmetric DA system in the steady state limit. As insightful as this study is it focused on the idealized, symmetric case which is not typical of what is found in nature. It is to this end that we study the interplay of weak incoherent light and a phonon bath in an asymmetric DA model, building upon the previous master equation approaches, to create a general description that incorporates the important features of the light, the effects of a thermal bath, energy relaxation, and the role of an energy bias in the system. The most obvious goal for this study is physical insight, since one would not expect most naturally occurring DA systems to be symmetric.
The layout of the paper is as follows. We first discuss, in Section II, a generalization to the master equations that incorporates an energy gap in the donor-acceptor system Hamiltonian as well as the metrics used to quantify the energy transfer. Two issues are emphasized: the relation of the trapping efficiency to the flux in the steady state limit, and the nature of the flux as an interference phenomenon. Numerical results obtained from these master equations by evolving them to the naturally relevant steady state limit are presented and discussed in Section III. Section IV summarizes our findings.
A. Model system for photosynthetic energy transfer
The model system used in this study is that presented in Ref. 9, i.e., a simple donoracceptor (DA) exciton system
with on-site energies k and site hopping coefficient J. However, unlike Ref. 9 the system here is asymmetric, i.e. 1 = 2 . The dipoleμ = µ 1 |g 1 1| is assumed to be aligned on the donor site, see figure 1. The DA system interacts with an incoherent radiation field and is coupled to a reaction center. The description of the reduced density matrix is that of a hybrid Bloch-Redfield-Lindblad approach, meaning that the light-matter interactions will be treated using standard non-secular Markovian 16 optical master equations 12, [17] [18] [19] and the effects of the bath, recombination of excitons, and the transfer of energy to the reaction center is treated as a dephasing of the system through the use of Lindblad operators. 14, 20 Under these assumptions the density matrix evolves according tȯ
where
The effects of the DA system coupled to the surrounding environment are given by
where [A, B] + is the anticommutator of A and B and γ d is the dephasing rate. The term
accounts for the recombination of the excitons, and transfer from the acceptor site to the reaction center |RC is described by
The Bloch-Redfield equations for the light-matter interaction are most easily described in the energy-eigenbasis 9,21 and are given by the following set of elements of L rad ρ:
where ∆ = ( 1 − 2 ) 2 + 4J 2 is the excitonic splitting, r e ± are the pumping rates of the energy eigenstates |e ± and γ e ± are incoherent emission rates. Table I summarizes the relevant parameters introduced, and the density matrix elements {ρ e + e + , ρ e − e − , ρ R e + e − , ρ I e + e − } describe the populations of the DA energy-eigenstates and the real and imaginary parts of the coherences between these eigenstates, respectively.
We seek working master equations in the eigenstate basis, where the light-matter interaction is most easily described. To do this the matrices that appear in the Lindblad operators need to be rotated into the eigenbasis via the unitary matrix 22, 23 
where θ is the diabatic mixing angle given by
where 0 ≤ θ ≤ π 2 . Four cases are relevant: local (in the site basis of the DA dimer) dephasing with either localized or delocalized trapping and global (in the energy-eigenbasis of the dimer) dephasing with either type of trapping. Specifically, diverse results can be obtained from each set of conditions, which are distinct from the often-studied localized trapping 3,24 and localized dephasing case. The two trapping cases considered correspond to two limiting regimes. The local trap arises when the dimer and the trap are close together allowing for spatial coherences to form between them. The delocalized trap corresponds to the opposite situation, in which there is a large spatial separation of the dimer and the trap that limits the formation of these coherences.
Consider first when dephasing is applied locally 14 (in the site basis of the DA dimer).
Combining the results of Eqs (3)-(7) with Eq. (2) the density matrix in this scenario evolves according to the following master equations:
where the ground state population ρ gg remains close to unity and where any eigenstate population density matrix elements multiplied by the pumping rates r e ± or the spontaneous emission rates γ e ± are negligible. Both of these approximations are justified due to the extremely weak interactions of the incoherent light with the system. κ is a parameter that is set to unity to describe localized trapping or to zero for delocalized trapping. Physically, this can be viewed as an alignment factor between the eigenstates and the reaction center.
[Note the extension from Ref. 9 where the DA system in that case was symmetric giving r e + = r e − = r. Here the rotation from the site basis to the eigenbasis modifies the pumping rates to r e + = 2r sin 2 θ and r e − = 2r cos 2 θ.] The second major case considered is when dephasing is applied globally (in the eigenbasis of the dimer). The system evolves similarly to Eq. (9) except now the dephasing rate affects both of the coherences equally and is absent in the evolution of the populations. 25 The master equations to describe this situation are given bẏ ρ e + e + = −[2Γ + 2Γ RC cos 2 (θ)]ρ e + e + + r e + ρ gg + κΓ RC sin(2θ)ρ R e + e − ,
The relationships between the density matrix elements in the energy-eigenbasis and in the site basis are given by ρ 11 = sin 2 θρ e + e + + cos 2 θρ e − e − + sin(2θ)ρ R e + e − ,
For the symmetric case 1 = 2 the transformation reduces to:
and the master equations given in Eqs. (9) and (10) which has been corrected here in Eq. (12) .] It is interesting to note that the imaginary part of the coherence in the dimer system only changes in sign between the site and eigenstate basis and is independent of the rotation angle. We emphasize that this result is particular to dimer models and would not hold for more complex systems. We also take note that the master equations presented above only depend on 1 and 2 through ∆ which contains the square of the energy difference. This has the effect that only | 1 − 2 | needs to be considered.
i.e. the system is insensitive to the sign of the energy difference.
The energy transfer efficiency is defined by the amount of population in the site neighboring the reaction center 10 appropriately weighted by the transfer rate to the reaction center and by the initial excitation rate. For localized trapping conditions it is given by
Equation (13) has a form suggestive of a coherent control scenario 26 for unitary dynamics, i.e.
two direct terms dependent on populations, and an interference-like term that is dependent on coherences. There is, however, a significant difference, emphasized later below. That is, in the unitary case the magnitude of the population terms do not depend on the coherences.
This is not the case here, where the time dependence of the population terms is tied to the coherences via Eqs. (9) and (10) . With this in mind η loc may be rewritten as
where η direct = Γ RC r (cos 2 θρ e + e + + sin 2 θρ e − e − ) and η inter = −2 Γ RC r sin θ cos θρ R e + e − . For the case of delocalized conditions the efficiency is given by
which reflects the fact that the coherences between different eigenstates have disappeared in this limiting case because of decoherence induced by delocalized trapping. Here, transfer happens directly between system eigenstates and the reaction center. This corresponds to an incoherent or Forster mechanism. [27] [28] [29] [30] [31] It is enlightening to explore the role of coherences in η loc in greater detail. For the pure dephasing model studied here the flux F 12 between the two sites is related to the imaginary component of the coherences (see Appendix A) as
where the second equality follows from the third line of Eq. (11), which only holds for the two level systems.
Equation (16) indicates that in an open quantum system subject to dissipation arising from pure dephasing conditions the flux between adjacent sites is entirely due to coherences.
Interpreted in terms of ρ 12 , this is a manifestation of spatial coherence between the two sites.
That is, zero overlap of the spatial contributions of the two sites would yield ρ 12 = 0. A non-zero value indicates overlap. In addition, in the case of two sites, this is directly related to the ability of the energy eigenstates |e ± to span both sites, reflected in ρ I e + e − . Further, ρ I e + e − can be viewed as arising from interfering pathways, that is, it arises because |e + can emit into the background bath and |e − can absorb from the bath, creating two pathways to the bath. (Similarly, |e − can emit and |e + can absorb).
As shown in Appendix B Eqs. (13) and (16) can be related in the steady state limit through the use of the general relationship
which is valid regardless of the trapping conditions, the basis the dephasing is applied in, the pumping rate, and the asymmetry of the dimer. Notice that Eq. (17) implies that one of the coherence terms in the energy-eigenbasis is linearly dependent on the other. Using
Eq. (17) with Eqs. (13) and (16) the localized trapping efficiency can be rewritten as
which shows that the trapping efficiency is directly related to the site-to-site flux.
The character of Eq. (18) is now clear. The first two terms are direct terms and the last term is an interference term. The magnitude of the interference contribution is modulated by the sum of the bath-induced dissipation rates, namely (2Γ + Γ RC + 2γ d ), and system characteristics (2J∆). Additional system dependence enters through the sin(2θ) term and the J contained within F 12 . The flux contains the explicit coherence contribution. A further comment should be made about the form of Eq. (18) . While it is tempting to refer to the first two terms as a classical contribution since these terms involve the populations, and the third term as a quantum contribution since the flux is related to the coherences, this attribution needs to be qualified since the steady state limit reached depends on the coupling of the population to the coherences, i.e., how Eqs. (9) and (10) are constructed. Thus, there is always an implicit dependence on the coherences in the first two terms in Eq. (18) and in Eq. (15) . This implicit dependence will be seen shortly for a series of limiting cases and explored further in the next section.
While it is possible to obtain analytic steady state solutions for Eqs. (9) and (10), the general forms are too cumbersome to extract physically meaningful results so we focus on special cases to obtain insight into the behavior of the efficiency expressions. The regime of easily interpretable results is that in which Γ RC γ d (for details see Appendix C). That is,, for Γ RC γ d :
(1) The first case of interest is that of a symmetric dimer with local dephasing and local trapping. In this limit the efficiency can be shown to be
= 1.
While the final result itself is uninteresting, as the dependence of η loc on any of the parameters has vanished, the second line in Eq. (19) is enlightening as it shows that the flux contribution to the population term is exactly canceled out by the flux contribution from the interference term. Even though this is a special case it highlights the dependence of the steady state population values on the flux/coherences and shows that the one cannot just look at η inter to determine the role that the flux plays.
(2) For the case of a general dimer with local dephasing and delocalized trapping the efficiency is given by
i.e., the efficiency is also unity and is also independent of | 1 − 2 |. In this limit the dependence of the efficiency on γ d cannot be inferred and Eq. (20) should be treated as an upper limit as it is an idealized case where the bath has minimal effects and the trapping rate is dominant.
(3) Considering again the symmetric dimer, but now with global dephasing and local trapping, the efficiency can be expressed as
revealing that, in a similar fashion to Eq. (19) , the flux term (which is equal to 2r) is canceled out and does not contribute in this limiting case.
(4) The final case is that of a general dimer with global dephasing and delocalized trapping, and the efficiency is given by
suggesting, along with Eq. (20) , that in the Forster trapping regime the details of the dimer, other than the energy bias, doesn't have a large impact in the energy transfer efficiency when the trapping rate is faster than the bath dephasing rate.
A major limitation of the above results is that they do not reveal the dependence of the efficiency on the dephasing rate γ d or the flux. They are meant to serve as guides to understand some of the complex behavior within this model. Hence, we now explore the stationary solutions of Eqs. (9) and (10) numerically to study the dependence of Eqs. (18) and (15) outside of the parameter regimes where simple analytical results are obtainable.
To obtain the stationary limit values of the density matrix either Eq. (9) or (10) (depending on the context) was integrated using the RK4 method with initial conditions Typical timescales for the steady state limit to be reached are within 10ps−100ps depending on the dephasing rate.
Consider first the localized trapping condition case with dephasing applied locally. Figure   2 shows a contour plot of the steady state efficiency as a function of γ d and | 1 − 2 |. The value of 1 − 2 = 0 corresponds to the symmetric dimer (i.e. ∆ = 2J) and as | 1 − 2 | increases the two sites become more asymmetric. The efficiency shows a monotonic decrease as the dephasing rate increases for small | 1 − 2 |. When the absolute value of the DA energy difference reaches a value of 1.3 ps −1 the efficiency begins to show a turnover. This turnover in the efficiency is a unique behavior found only in asymmetric systems. An example of this turnover is presented in Fig. 3 , which also shows η direct and η inter , to be discussed below. A similar turnover is seen in the imaginary part of the coherences. In fact, the enhancement of the efficiency appears when there is a minimum in ρ I e + e − , corresponding to a maximum in the flux (see Eq. (16)). Despite the difference in scale between Fig. 3 (a) and Fig. 3(d) the flux and η loc appear to be proportional to each other. This maximum occurs due to an enchantment in the site-to-site flux arising from the environmental dephasing rate γ d
and modulates the rate at which the interference grows. Prior to the maximum, η inter grows slower than η direct and is accompanied by a rising overall efficiency. After the maximum, η inter now grows faster approaching the value of η direct allowing the two terms to become competitive, leading to a vanishing efficiency. These effects are subtle, but demonstrative.
If the dephasing is too slow the difference between the two contributions to the efficiency stays roughly constant due to the flux staying constant, and if the dephasing occurs too quickly the flux vanishes due to the coherence being destroyed between the two eigenstates which causes η inter to approach −η direct . An important observation is that this turnover behavior is only found for small values of J for the dephasing rates explored here. For example, Fig. 4 shows that with increased J the efficiency shows an overall increase in magnitude relative to Figure 2 . In this regime the system is so strongly coupled that dephasing has a much smaller effect on the flux as compared to the small J case, and the rate at which the efficiency changes is monotonic, i.e. no maximum occurs.
Focusing now on the dimer with a delocalized trapping mechanism and dephasing applied locally we see from Fig. 5 a notable difference from the previous case in that the long term efficiency for small splittings under this mechanism is almost completely insensitive to the magnitude of the dephasing rate and only depends on the splitting of the two eigenstates.
This agrees with the result in Eq. (20) which states that the efficiency should be at or near unity. Again, around | 1 − 2 | = 1.3 ps −1 the same enhancement in the efficiency is observed for similar values of γ d . Figure 6 shows the same relation of the imaginary part of the coherences and the transfer efficiency as was seen in Fig. 3 . Here, the turnover is manifest in an interplay of the two contributions in η direct , for smaller values of γ d the cos 2 θρ e + e + contribution grows at faster rate than the sin 2 θρ e − e − contribution. The general relationship of the efficiency and the flux/coherences is not as clearcut for this case as it was for the previous one. That is, as mentioned previously, although Eq. (15) doesn't explicitly depend on the flux/coherences their presence is still seen in the modulation of the efficiency through the individual direct contributions. This implicit dependence is stronger in this case than the previous one and could be due to the fact that the κΓ RC sin(2θ) terms are absent in Eq.
(9) when κ = 0, which could allow for more coherent behavior to manifest in η direct due to the absence of this dissipative term.
Turning attention to localized trapping conditions when dephasing is applied globally, not be revealed if one focused on symmetric or nearly symmetric dimers as the efficiency is nearly independent of γ d for small site energy differences in both cases. This suggests that the trapping mechanism can play a significant role in modulation of the efficiency as the asymmetry of the sites increases.
The final case of global dephasing applied to a dimer with a delocalized trapping mechanism is shown in Fig. 8 . Similar to the case of the dimer with a delocalized trapping mechanism and local dephasing (Fig. 5 ) the efficiency depends on the DA energy difference more than any other parameter. However, a unique feature of this case is that it appears that the efficiency only depends on the site energy difference. This behavior is in excellent agreement within a factor of proportionality to Eq. (22) , which is rather surprising given that Eq. (22) was derived for a special case. Here, then, the efficiency for this case is completely determined by the site energy difference. 
We have generalized a set of master equations describing a dimer coupled to an environment comprised of a phonon bath and natural incoherent light to consider the effects of asymmetry on the energy transfer efficiency. It was seen that asymmetry reduces the overall efficiency as compared to the fully symmetric case regardless of trapping mechanism or dephasing model. Further, the results only depend upon | 1 − 2 | as opposed to 1 − 2 .
We have also derived a general relationship between the ratio of the real and imaginary part of the coherences in the steady state limit. This result has allowed us to directly connect a contribution of η loc to the site-to-site flux and to the coherences. For a series of limiting cases we were able to relate the flux to the steady state population values, although in these cases the flux dependence was canceled out. We then went on to numerically demonstrate the role flux and bath induced dephasing play in energy transfer process for more general parameter regimes. Our results demonstrate that within the localized trapping models the coherence is closely related to the steady state efficiency. This suggests that if one finds a way to control and tune the coherences that the efficiency can likewise be controlled.
Under both local dephasing models it was found that a turnover behavior appears in the efficiency as a function of the dephasing rate γ d as the difference in site energies of the dimer become larger, regardless of the assumed trapping mechanism. This effect can only occur in asymmetric systems and reveals that in some cases there is an optimal amount of dephasing required to reach maximum efficiency for a given splitting. The results suggest that a turnover occurs regardless of the distance between the dimer and the reaction center.
Rather, the importance lies in whether the bath couples to the individual dimer sites or to the eigenstates of the dimer. Further, we found that whenever this turnover occurred it was the coherences play a role in the steady state efficiency. A summary of results is presented in Table II .
One key issue that needs to be overcome in future modeling is the use of pure dephasing rates. The Lindblad operators from which these dephasing rates stem generally have little physical meaning as they have no origin from a proper Hamiltonian and they also introduce many phenomenological parameters. Future work in this area must include a model
Hamiltonian from which specific information can be derived. To address these shortcomings more rigorous modeling would need to be done to develop a deeper understanding of energy transfer and the role of asymmetry and coherences. Such work is ongoing in our laboratory.
The results contained here demonstrate that within the localized trapping models the coherence is closely related to the steady state efficiency. As a consequence, the efficiency can be controlled by tuning the coherences, This is a direction that warrants further investigation. 
where 
and
Note that the interaction Hamiltonian is proportional to the occupation on the sites i.e.
H I ∝ |n n|. This is the pure dephasing condition assumed, for example, in the Lindblad operators.
The rate of change of the population in each site is given by
Since the total population is conserved i.e. n P n = 1 this implies the continuity condition
where F mn is the flux from site m to site n. Comparing Eq. A5 and Eq. A6 we see that
which states that the flux from one site to another is solely dependent on the coherences between them. This statement is only true when H I ∝ |n n|. To demonstrate this point consider changing the interaction Hamiltonian to
Note that the interaction Hamiltonian is still dependent only on system operators but it no longer dependent on the population of each site. Repeating the calculation above with this new interaction form yields
Notice the additional term that has appeared which is a direct result of new form of H I .
It is also worth noting that if the form of the interaction is taken to be
Then the flux will still be given by Eq. (A7) due to the fact that T r B ρ eq (a † k + a k ) = 0 for a harmonic bath.
Appendix B: Derivation of Eq. (17) We begin by rewriting Eq. (9) in matrix form aṡ
where ρ(t) = [ρ e + e + (t), ρ e − e − (t), ρ R e + e − (t), ρ I e + e − (t)], B = [r e + , r e − , √ r e + r e − , 0], and
with Γ + = 2Γ+2Γ RC cos 2 θ, Γ − = 2Γ+2Γ RC sin 2 θ, Γ * RC = kΓ RC sin(2θ), γ s = 2Γ+Γ RC +2γ d , γ * d = γ d sin(2θ) cos(2θ), and d = γ d sin 2 (2θ). The steady state solutions can be obtained from
which yields for the coherences
and
which is Eq. (17) . One can follow similar steps to verify that if Eq. (10) were used instead of Eq. (9) then the end result would be the same, showing that the result is independent of how the dephasing is applied. Equation (17) is also independent of κ, meaning that the trapping conditions also do not change the ratio of the real and imaginary parts of the coherences in the steady state limit.
Appendix C: Efficiency expressions in the limit of Γ RC γ d
Equation (19): The steady state solutions of Eqs. (9) and (10) 