Landscape, kinetics, paths and statistics of curl flux, coherence, entanglement and energy transfer in non-equilibrium quantum systems

We develop a population and flux landscape theory for general non-equilibrium quantum systems. We illustrate our theory by modelling the quantum transport of donor-acceptor energy transfer. We find two driving forces for the non-equilibrium quantum dynamics. The symmetric part of the driving force corresponds to the population landscape contribution which mainly governs the equilibrium part of dynamics while the anti-symmetric part of the driving force generates the non-equilibrium curl quantum flux which leads to the detailed-balance-breaking and time-irreversibility. The multi-loop structure of the flux emerges forms the flux-landscape. We study the trend of changes in population and flux-landscape with respect to the voltage (temperature difference induced by environments) and electronic coupling. Improving the voltage and electronic coupling in general facilitates the quantum transport by reducing the population landscape barriers between major states and increasing the mean value of the flux. A limit-cycle mode emerges when the underlying flux-landscape becomes funnelled with a significant gap between the largest flux loop and the rest of them. On the kinetic level, we find that multiple kinetic paths between quantum states emerge and illustrate the interference effects. The degree of interference is determined by the landscape and flux. Furthermore, we quantify kinetic rate which strongly correlates with the population landscape and flux. For quantum transport, we demonstrate that as the coherence or the quantum entanglement is enhanced, the flux and energy transfer efficiency are increased. Finally it is surprising that the non-equilibriumness quantified by voltage has a non-trivial contribution on strengthening the entanglement, which is attributed to the non-local feature of the quantum curl flux.


Introduction
The non-equilibrium system is everywhere around us, ranging from plants, animals to our Earth, Sun and galaxies. The non-equilibrium quantum processes are important in physics, chemistry, biology, even sociology and economics [1][2][3][4][5]. On the large and small scales, there are abundant examples for non-equilibriumness in acton: baryon genesis in early universe [6][7][8], transport and phase transition of quark-gluon-plasma [9][10][11][12][13][14], the black hole evaporation [15,16], transport in stars etc. On the mesoscopic level at the nano-scale, the electron and thermal transport in quantum dots and single molecule reactions have been recently investigated in both theory [17][18][19] and experiments [20][21][22][23][24][25]. The energy transfer in photosynthesis from light harvesting complex or Fenna-Matthews-Olson (FMO) complex to reaction center is with remarkably high efficiency, sometimes reaching >90% [26][27][28], which is still a mystery despite numerous studies on this subject [29][30][31][32]. However, the poor understanding of the global nature of the non-equilibrium quantum systems posts a challenge on the investigation of the underlying mechanisms of these problems. equilibriumness induced by environments leads to the enhancement of entanglement. We attribute this environment-assistant improvement to the non-local feature of the quantum flux.

Stochastic dynamics of non-equilibrium quantum systems from QMEs
For the purpose of illustration of our non-equilibrium quantum population and flux landscape theory, we use a generic energy transfer model in dimer systems [63], as modelled by multi-molecule system here. In this section our aim is to derive the QME for this general system. To describe the energy transfer between donor-acceptor molecules in dimers, we use the representation of localized excitons, which can be used to describe the excitations of a single molecule from ground state. When the hopping between different molecules is considered, the motion of the excitons becomes delocalized. This leads to a new quasi-particle: Frenkel exciton. It can be interpreted as the superposition of the localized excitons. To realize the energy transfer, the quantum system is now coupled to two environments, consisting of harmonic oscillators with the linear dispersion relation, namely, ω = vk k where v is the speed of light (sound wave) in vacuum for photon (phonon) and k is the wave vector. One bath serves as the high-temperature heat source pumping energy into systems while the other one serves as the low-temperature heat source absorbing energy from systems. The scheme of the system is shown in figure 1, where the molecules in donor and acceptor are connected to high-and low-temperature environments, respectively.
The Hamiltonian of the system and environments reads as  where the molecules in donor and acceptor are labelled from 1 to N and + N 1 to M, respectively. η and βʼs are the annihilation and creation operators for excitons for the donor and acceptor respectively which obey the standard fermionic anti-commutation relations. The annihilation and creation operators a(b) and a b ( ) † † for the radiation or phonon environmental baths satisfy the standard bosonic commutation relations. Two distinct environments are coupled to donor and acceptor with different temperatures, the difference of which defines the non-equilibriumness quantified as non-equilibrium voltage. p and s denote the polarization of photons or phonons. εʼs are the on-site energy of excitations and Δʼs are the electronic coupling between different sites, which are induced by dipole-dipole interactions between molecules. Here we only include the diagonal disorder rather than the off-diagonal disorder, for approximation, in that the fluctuation of energies of atom or molecules are always dominant [64]. γʼs stand for the electronic coupling between the donor and acceptor molecules. In Figure 1. Schematic of energy transport in dimers, where the donor and acceptor connect to two radiation environments with temperatures T 1 and T 2 , respectively. Each sphere represents a site in the dimer. equation (2) the rotating wave approximation [65,66] was applied since we are interested in the long-time dynamics of quantum systems so that the energy broadening is small and then the virtual-quanta process can be ignored.
In order to perform the further derivation for the QME we first diagonalize the system Hamiltonian by Bogoliubov transformation [67].
where the unitary matrix introduced in equation (3) is to diagonalize the single-exciton Hamiltonian In terms of transformed creation and annihilation operators, c and c † , we introduce a new operator to write the interaction in compact form For simplicity, we focus on the first part of the interaction in equation (2), namely, the coupling to hightemperature bath, where the environments are assumed to be at thermal equilibrium. The procedure for the remaining part of equation (2) is formally the same as the first part. From equation (8) the correlation function of the environment can be obtained ⎦ ⎥ is the Bose average occupation on frequency ω at temperature T. By substituting the equation (7) and (9)) into QME, we obtain the density matrix evolution of the system after summing over the environmental influences (denoted by R 1 , 1 e Tr d˜( ),˜( ),˜( ) (0) e (10) we see that a dissipation term emerges up to the second order of the coupling strength after summing (trace) over the environmental influences. In local exciton basis, this reads where the relation in equation (3) was used and the damping rate is of the form: for ω > ν 0. The Weisskopf-Wigner approximation that the upper limit of the integral over time can be extended to ∞ was used since the integrand oscillates very rapidly for ≪ s 0 while the most dominant contribution originates from times ≃ s 0. In addition, the rotating wave approximation was also applied in deriving equation (11). Now let us include the second environment with low temperature T 2 , the procedure for dealing with this bath is the same as that for the bath with temperature T 1 , except the dependence of coupling constant on the site index, as shown in equation (2). Then in the interaction picture the coupling term takes the form of 2 2 for ω >  and therefore the QME for the system after tracing out the environmental influences is finally reduced to In what follows the reorganization energy that measures the interaction between system and environments is set to 6.37 cm 1 in the regime of weak coupling, in which the perturbation expansion is applicable. Nevertheless the perturbation approach will fail for the large coupling.

Landscape and curl flux decomposition of the driving force for non-equilibrium quantum dynamics
To solve the QME, it is convenient to project the equation of motion into Liouville space, in which the matrix is represented as column (vector) and the inner product is defined as 〈 | 〉 = A B A B Tr( ) † . Thus the QME obtained in the last section can be projected as: ρ ρ ∂ | 〉 = | 〉  t . In our system, the many-body state of excitons reads | ⋯ 〉 n n n , , , M But the features of our system make it possible to reduce the dimension of Liouville space. This is because that the coherence between many-body states in Fock space with different total occupation n are decoupled from those terms with the same total occupation n, in the QME in Liouville space. In other words, the evolution of coherence with different n are independent of those (including both population and coherence) with the same n. Hence the matrix  in our case is of the block form, which consequently reduces the dimension of Liouville , which is much smaller than 4 M for large M. Physically, the probability of multi-particle excitation is much lower than the single-quanta excitation, roughly estimated by Fermi's golden rule [66]. Thus the space of our system can be truncated up to two-quanta excitations, which leads to the further dimension reduction to be + 2 . This can be obtained by setting the upper limit in the summation in L dim( ) to be 2. After some calculations the expressions for the elements of  can be obtained from equation (14).
To study the non-equilibrium quantum dynamics in the follows we need to further explore the QME in Liouville space, by separating the population and coherence terms as where the forms of elements in the block matrices ʼs can be obtained from equation (14), (15) and here we skip the details to avoid redundancy. ρ p and ρ c represent the population and coherence components of density matrix, respectively. Then we can project the QME into population space, by eliminating the coherence components in equation ((16) from Laplace transform, and one obtains which indicates that the quantum dynamics leads to a memory effect that is independent of the random collision. For the long time behavior at steady state, however, the reduced QME in population space has a simple form, by evaluating the integral in time domain where the assumption of negativity of the eigenvalues of matrix  c is essential, in order to ensure the convergence of the integral in time domain [63]. For convenience we define ≡ − −      p p c c cp 1 and the matrix  serves as the steady-state effective Liouvillian in population space. The second part of matrix  originates from the coupling between population and coherence terms. Because the reduced QME in population space at steady state is of the same mathematical form as the classical master equation (CME), we can properly define the transfer matrix: . Obviously, the symmetric part obeys the detailed balance condition, which indicates the equilibrium, while the anti-symmetric generates the non-zero steady state flux, giving detailed-balancebreaking. Hence in open quantum systems the populations here quantify the probability of the states of the system and therefore the population landscape, which mainly governs the equilibrium part of dynamics. The quantum flux has rotational curl property due to the nature of the anti-symmetric part of the transfer matrix and can be further decomposed into a sum of flux loops [63].
The total transition rate matrix that describes the net non-equilibrium quantum flux between different pairs of states, can be defined as follow which is closely correlated to the anti-symmetric matrix T mn a , as defined above. c mn in equation (19) actually describes the net probability current from the state | 〉 m to | 〉 n , of which the form in equation (19) ensures the non-negativity of the net probability current.
is the ith closed curl matrix (circular) with each describing a closed loop flux, if the following conditions are satisfied: This flux-decomposition theorem was mathematically proven at the classical level [69] and applied to synchronized dynamics [70]. Because of the same mathematical formalism as the CME, it is still applicable for the quantum case, with the replacement of the classical rate to matrix  where the quantum effect is contained in the 2nd term in , as reflected in equation (18). Therefore the transition rate matrix introduced in equation (19) at steady state is the superposition of closed loops. As we will show later, the distribution of these loops plays an important role on characterizing the intrinsic non-equilibriumness and the macroscopic features, such as ETE, heat current and entropy production. As illustrated in the theorem above, the condition (2) actually means the unidirectional feature of the flux that breaks the time-irreversibility between each pair of states (i.e., = c c 0 12 21 reflects that only one direction of the flow can survive !). The condition (3) reflects the constraint on the summed flux being zero (not individual flux) at steady state and also the probability conservation, since the flux c mn describes how much probability of excitons propagates from state | 〉 m to state | 〉 n , based on equation (19).

Population landscape and its relationship to non-equilibriumness and coherence
As mentioned above, the state space of the system is of M-dimension, in which each quantum many-body state is labelled by a set of M numbers, namely, { ⋯ n n n , , , M 1 2 }. In order to visualize the population landscape, the Principal Component Analysis (PCA) [33] can be applied to project the data points into two-dimensional space, spanned by the first two principal components. We can do so due to inhomogeneous distribution of the data points from the truncations up to two-exciton population in the state space of our system. The population landscape of our molecular donor-acceptor system is illustrated in figure 2, by the definition of = − U P ln ss ss , where P ss is the steady-state probability distribution (the diagonal elements of steady-state density matrix), as in the similar form as the classical non-equilibrium dynamics [34,35]. Notice that the color in figure 2 from light blue to red indicates that the probability becomes smaller and the corresponding population potential landscape becomes higher. The area in dark blue is the background with very low probability. Additionally, the X and Y coordinates in figure 2 represent the first two principal components generated by the inhomogeneous data points based on PCA method and they are independent of the state populations, since the PCA method used here only deals with the vectors in Fock space, which are uncorrelated to the dynamics.
It should be pointed out that here the population does depend on the coherence which is a unique quantum nature as compared to classical zero coherence case based on the QME in Liouville space, even though it has the same expression or form as in classical non-equilibrium systems. The main difference between classical and quantum case lies in the coefficients. The distribution of populations reflected by the underlying landscape is quantified in figure 2. In general, the landscape provides a global picture to the system dynamics, mainly the equilibrium part of the dynamics. For example, two basins in the topography of landscape indicates the bistability of the systems while the topography of Mexican-hat shape indicates a limit cycle oscillation, in the classical non-equilibrium dynamics [34]. The landscape in figure 2(a) shows that in the coherent regime (electronic coupling ≫ exciton-photon coupling) there are several major basins (local minimum), in green with lighter color surrounded, which consists of the ground state and the single-exciton states where all excitons stay in one ring (donor or acceptor). The major state transitions can be seen as transitions from one basin to another. However the feature of several major basins (local minimum) becomes obscured in the incoherent regime (electronic coupling ≪ exciton-photon coupling) , as shown in figure 2(b), where only one major basin seems to appear. In this case, major state transitions are not expected. This reflects globally that the coupling between donor and acceptor is required if the efficient energy transfer is expected to be achieved. We will study this issue again in the following sections where we discuss the quantum flux and transport in details.
It is worthwhile to point out that the landscape is not only for visualization purpose. In fact, the number of the basins and barrier heights are global and physical quantities characterizing landscape topography, which can be quantified to explore their relationship with respect to the state of the system, controlled by, i.e., electronic coupling, non-equilibriumness and coherence. We will touch this issue in the following discussion.
By comparison between the topographies of the population landscapes for coherent transport, under different non-equilibrium conditions as shown in figures 2(a)-(c), we can see that the barriers between the basins are suppressed as the transition from near-to-equilibrium to far-from-equilibrium. This enables the farfrom-equilibrium to realize the energy transfer with higher efficiency than near-to-equilibrium case. In contrast, for incoherent transport under different non-equilibrium conditions, the major basins becomes obscured on the topography of landscape from near-to-equilibrium to far-from-equilibrium cases and single basin becomes dominant, shown in figures 2(b)-(d). Thus the energy prefers to flow into the global minimum of the landscape instead of flowing between different states at local minimum, which results in the low-efficiency energy transfer. Therefore it is further confirmed that the coherent transport is essential to optimize the energy transfer between donor and acceptor. In the moderate non-equilibrium conditions, there are more states but more difficult to reach.
On the other hand, in figure 2 the distribution of quantum coherence is shown to become more spread with larger mean values and fluctuations at strong electronic coupling regime between donor and acceptor, compared to that at weak strength of coupling. Since the probability landscape is coupled to coherence, therefore, the enhancement of quantum coherence correlated to the optimization the energy transport shown later is based on the topography of landscape mentioned above as shown in figures 2(a)-(b).

Flux landscape 3.2.1. Statistics of curl fluxes
As shown earlier, the detailed balance breaking and time-irreversibility are reflected by the anti-symmetric part of transfer matrix defined above, giving the non-equilibrium quantum flux. In order to study the flux further, we adopt the definition of transition rate matrix in equation (19). At steady state, it is straightforward to verify that the transition rate matrix c satisfies the conditions in theorem 1, hence the non-equilibrium flux can be decomposed into superposition of the components of loops. Physically the loop feature here is due to the zero total flux with no sink or source condition at steady state. Here, we will investigate the flux landscape, reflected by the distribution of the magnitude of loops. On the other hand, the quantum entanglement will be briefly mentioned in the following discussion, but a detailed discussions will be given in next section rather than here, in order to avoid the redundancy.
From figures 3(a), (c) and 11(a) we can see that the non-equilibrium voltage from the difference in the temperatures between the two environments plays a significant role on the determination of statistics of the curl flux, namely, the mean values, the fluctuations and the whole distribution. To quantify the statistical fluctuations of the curl flux we define the relative fluctuations with respect to the mean of curl flux where  DA is the component of loops which passage the states in both the donor and acceptor rings.
2 represents the second order fluctuations of curl fluxes. Furthermore, both of the mean values and the fluctuations of quantum curl flux increase when the effective non-equilibrium voltage increases. This shows that the degree of non-equilibriumness drives the quantum transport. In other words, the distribution of curl flux that passes through both donor and acceptor becomes broader as the system deviates from equilibrium. This is shown in particular in figures 3(e) and (f) where we see the distribution of curl flux passing through both donor and acceptor in three regimes: near-to-equilibrium, intermediate and far-from-equilibrium. Here we focus on the curl fluxes which pass through both the donor and acceptor since they are physically relevant to the energy transport. Figure 3(c) shows that (i) the relative fluctuation decreases as the transition from the regime of incoherent to coherent transfer at far-from-equilibrium conditions. This indicates that the noise, quantified by relative fluctuation, is effectively suppressed at larger electronic coupling (coherent regime). (ii) The optimal energy transfer occurs in the regime of coherent transport, owing to the monotonic decreasing of noise with respect to voltage. This also reveals the non-trivial quantum effect on the energy transfer process. In contrast, the randomlike behavior of relative fluctuations with respect to the voltage in the regime of incoherent transport leads to the large noise and the resulting low-efficient energy transport. Moreover, as seen from figure 3(d), the entanglement is improved as the noise is suppressed at large electronic coupling, which means that the quantum entanglement is strongly correlated to the optimization of noise from incoherent to coherent transfer.
In figures 3(b) and (d) we illustrate the correlation between curl flux and quantum entanglement. As we can see, the intrinsic non-equilibriumness and quantum transport, quantified by the curl flux can be promoted by the entanglement, by controlling the voltage. On the other hand, it is also shown that environment-assistance is essential for quantum effects on improving the steady-state quantum flux and transport. This is particularly interesting even though the environment led to the dynamical dephasing as already pointed out before, while enhancement of quantum transport here is at the steady state level. show that the non-equilibriumness quantified by the voltage drives and increases the oscillation frequency of limit cycle, in both coherent and incoherent regimes. This also reflects that the larger energy pump provides the maintenance of limit-cycle oscillation mode with higher frequency. On the other hand, the coexistence of the limit-cycle oscillation and the energy transfer is also shown since the vanishing coupling between donor and acceptor leads to the absence of energy transport, as illustrated in following section. Without the flux, the limit cycle oscillation can not sustain (flux is correlated with the coupling as shown earlier).
As a limit-cycle mode, it implies an underlying funnelled flux landscape, as shown in figures 4(c) and (d). Here we use the convention of negative flux for showing the flux spectrum in analogy to the energy spectrum (i.e., hydrogen atom). To quantify the degree of funnelness of the flux landscape, we define the robustness ratio of funnel as the ratio of gap measuring the difference between the dominant flux and mean of the other flux loops, versus the fluctuation of flux distribution in the flux landscape From the picture of the landscape in general, it is proper to define the spectrum or statistical distribution of the magnitudes of curl fluxes as flux-landscape, in order to illustrate the global behavior of multiple curl fluxes and be consistent with the language used before. Our results show a single loop flux dominating the other ones, which can lead to a sharp peak observed in the power spectrum experiments. This signature further demonstrates the existence the limit cycle oscillation.
The physical insights obtained through this approach give a deeper understanding of the non-equilibrium nature of the quantum systems. It goes beyond the direct quantification of the detailed balance condition by the flux at the mean level. With the consideration of higher order statistics and the flux-landscape, the statistical nature of the non-equilibriumness beyond the mean description is provided. Furthermore, a limit cycle mode is a typical example of the non-equilibrium systems. Both landscape and flux drive the limit cycle at the mean level. The flux-landscape gives a further statistical description of the limit cycle oscillation as the single dominating flux loop among the existing others, which can be probed from the experiments.

Kinetic paths
Here we discuss the kinetic pathways in the similar spirit as the work given in [49,50]. Figure 5 presents the dominant kinetic paths for the energy transfer between donor and acceptor, which are quantitatively obtained using the maximum flux between a pair of states to determine the next state, after fixing the starting state each time. This is in the similar manner as given in [50]. In the energy transfer processes such as that in light harvesting complex in photosynthesis, there could be several molecules simultaneously performing the function. Therefore, the corresponding kinetic paths may start with different starting positions. The similar case for single-excitation manifold was studied in [50]. Here, we only show some of the dominant paths in the figure which are enough to reflect the multi-path feature and interference effect. In particular, for the set of the parameters γ γ , with the starting state being site 1 as the example shown in figure 5( (2) 11 s, respectively. Notice that the time spent on the paths above is estimated by the inverse of the flux on the path. Therefore one can conclude that as the barriers between basins on the topography of landscape become lower from incoherent transport to coherent transport (shown in figures 5(b)-(d)), the passage time on the dominant paths drops down which on the other hand, leads to higher efficiency on energy transfer. In other words, the smooth landscape leads to the shorter path time while the rough landscape leads to longer path time. On the other hand, by comparison to the statistics of curl fluxes shown in figures 3(a) and (c), the mean value decreases and noise of curl fluxes increases when the system transits from coherent transport to incoherent transport. Hence it illustrates that the sharp-distributed curl fluxes with large mean value and small noise leads to the shorter path time, namely, higher efficiency on energy transfer.
For each set of parameter γʼs, as we can see clearly in figure 5,   . This is consistent with 2D electronic spectroscopy, which confirmed the coherent multi-path mechanism in energy transfer in FMO [68]. The multi-path structure helps to keep the efficient energy transport, even if a particular path is destroyed, since the pathways have different starting positions so that the function of energy transfer can still hold with the destruction of one pathway.
From the QME derived before one can see that the electronic coupling between excitons provides a measure on the coupling between population and coherence [54], which alternatively leads to the term −    pc c cp 1 . In particular, the magnitude of the coupling coefficient between population and coherence increases with the improvement of electronic couplings. Coherence is quantified by the presence of the non-diagonal elements of the density matrix, giving the quantum contribution distinctly from the pure classical case where only diagonal elements of density matrix as populations are needed. So the presence of the off-diagonal in addition to the diagonal elements of the density matrix is purely originated from the quantum nature. The interference is analogous to coherence in the sense that they both represent the pure quantum nature. Conventionally the origin of the interference comes from the additive of the amplitudes but not at the intensity level, leading to the intensity interferences from one path to another. The total intensity in our system is not only made of the sum of the individual intensities, but also the interferences among them. It is expected that the absence of coherence will lead to the death of interference. When the coupling between the population and coherence increases, the interference between the paths should be strengthened. We can see this through the sum of the fluxes ∑ ≠ c m n mn . The flux itself reflects the transition rate flux as a combination of populations and rates between states. We can illustrate and quantify the interference from the following two arguments: first, the individual flux can be decomposed of two contributions in matrix  defined before. Second, for each path, we can sum over the flux on the path. For multiple paths, we can sum over all the fluxes on the different paths. The sum of the fluxes provides a measure on the total flux flowing into the specified state | 〉 n from different pathways. This is in the same spirit when discussing the interferences. Then the interferences among pathways occurring at the state | 〉 n are governed by the total coherence contribution in the sum of the fluxes. The larger the coherence is, the larger the interferences are. For example, in the mesoscopic electron transport the interference between electron waves occurs among the nodes of the circuits, where the interference is contained in the current conservation reflected by the sum of currents in all circuits. In a limiting case that the electronic coupling vanishes, it is straightforward to show that in QME the population and coherence terms are decoupled. In other words, this indicates the independent dynamics between population and coherence, so that coherence has no contribution to the non-equilibrium flux. The interference is not expected to be present as discussed earlier. To quantify the interference between the kinetic paths in the flux network, we can introduce the following constant similar to [50]  Before leaving this section, it should be pointed out that the quantum entanglement is strongly correlated to the flux landscape and curl fluxes statistics, as illustrated before. In the next section we will provide an elaborate discussion on entanglement, described by the so-called entanglement entropy.

Entanglement entropy
In this section we will study the quantum entanglement, which is quantified by entanglement entropy. This is not only a significant feature of quantum systems in general, but also is crucial for many-body correlation, due to its absence in single-particle case. In other words, one could say that in quantum many-body systems, quantum signatures are reflected by two features: coherence (interference) and entanglement. In general, by assuming A and B components of the whole system, the von Neumann entropy of the subsystem A is of the form ρ ρ = − S k Tr ln

A B A
A where ρ A is defined as density matrix of the subsystem A after tracing out the subsystem B, namely, ρ ρ = Tr , after partial trace, which gives As is shown in equation (27), for ψ AB e the system A is in mixture ensemble after partial trace while for ψ AB d it is still in pure ensemble. Thus the von Neumann entropy of these two states are = S k ln 2 e B and S d = 0, respectively. Therefore for ψ AB e we can say the system is entangled while for ψ AB d we say the system is unentangled, so that the entanglement entropy provides a quantitative measure of the degree of quantum entanglement. In our case, the density matrix of excitons at steady state reads  equation (19), in order to study the behavior of quantum flux. Therefore the flux captured by the molecules in acceptor ring can be written as  (19). Figures 8(a) and (b) show the effect of voltage, namely, the temperature difference between two environments, and the quantum entanglement quantified by entanglement entropy on the non-equilibrium quantum flux flowing into the acceptor, respectively. As we can see, the flux is significantly promoted in the regime of far-from-equilibrium where the voltage becomes large. This is due to the fact that the more effective energy pump is generated as the system becomes far from equilibrium.
On the other hand, the multi-excitation and delocalization of excitons in our system would drive the quantum states to be entangled and this is quantitatively described by the entanglement entropy. In distinction from the classical physics, the quantum entanglement has no classical analogy and it undoubtedly affects the behavior of non-equilibrium quantum flux. The quantum entanglement leads to a great improvement of flux, as shown in figure 8(b). This is directly owing to the enhancement of entanglement by environments. Here the conclusion that we could first draw is that the environments do play an important role on the enhancement of quantum entanglement, which then is crucial for the promotion of quantum flux.
Next we will study the ETE which is governed by the flux captured by acceptor discussed above.

Transfer efficiency
From the definition and discussion of curl flux we know that it provides a measure on how much energy (excitons) is transferred from donor to acceptor. Therefore the ETE can be introduced in terms of flux captured by the acceptor, so that represents the decay rate back to ground state from the excitations in donor. Based on equation (35) we know that there is no energy transport at equilibrium.
In order to reveal the behavior of ETE in terms of non-equilibriumness as well as entanglement, we performed the numerical simulation based on equation (28), which results in the figure 9. In particular, figure 9(a) shows that high quality of energy transfer with the efficiency of 70% is always realized in the regime of far-from-equilibrium, instead of near-to-equilibrium case studied before by the application of non-equilibrium Green's function method. Figure 9(b) illustrates the strong correlation between quantum entanglement and the enhancement of ETE, through the bridge of environments.
As we can see in figure 9, on the other hand, the transition from incoherent to coherent transfer leads to the increase of ETE as well, which we attribute to the improvement of mean value and fluctuation of nonequilibrium curl flux by the electronic coupling, as discussed on the smoother landscape and the larger spreading flux before. Physically this behavior can be understood by a limiting case: electronic coupling = 0, which indicates that the donor and acceptor are decoupled, so that no energy transfer will occur.

Macroscopic energy dissipation
On the macroscopic level, the possibility of the experimental realization makes it essential to measure the energy dissipation (heat current) and entropy production rate (EPR) reflecting the degree of non-equilibriumness and transport features. Therefore we need to explore the heat current, EPR and associated quantum thermodynamics. From the 1st and 2nd laws in thermodynamics where the former describes the energy conservation and the latter indicates the positivity of total entropy production in the non-equilibrium process, one has where based on QME in equation (14) and (15) ⎦ are the energy flowing into system per unit time from the high-temperature and low-temperature environments, respectively. The Q 1 and Q 2 can be exactly evaluated using equation (14). S˙and S˙t are the rate of system entropy and total EPR, respectively. Due to the assumption of large environments which lead to negligible back influence from system to environments, the entropy flux flowing from system to environments is σ = − + . Here Q 1 is called the energy current. The behaviors of heat current and EPR at steady state with respect to non-equilibrium voltage are shown in the small and large figures in figure 10(a), respectively. It is found the asymptotic behavior of 2 which coincides with the Fourier's law. Moreover we can see that the large electronic coupling in the regime of coherent transport is essential to enhance the heat current and EPR, besides non-equilibriumness quantified by voltage. This is owing to the similar behavior of mean of curl fluxes discussed in section III B(1). On the other hand, as is known that the two-exciton process was included in our model, the quantum entanglement must have strong correlation with the macroscopic heat current and EPR here, which is confirmed in figure 10(b). Both the heat current and EPR are increased when the entanglement entropy increases by controlling the voltage.

Summary and remarks
In this work, we developed the population and flux landscape theory for non-equilibrium quantum systems at steady state. We illustrate our theory with an example of donor-acceptor energy transfer. We developed a symmetric-anti-symmetric decomposition of the driving force in non-equilibrium quantum dynamics. We showed that the symmetric part of the driving force generated the population landscape which mainly governs the equilibrium part of dynamics, while the anti-symmetric part generated the steady-state quantum flux, leading to the detailed-balance-breaking and time-irreversibility. We found that the large electronic coupling optimizes the topography of population landscape for energy transport, by controlling the voltage quantified by − T T 1 2 in experiments. The non-equilibrium quantum flux generated by the anti-symmetric part has the curl nature due to its close loop nature. We study the trend of changes in population and flux landscape with respect to the nonequilibrium voltage and electronic coupling. Enhancing the non-equilibrium voltage and electronic coupling in general facilitates the quantum transport by reducing the population landscape barriers between major states and increasing the average flux.
The curl flux can be further decomposed to the sum of closed loops (curl fluxes). Due to the emergence of multi-loops, in contrast to single loop [63], we investigated the flux landscape, generated by the distribution of curl fluxes. In our case both the coherence and quantum entanglement are strongly correlated to the curl fluxes as well as the quantum transport, due to the multi-excitation. In particular, the distribution of curl fluxes are broadened. This is reflected by the fluctuations and the mean values. They increase when the coherence or entanglement quantified by entanglement entropy is increased by the increase of the non-equilibrium voltage. Moreover, the significant suppression of relative noise of curl fluxes which optimizes the energy transport occurs at the large electronic coupling and far-from-equilibrium regime, where the coherence and entanglement are enhanced. To our surprise, the non-equilibriumness induced by the environments leads to the enhancement of entanglement at steady state, despite of the dynamical decoherence. We attribute this novelty to the non-local feature of quantum flux. On the other hand, we found the dominant curl flux cycling inside donor with a gap reflected in the funnelled flux landscape implies a limit-cycle oscillation. The oscillation period is strongly correlated to non-equilibrium voltage and electronic coupling. Faster oscillations emerge with large electronic couplings and large non-equilibrium voltage. We quantified kinetic paths. The multi-path structure was uncovered. We also quantitatively showed that the interference between pathways due to the term −    pc c cp 1 originated from the quantum coherence effect. Moreover, we found that the transition rate for single-exciton transfer increases (after a local minimum) as the voltage increases while that for double-exciton transfer the rate shows a decreasing trend as the voltage increases.
Finally, we study the quantum transport such as ETE and dissipation quantified by the EPR. Both ETE and EPR have a significant increase when the coherence and entanglement are enhanced by the non-equilibrium voltage. This is due to the strong correlation of ETE and EPR to curl fluxes.