Efficient quantum transport in a multi-site system combining classical noise and quantum baths

We study the population dynamics and quantum transport efficiency of a multi-site dissipative system driven by a random telegraph noise (RTN) by using a variational polaron master equation for both linear chain and ring configurations. By using two different environment descriptions—RTN only and a thermal bath+RTN—we show that the presence of the classical noise has a non-trivial role on quantum transport. We observe that there exist large areas of parameter space where the combined bath+RTN influence is clearly beneficial for populating the target state of the transport, and for average trapping time and transport efficiency when accounting for the presence of the reaction center via the use of the sink. This result holds for both of the considered intra-site coupling configurations including a chain and ring. In general, our formalism and achieved results provide a platform for engineering and characterizing efficient quantum transport in multi-site systems both for realistic environments and engineered systems.


Introduction
Recent research has shown that quantum transport efficiency can be enhanced by the help of environmental noise, an effect known as environment-assisted quantum transport (ENAQT) [1] or dephasing-assisted transport [2,3]. The phenomenon has received particular attention in light-harvesting systems, such as photosynthetic complexes, but is also relevant in engineered systems, e.g. superconducting circuits and trappedion systems.
For instance Nesterov et al [7], have studied the dependence of efficient energy transfer (EET) on the correlation properties of the random fluctuations of the protein environment by modeling those fluctuations by RTN for a donor-acceptor system (i.e. a two-level system). They found that in case of strong-electronic coupling, the independent noise fluctuations on the site energies may be more effective in helping EET than collective noise.
The effects of RTN on two-level systems are also studied in the field of quantum information [26][27][28], with particular attention to solid-state devices [29]. Uchiyama et al [17] have analyzed the effect of spatial and temporal correlations on EET in a multi-site model by using a Ornstein-Uhlenbeck noise process to describe the environment, and observe that negative spatial correlation of the noise is the most effective in helping the EET. The effect of RTN on transport via continuous-time quantum walks has been studied on lattices [30,31], also in presence of spatial correlations [32]. Its effect on the non-Markovianity of the dynamics of the spin-boson model has also been considered [33].
Besides theoretical studies, ENAQT has been reported in experiments with superconducting quantum circuits [34], photonic setups [35], classical oscillators [36] and trapped-ions [37]. All these works show that in certain parameter regimes the efficiency of the energy transfer is enhanced by environmental noise.
The assumption that the environment can be modeled as a collection of independent harmonic oscillators might not be adequate for various situations where there exist different types of motions that can not be handled by a harmonic approximation. For example, in light harvesting systems such as the Fenna-Matthews-Olson (FMO) complexes, the transfer of exciton might be affected by large amplitude modes, e.g. motion of the protein scaffold. Besides, the concept of the so-called real environment as both a quantum mechanical thermal bath and a classical noise source representing anharmonic effects for the multi-site model in the open quantum system is lacking [38], and the important question remains: can one observe ENAQT-like phenomena in an environment treatment consisting of both a stochastic noise and quantum thermal bath?
To answer this question, we study numerically an exact noisy multi-site spin-boson model in the variational polaron framework [39] deriving the corresponding master equation, and also add RTN on the site energies. We then study the transport efficiency of a single electronic excitation, and compare the effect of classical and quantum noise, and their interplay. Our results show the presence of a non-trivial interplay between the quantum thermal bath and the classical noise-displaying enhanced transport efficiency. The used formalism and results go beyond a recent study [17] which demonstrated the improved transport efficiency when using (anti)correlated Ornstein-Uhlenbeck noise only. Our study does not necessarily model a specific real system, but considers ingredients that are in general plausible in the different scenarios described above, in order to shed light on the ENAQT phenomenon and to provide a platform for engineering and characterizing efficient quantum transport in multi-site systems both for realistic environments and engineered systems.
The article is organized in the following way. In section 2 we discuss the model and derive the variational master equation. Then in section 3 we present the results of the study, by first looking at the dynamics of the system followed by detailed results on the efficiency of the transport and exciton trapping time. Finally, we conclude with discussion and outlook in section 4.

Model
We consider a noisy multi-site spin-boson model where each site interacts with a separate thermal environment, and the site energy is modulated by RTN, as represented in figure 1(a). The total Hamiltonian of the system is given by where H S (t) is the (time-dependent) Hamiltonian for the spin system, H E is the Hamiltonian for the bosonic environment and H I is the interaction part:  , is the annihilation (creation) operator for the kth oscillator of the nth site whose state is described by | ñ n , V nm is the electronic coupling between nth and mth sites, and g n,k is the interaction strength between the exciton on the nth site and kth harmonic oscillator of its bath. We assume that the energy of each site fluctuates with RTN, i.e. ò n (t)=ò n0 +Ω n α n (t), where ò n0 is the bare site energy of nth site and Ω n is the noise amplitude at nth site. The RTN is a stochastic process that describes a bistable fluctuator that jumps between two values α=±1 with a certain switching rate ν. It is characterized by zero average, ( ) a á ñ = t 0 n , and an exponentially decaying auto-correlation function where the correlation time of the noise is τ c =1/ν.
The bosonic environment of the system is treated as a collection of non-interacting harmonic oscillators. The system-bath interaction is characterized by the spectral density, defined by ( ) In the continuum limit, the strength of the coupling of the system to the environment is measured by site specific reorganization energy ( ) . For the present study, we consider a spectral density J com (ω) consisting of two parts. The first part defines the broad-range background modes (overdamped) [40], while the second one describes a discrete underdamped vibrational [41] mode which interacts with an Ohmic environment with cut-off frequency Λ: is the exponential integral function while parameters S and X determine the peak magnitude of J bg (ω) and J vib (ω), respectively. Moreover, ξ and Λ act as damping factors for the discrete oscillator which corresponds to underdamped mode, whilst ζ governs not only the position of the underdamped mode, but also the magnitude of J vib (ω) in the present problem. Each site is assumed to have its own independent environment which is described by site specific spectral function parameters in equation (2). Figure 1(b) shows the spectral densities of the quantum baths.
When the ratio between the electronic coupling V nm and the bath reorganization energy E n r is small, the use of Redfield master equation is adequate [42,43], while the full polaron master equation [44] is useful when the ratio is large. In many photosynthetic systems, such as the FMO complex, the convenient coupling regime corresponds to the intermediate regime [45]. In order to describe the dynamics of the reduced density matrix of the current system, we adopt the variational polaron master equation for the multi-site spin-boson model derived by Pollock et al [39]. This is based on using the variational polaron transform , a n d , ,

*
In the transformed interaction HamiltonianH I , we have two types of interaction terms;H L is similar to the weak interaction term, whileH D is the full polaron term. B nm is the environmental displacement operator which includes the expectation values  n , defined as The electronic coupling V nm and site energy ò n (t) are renormalized respectively by  n and R n , the latter being defined as: , , and f n,k are left as free optimization parameters which are determined numerically by minimizing the contribution ofH I to the free energy [39].

The interaction HamiltonianH I in the variational polaron frame is assumed to be˜å
N is the number of sites, S i and E i are the system and environment operators, respectively. One can use the projection operator method to obtain a master equation for the system density matrix˜( ) with the Hamiltonian in equation (4) in the Schrödinger picture as [39]: [1,22,23]-an anti-Hermitian operator that is used to describe the excitation leaving the system permanently. Below, for the studies of population dynamics, this part is excluded, while it is included for transport efficiency and trapping time studies. In the latter case, the trapping is assumed to occur only on site n with rate κ. Above, are the bath correlation functions having four non-zero types (see equations (17) is the jth system operator in the interaction Hamiltonian at time t which is obtained by using the time evolution operator presence of RTN on the site energies, it is not possible to derive a simple expression for U(t) for the given system Hamiltonian.
We adopt the following procedure to calculate U(t) for each RTN noise realization. Let the state flipping times of the RTN be {0, t 1 , t 2 , K, t F } in the time interval [0, t F ]. Since the magnitude of α(t) is either 1 or −1 in between any two of those flipping times, H S (t) is constant in that sub-internal and U(t) can be easily calculated as ( ) -H t exp i S . The time evolution operator in the interval [0, t 1 ] would be either U + (t) or U − (t) depending on whether α(t) is in state +1 or −1. In general, if t n <t<t n+1 , the time evolution operator U(t) would be equal to if the noise starts from negative values and the number of noise flips in the interval (0, t n+1 ] is odd. The master equation (7) is derived for an arbitrary number of sites. In the following, however, we will focus on a three-site system, in different network configurations, to study the population dynamics, transport efficiency, and exciton trapping times.

The dynamics and steady-state behavior in three-site system
In this section, we study the dynamics of the system and its steady state in the absence of a sink (˜= H 0 trap ). We calculate the time evolution of the system numerically by using the ensemble averaging method based on the average over the noise realizations for a system in which site energy differences are comparable with the coherence electronic site couplings at low temperature (k B T=1 ps −1 ). The system-bath couplings are assumed to be site dependent, with reorganization energies that are also comparable to the electronic couplings. Although it is customarily assumed that the ratio between the reorganization energy (E r i ) and the electronic coupling (V i ) describes the strength of the system-bath coupling [46], this is not straightforward for the system where each site has a spectral density with different interaction constant. In the present framework, the system-bath interaction strength can be gauged by analyzing the tunneling renormalization factor B i , polaron shift R i , and the sitedependent variational transform functions. As the system-bath coupling increases, B i tends to zero, while R i approaches negative of the reorganization energy.
The system Hamiltonian and the spectral parameters of the thermal environment are given in table 1. As can be seen from the table, site 2 has the strongest bath coupling, while the bath couplings of site 1 and site 3 are intermediate based on values of B i and comparison between R i and E r i . We carried out a detailed study of uncorrelated and anti-correlated RTN affecting the three sites. Uncorrelated noise proves to be detrimental for the excitonic transport (see also [7,17]), hence in this paper we present the results obtained by a collective noise model. We assume that each site energy level is fluctuating due to the effect of collective RTN with the form {1, −1, 1}Ω: this means that the noise at site 2 is anti-correlated with those at sites 1 and 3, which are fully correlated. The results presented below were obtained by averaging the solutions of the master equation (7) over 100 realizations of the RTN process. The convergence of the averages as function of number of noise realizations has been checked carefully.
We consider two different environmental models-RTN only and bath+RTN-and the dynamics of the populations for the three-site system using both linear chain and ring configurations are given in figure 2. We use high noise frequency and intermediate noise amplitudes, and initially the exciton is localized at site 1. In figures 2(a) and (c) we display the time evolution of the populations for the noise only model, which approaches the maximally mixed state for both the linear chain and ring geometry networks in the long time period, even if those populations become qualitatively different in the short time period. This result is expected due to the nature of RTN noise: it is symmetric, hence it corresponds to the infinite temperature [9]. However, the effect of the RTN on the system in contact with the thermal environment has two components. Since the site energies are changed in a time dependent way due to the external noise, the system exists in a non-equilibrium state. External noise, depending on its amplitude and frequency, might change the energetic order of the states of the system. This can be seen in figures 2(b) and (d), where it is found that the exciton is delocalized over all three sites with comparable populations for the linear chain. For the ring geometry, the steady-state is quite different; here the exciton is shared by the sites 2 and 3. Comparing the steady-state populations of site-3 for the RTN only and bath +RTN models, for both interaction geometries, one might deduce that the excitation transport is enhanced with simultaneous action of the quantum bath and the external noise.
To further explore the noise dependence of the system dynamics, in figure 3 we display the steady-state populations as function of noise amplitude Ω and noise frequency ν for the RTN only and bath+RTN cases using again linear chain and ring geometries. In the RTN only case (figures 3(a) and (c)), the steady-state populations tend to approach the maximally mixed state for both the linear chain and ring geometries at large ν and Ω, which is similar to what we observed in figures 2(a) and (c). Besides, it is clearly observed that at the low frequency limit (where the correlation time of the noise is larger than the dynamics time), the transport is prevented in the linear chain geometry in figure 3(a). This result stems from the fact that the anti-correlated noise increases the energy difference between the sites with increasing Ω, which prevents the transfer of the exciton to the other sites above a certain noise amplitude. Table 1. Hamiltonian parameters, the calculated renormalization constants of the three-site system, equation (4), and the parameters of the site specific bath spectral functions, equation (2). V 1,2,3 correspond to V 12 ,V 23 and V 13 . Also, V 13 has two values (0 or v) depending on the type of coupling configurations. B i , σ i , ξ i ,S i and X i are dimensionless, while all the other parameters are in ps −1 . Note that w c i , , σ i , S i are the parameters of J bg (ω), while X i , ξ i , Λ i and ζ i describe the discrete mode. v is set to 10 for the calculations in the present study.
On the other hand, when having thermal bath in addition to RTN, the steady-state populations have very interesting behavior depending on the site coupling configuration, see figures 3(b) and (d). Here, the Ω=0 curve displays the steady-state populations for the system in contact with the thermal bosonic bath only-these  dynamics have been studied by different authors over the years [2,3,47]. The system Hamiltonian in the variational frame without the external RTN indicates that the site-2 has the lowest energy in the thermal equilibrium (see table 1) and also the energy of site-3 is close to that of site-2. This explains why site-2 has initially (Ω=0) the highest steady state population for the linear chain.
Moreover, we observe that for both the linear chain and ring configurations, the exciton is delocalized over site-2 and 3 with comparable probabilities when the external noise has absent or small amplitude, see figures 3(b) and (d). As the noise amplitude increases, figure 3(b) displays that P 1,ss , P 2,ss and P 3,ss eventually relax to the maximally mixed state for the nearest-neighbor coupling network, as is the case for the RTN only modeleven if site-3 has the maximum value at around Ω=10. However, when going to ring configuration with coherently coupling also site-1 and site-3, one can see that the steady-state population of site-3 reaches the maximum value after crossing the resonance boundary of the bare system for the ring geometry in figure 3(d).
Understanding the interplay of coupling geometry (e.g. the additional pathways for interference, as well as different structure of the energy levels introduced by the ring geometry), thermal bath and RTN is a non-trivial task that is worth investigating in future work. Thereby, we find that site-3 is clearly the most populated one when having high amplitudes of the noise for the ring geometry and the model includes both the RTN and the bath.
In general, this means that RTN+bath model seems to be very efficient in terms of reaching the target site-3. Does this mean that high occupation of site-3 also corresponds to the most efficient transport and shortest trapping time in our proposed model? To explore these questions and account for the excitation leaving the system, we add a sink to the target site-3. This is associated to having a reaction center where the energy is converted to chemical compounds. With all the above ingredients, we can then also see whether the transport can be enhanced even further when comparing to recent results of [17] and having also the action of the bath accounted for in addition of RTN.

Transport efficiency and average trapping time
We now investigate in detail the dependence of transport efficiency η and average trapping time á ñ t in the above described three site dissipative systems including also the influence of sink site which transfer the exciton to the reaction center. Here, 200 RTN realizations are used to average the reported quantities. The transport efficiency η essentially quantifies which fraction of the population has been transferred to the reaction center within a given interval of time. The average trapping time á ñ t , in turn, describes how fast is the transfer of the excitation to reaction center. Subsequently, these two quantities can defined and expressed in terms of the site populations in the following way [17,19,20]: Here, n=3 indicates the trap site and we have chosen tu=7 ps following similar arguments as in [17] . The trapping rate κ in equation (7) is chosen to have a relatively low value, 0.5 ps −1 , to be reasonably distant from the localization limit [48].
Notice that recently there has been a study properly accounting for the intensity of the incoming energy (sunlight) for realistic photosynthetic systems [12]. Here, however, we focus on engineered systems for excitation transport, therefore we employ the commonly used figures of merit defined above.
We first present the transport efficiency and average trapping time for the RTN only model, which neglects the thermal environment completely. The noise only model considered here is similar in spirit to the one studied by Chen and Silbey [5] and simulates the dynamics of the system H S (t) by ensemble averaging over different realizations of the RTN processes. Figure 4 displays the behavior of η and á ñ t as function of noise amplitude Ω at various noise frequencies ν for linear chain and ring configurations. Efficiency plots for both the linear chain ( figure 4(a)) and the ring (figure 4(c)) configurations show similar resonant enhancement of η with increasing noise frequency ν. The increase in transport efficiency with noise amplitude Ω is somewhat greater for the ring geometry compared to the one in the linear chain, as expected, probably due to higher number of decay channels in the former geometry. The enhancement is most pronounced when the noise amplitude Ω values are close to the site energy differences.
The Ω-dependence of η is noticeably different for both geometries at very low noise frequencies compared to that at intermediate and fast ones. At very low noise frequencies, the noise does not flip in the considered time interval and the averaging over the two solutions of equation (7) with α(t)=±1. For the linear chain in figure 4(a) the transport efficiency decreases to very low values with increasing noise amplitude indicating that high amplitude noise prevents the transport which is not observed for the ring geometry (see figure 4(c)). The numerical trends in average trapping times (figures 4(b) and (d)) are almost opposite of those observed for the efficiencies; at those Ω values where η is enhanced, á ñ t is strongly suppressed demonstrating fast transport of the excitation to the reaction center. This may be partially expected since efficient transport is coinciding here with fast transport-even though the two concepts are slightly different. In general, the results presented above show that a simple classical noise might enhance transport efficiency in a quantum setting similar to ENAQT phenomena.
Let us turn our attention now to the bath+RTN case. Figure 5 shows the behavior of transport efficiency and average trapping time with the same RTN parameters as in figure 4 but including now also the influence of the bath. Here, the Ω=0 line indicates only a thermal bath model which does not contain the effects of the RTN on the system-this framework has been also studied by some authors [2,3]. Without any external RTN noise, the transport efficiency has almost equal values between the two site configurations, see figures 4(a) and (c), and 5(a) and (c) at the values Ω=0. However, in the former case we have η≈0.50 and in the latter case η≈0.75. This difference in transport efficiency may be based on the presence of thermal fluctuations in the latter case which increases the probability of shifting the exciton to reach the trapping site.
In general, the set of results for RTN only, figure 4, and the bath+RTN model, figure 5, show that the dependence of both the efficiency and trapping time on the noise parameters is qualitatively different for different site configurations, and that the bath+RTN model has very rich ENAQT behavior. Typically, simulations of the ENAQT in FMO has found that there is an optimal region where the transport efficiency is enhanced by noise parameters depending on the proposed model. However, we show that the behavior of transport efficiency for RTN only and RTN+bath model is different from the results obtained in both a thermal bath [2,3] and the pure dephasing [23,24] cases. That is, in the RTN-only model, there is a single maximum in the dependence of the efficiency on the noise amplitude and noise frequency depending on the site configurations. The structure is richer for the bath+RTN model, where one can see resonance peaks.
It is also interesting to note that for both models the enhancement rate in the transport efficiency η in the linear chain configuration is relatively higher than in the ring geometry for the parameter range explored here. The enhancement rate is obtained by using the values at Ω=0 and the overall maximum value for any ν and Ω. For the linear chain, the rate is about 80% (for RTN only, 0.50-0.90 ) and 39% (for bath+RTN, 0.72-1.00), while for the ring geometry it is about 63% (RTN only, 0.54-0.88) and 30% (bath+RTN, 0.74-0.96). The enhancement is more prominent for the noise only model, which is expected because the presence of the quantum bath already helps the transport. On the other hand, it can be seen clearly from figures 4 and 5 that the noise frequency ν dependence of η is similar for both models; the transport efficiency increases with the increasing ν. This tendency is, however, not observed in the results corresponding to the steady-state populations of the system (see figure 3). This means that, although one could expect that high steady-state populations indicate efficient transport to the reaction center [23,49], this expectation is not always correct.
In order to show that the bath+RTN model significantly improves the efficiency and transport time compared to RTN only model, we display the differences Δη and Dá ñ t between the former and the latter in figure 6. The results show that the bath+RTN model is constantly dominant against the RTN only model for linear chain in terms of both measures, see figures 6(a) and (b) and the corresponding gray colored areas. For ring configuration in figures 6(c) and (d), the bath+RTN model is superior for the average trapping time (figure 6(d)) while for efficiency there are specific limited areas of noise amplitude and frequency, where RTN only dominates. In general, we conclude that when looking for efficient and fast transport, having the influence of both classical RTN noise and quantum bath provides the most useful avenue.

Conclusions
In this work, we have developed a formalism for and considered multi-site dissipative systems influenced by both the classical RTN noise and quantum baths. Technically, this requires the use of the variational polaron transformation for the derivation of the corresponding master equation, generating stochastic realizations of open system density matrix evolutions and using the ensemble averaging over noise realizations to obtain the final result for multi-site density matrix.
We used the method to study the efficiency of the energy transport in a three-site system by focussing and comparing the results from two different models-RTN only and bath+RTN model, where RTN contains full (anti)correlations between the three sites. One of the motivations was to go beyond a recent observation that correlated stochastic noise improves the transport efficiency [17]. The results for steady state populations, by using a chain and ring configurations, clearly demonstrate in both cases that there exists a considerable parameter region where the target site population is significantly increased when combining stochastic noise with a quantum bath. To take into account the transfer of an excitation from the target state to the reaction center, we consider similar models as above but now with a sink, and calculate the transport efficiency and average trapping time-which are commonly used to quantify the efficiency of the transport [17]. Here, the ENAQT behavior has quite a rich structure in terms of the frequency and amplitude of the random noise. The results clearly show that in almost all regions of the used parameter space, the combined influence of bath+RTN is superior to RTN only model when trying to achieve efficient and fast transport. However, many questions remain open for future studies. How does the efficiency change when increasing the number of sites and the complexity of the coupling configurations between the sites within our framework? What is the role of changing the temperature of the environment and the spectral character of the local baths of the sites? In general, our result opens the avenue to study, e.g., the questions above and when looking for ways and characterization of efficient energy transport in engineered and ambient systems.