Environmental engineering for quantum energy transport

Transport phenomena are ubiquitous throughout the science, engineering and technology disciplines as it concerns energy, mass, charge and information exchange between systems. In particular, energy transport in the nanoscale regime has attracted significant attention within the physical science community due to its potential to explain complex phenomena like the electronic energy transfer in molecular crystals or the Fenna-Matthews-Olson / light harvesting complexes in photosynthetic bacteria with long time coherences. Energy transport in these systems is highly affected by environmental noise but surprisingly not always in a detrimental way. It was recently found that situations exist where noise actually enhances the transport phenomena. Such noise can take many forms, but can be characterised in three basic behaviours: quantum, coloured or nonlocal. All have been shown potential to offer an energy transport enhancement. The focus of this work is on quantum transport caused by stochastic environment with spatio-temporal correlation. We consider a multi-site nearest neighbour interaction model with pure dephasing environmental noise with coloured and nonlocal character and show how an accelerated rate for the energy transfer results especially under anti-correlation. Negative spatial correlations provide another control parameter to help one establish the most efficient transfer of energy and may provide new insights into the working of exciton transport in photosynthetic complexes. Further the usage of spatio-temporal correlated noise may be a beneficial resource for efficient transport in large scale quantum networks.


Introduction
Energy transport is a fundamental primitive in our world operating on length scales ranging from the atomic scale to cosmological ones. It is at the core of natural life as well as our current technology. Any, even slight, improvements in transport efficiency can bring profound effects. In the natural world this can lead to a species dominating another, while in the technological world it can lead to lower energy consumption devices. Many of these improvements in the technological arena have been achieved by better device engineering reducing noise and imperfections. While this seems a perfectly logical approach, recently however it has been found, though nature already knew, that energy transport can be enhanced by adding environmental noise. This counter intuitive behaviour has shed new light on the conventional thinking that noises in transport phenomena should be removed, and it was indeed triggered by energy transport discussions of light harvesting complex of photosynthetic bacteria [1].
The intensive research on four wave mixing in the light harvesting complex, such as the Fenna-Matthews-Olson (FMO) complex in green sulfur bacteria [2][3][4], light harvesting complex in marine algae [5] or molecular crystals [6], has given us clues to understand the mechanisms for efficient quantum energy transport. The striking long lived coherence observed in these experiments strongly suggests the positive effect of environmental noise. Extending this concept to the dynamics of exciton motion with delta-correlated stochastic noise (white noise) proposed by Haken and Strobl [7], a simple theoretical model using pure dephasing was introduced to assist the excitation transport, referred as the environment-assisted quantum transport (ENAQT) [8] or dephasing-assisted transport [9]. To describe the longlived coherence observations, the model was further extended to include finite noise correlation times and lengths, including a coloured noise approach using dichotomic stochastic process [10], a reservoir modeling with an infinite number of quantum harmonic oscillators [11][12][13][14][15], and an analysis with finite correlation lengths [16][17][18][19]. Positive spatial correlations were shown to reduce the environmental effects, assisting the long-lived coherence. Both effects of temporal and spatial correlation together [20,21] were shown to be beneficial, and recently an artificial realization of the ENAQT model has been proposed [22,23]. The long-lived coherence implies the lasting energy exchange between sites, which however may affect the efficiency on the energy transfer rate. If the long-lived coherence can contribute to a higher transport efficiency, these contradicting evidences need to be somehow resolved. Spatial correlation might give us a solution to this [17-19, 24, 25]. Cao and Silbey [17] extended the interaction between the sites to include a phase relation and found the dependency of transfer efficiency on the phase, while the spatial correlation by propagation of environmental phonon modes [24] and the application of the extended ENAQT model with the spatial correlation to the photosynthetic bacteria [18,19] have been considered. Anti-correlation in noise has also been taken into account as the effects of anti-correlated (anti-phase) motion of two harmonic oscillators were investigated [25], and the anti-spatial correlations between the bacteriochlorophylls in the FMO complex could both positively and negatively influence the energy transport [19]. These results imply that the phase relation can affect the quantum energy transport, which leads us to the natural question of what is the best spatio-temporal correlation function.
In this manuscript, we consider a simple multi-site model for the environment-assisted quantum transport using spatio-temporal correlated stochastic noise processes and show how we can engineer the environment to enhance quantum energy transport. The fundamental difference with the conventional treatment is to exploit the negative parameter regime (anti-correlation) of the spatial correlation as well as including finite temporal correlations with exponential decay. For the two site model, using the 2nd order of time-convolutionless (TCL) master equation, we analytically find that the time evolution of the population under negative spatial correlation is significantly faster than uncorrelated and positively correlated spatial noise. Extending the spatial correlations to include anti-correlations provides another control parameter to find efficient energy transfer regimes. Using two independent measures for transfer efficiency, we find a significant improvement in transported quantity and elapsed time for negative spatial-correlation compared to noise with delta-correlation both in space and time as in the original ENAQT model. Extending to the three site model, we find that the negative spatial correlation between every other site shows the best transfer efficiency, while the negative correlation between the next nearest neighbors is worse.

Results
Model: Let us begin by describing a simple multi-site model for the environment-assisted quantum transport whose energy level diagram we schematically depict in Fig. (1). For N sites, the total Hamiltonian can be written as where is the cumulant expansion parameter of the time-convolutionless master equation with |n being the n-th excitation site, ω n the n-th site Larmor frequency, V nm the transition frequency between the n-th and m-th site, and f n (t) the fluctuating frequency on the n-th site which we consider as a stochastic process. We assume that the average of the frequency fluctuation for each site n is zero as f n (t) = 0, and the correlation function f n (0)f m (t) of these fluctuations can for convenience be described by a simple exponential decay as where to allow for both positive and negative (anti) spatial correlations we introduce the quantity c n,m with n, m = {1, 2 . . . N }. It is defined over the range by −1 ≤ c n,m ≤ 1 with the extremal values −1 (+1) corresponding to perfectly anti-correlated (correlated) noise respectively. ∆ n,m is the amplitude of the fluctuation, whereas τ c,{n,m} is the correlation time of the fluctuation. We also include damping with the rate κ at the end of the linear chain to trap the system in its ground state [17]. This corresponds to the excitation leaving the chain. Now averaging the density operator over the fluctuation using a time-convolutionless decomposition approach, we obtain the master equation for an arbitrary operator X withL 1 (t) = e iL0t L 1 (t)e −iL0t . The master equation given by (4) is inherently non-Markovian and incorporates both spatial and temporal correlations. In the limit that such correlations vanish, it reduces to a Markovian master equation(see Methods). Let us now consider a simple example. Two-site model : As the description of our multi-site model has been completed, we will now consider the simplest situation, N = 2. In the two-site model, when initially only the site 1 is fully excited, the time evolution of the probability of finding the second site excited is shown in Fig. (2). , we find that the anti-correlation between the sites makes transport to finish faster than the uncorrelated case. This is not unexpected when one examines the two time noise correlation function given by if the noise is the same for each site, while for the anti-correlated noise, the quantity becomes, φ(t) = 4 f 1 (0)f 1 (−τ ) , to be twice the size of the uncorrelated case.
To quantify this improvement in more details, let us examine two independent measures: (1)  and (2) ratio of transported quantity. The average trapping time is defined by [17] The average trapping time indicates how long the population of the initial excitation remains in the system. Since the excitation trapping only at the second site with rate κ leads to ∞ 0 dtρ 22 (t) = κ −1 , we focus on the average trapping time t minus an offset as κ −1 . This is drawn in Fig. (3 a) for varying degree of spatial correlation c 1,2 = c 2,1 = c. This indicates that the transport feature depends on the correlation between sites. Especially, as c approaches to −1, the average trapping time decreases, i.e. the transport finishes faster. Moreover, the quantum yield becomes larger as c approaches to ∼1, since the average trapping time is found to be inversely related to the quantum yield [17] given by q ≈ 1/(k d t + 1) where k d is the decay rate of recombination (see Methods).
Next, we introduce a quantity η to define an integrated probability that the excitation has been transported to the second site upto a time t u as Since the trapping rate κ is given by ∞ 0 ρ 22 (t )dt = κ −1 , η indicates the ratio of transported quantity between upto a finite time t u and completion of the transport. In Fig.(3 b), we chose t u to be 2000∆ when the transport to the second site is 98% completed for α = 1 and c = −1. We used it for all other evaluation in the figure. It is straightforward to observe from Figs. (2) and (3) the dependence of the energy transfer on both the degree of spatial correlation c and temporal correlation time α = ∆ · τ c . In both measure, there is a clear monotonic dependence on the degree of spatial correlation, with the best results occurring as we approach perfect anti-correlation (φ(t) is larger in this case). We also observe a non-monotonic dependence of the populations time evolution and the measures of transfer efficiency for temporal correlations. This indicates that there will be a condition for optimal energy transfer on the correlation time. For small and decreasing α, the energy transfer takes longer time due to the fact that fluctuation becomes too fast making it difficult for the energy gaps to be close and the transition probability between sites becomes smaller (this moves us towards the Markovian limit). For α larger than the optimal value, the correlation time becomes larger making the transition probability smaller again. Thus for the best performance one should operate near this optimal value. Now what happens when we have more sites?
Three-site model. The two-site model has shown how spatial correlations are potentially an important resource for efficient energy transfer. The natural question which follows to this would be whether this is true for when the system has more than two sites. We extend the model to three-site linear chain under the nearest neighbor interaction with the interaction strengths given by and two spatially-anti correlated noise cases. Fig.(4 a) corresponds to the uncorrelated case with c n,m = δ n,m for n, m = {1, 2, 3}, while Fig. (4 b) represents the case for the anti-correlated noise between the nearest neighbor sites with c 1,2 = c 2,1 = c 2,3 = c 3,2 = −1 and correlated between the end sites c 1,3 = c 3,1 = 1. Finally, Fig. (4 c) represents the case for the anti-correlated noise between the sites 1 and 3, and the sites 2 and 3, with c 1,3 = c 3,1 = c 2,3 = c 3,2 = −1 and correlated noise between the sites 1 and 2, c 1,2 = c 2,1 = 1. In each sub-figure, we find that the dependence of the dynamics on the correlation time for the two-site model remains for the three-site case: transport finishes faster as the correlation time becomes longer particularly in the short correlation time regime, α 1. Comparing Fig.(4 a, b, c), we find that the anti-correlated between the nearest neighbor sites typically shows the most efficient transport.
To explore these features systematically, we evaluate the average trapping time t . In Fig.(5), we show the contour plot of the dependence of t − κ −1 on the degree of correlation c n,m . Comparing Figs.(5 a,b), we find that "antiferromagnetic" correlation where the fluctuation is anti-correlated between the adjoint site shows the most efficient FIG. 5. Average trapping time t −κ −1 for three-site transport depending on the degree of spatial correlation for α = 0.3. a) shows this quantities dependency onc2,3(= c3,2) and c1,2(= c2,1) while setting c1,3 = c3,1 = c2,3c1,2. Similarly b) shows our quantities dependency on c1,3(= c3,1) and c1,2(= c2,1) while setting c2,3 = c3,2 = c1,3c1,2. Other parameters are the same as in Fig. (4). The shorter values of average trapping time occurs at the left hand bottom side corner in a) and right hand bottom side corner in b). This means the "anti-ferromagnetic" correlation case is the best choice for acceleration of transport in this case.
transport. The origin of the efficient transport lies at the higher transition probability between the adjoint sites with anti-correlation, which in the three site model corresponds to the"anti-ferromagnetic" correlation. There the average trapping time t − κ −1 for the uncorrelated case is ∼ 12 % longer than the "anti-ferromagnetic" case. This clearly shows the acceleration in the energy transport occurs with "anti-ferromagnetic" correlations, while the difference is smaller as α increases (5 % for α = 0.5). The question follows is how robust this efficiency in energy transport.
The most important feature is the stability against inhomogeneity in coupling between adjoint sites. To see this, we deviate V 23 /∆ from 0.1 to be 0.15 and 0.3, keeping V 12 /∆ = 0.1 in the three site model. The average trapping time decreases for any spatial correlation as increasing V 23 /∆, due to the stronger coupling strength accelerating the energy transfer. The difference of the average trapping time t between uncorrelated and "anti-ferromagnetic" case decreases from ∼ 12 % to ∼ 8.5 % for V 23 /∆ = 0.15 and ∼ 4 % for V 23 /∆ = 0.3. This indicates a significant robustness in the acceleration mechanism, yielding the robustness around ∼ 10 % with the inhomogeneous coupling upto 150%.

Discussion and Conclusion
Energy transport in molecular complexes / light harvesting complexes in photosynthetic bacteria has raised many fundamental questions on how our nature operates and our understanding of noise effects in such processes has been challenged. The recent research indicates that noise can enhance transport rates through temporal correlations [11][12][13][14][15] or spatial correlations [18]. In this work we have considered the transport of excitation for a multi-site linear chain model whose energy levels are affected by spatio-temporal correlated stochastic noise processes [16,17]. We find that the energy transport can be accelerated by extending the spatial correlation into the negative region (anti-correlation). In the two-site model, our numerical analysis showed the significant acceleration in the energy transfer with the negative spatial correlation. By extending the model to three sites, we numerically demonstrated that the optimal efficiencies can be obtained with the "anti-ferromagnetic" correlation. The difference in transport with anti-correlated and uncorrelated noises in the four site model becomes even larger. For the temporal correlation dependencies, we explored the energy transfer dependence by changing the parameter α = ∆ · τ c . A non-monotonic dependence on both the population time evolution and the transfer efficiency measures were observed, and hence the correlation time needs to be chosen to achieve the optimum energy transport. These results show new possibilities to understand efficient energy transport in nature and engineer it to our technologies.

Methods
An open systems approach: As shown in Fig (1) our multi-site model contains fluctuating excited state energy levels as well as an energy trap on the last site (energy decay). In principle this means an open systems approach must be used -especially as our energy level fluctuations are stochastic in nature. It would thus be natural to write down a master equation in Lindblad form [26,27] (which can easily handle the last site energy decay), which would force us to use a white noise model where the fluctuations are not correlated in time. However in this case we want to examine temporal and spatial correlation effects. This means the typical master equation is not appropiate, however a master equation can be derived using time-convolutionless decomposition techniques [28][29][30][31][32][33][34][35][36].

Derivation of the time convolution type of master equation:
The master equation given by Eq.(4) is obtained by extracting the time evolution of the excitations in each site from the one of the total system which is written by the Liouville-von Neuman equation as where W (t) is the density operator for the total system and L(t) is the Liouville operator defined as L(t)X = 1 h [H(t), X] for an arbitrary operator X. In such a case, the extraction averages W (t) over the stochastic process of the fluctuation. Our purpose is to obtain the time evolution of the reduced density operator W (t) (≡ ρ(t)) where we denote the average operation as · · · . For this purpose, it is convenient to use the projection operator method [28][29][30][31][32][33][34][35][36]. Introducing a projection operator, P, which is an idempotent operator with a property as P 2 = P, we describe the reduced density operator as PW (t) ≡ W (t) (= ρ(t)). To obtain the time evolution of PW (t), we use the formal solution of Eq. (7) as where T + is an increasing time ordering operator from the right to the left. The procedure to obtain the master equation, Eq.(4), is roughly divided into the following six steps: 1. First we define the relevant P part and the irrelevant Q(≡ 1 − P) part of the time evolution operatorÛ + (t, t 0 ) as where we use interaction picture with the definition asÛ + (t, t 0 ) = e −iL0(t−t0) U + (t, t 0 ) = T + exp[− t t0 dt iL 1 (t )dt ] andL 1 (t) = e iL0(t−t0) L 1 (t)e −iL0(t−t0) .
3. Next the formal solution of the irrelevant Q part can be written as, withV We then rewrite x(τ ) in Eq. (11) with x(t) and y(t) using the relation whereÛ − (t, t 0 ) = T − exp[i t t0L 1 (τ )dτ ] with T − an increasing time ordering operator from the left to the right. The formal solution of y(t) has the form where we define 5. Next we substitute the formal solution of y(t) into Eq. (9) and multiply W (t 0 ) from the right on the both hand sides of Eq. (9), d dtρ (t) = P(−iL 1 (t))ρ(t) + Ξ(t, t 0 )ρ(t) + I(t, t 0 )W (t 0 ), (15) withρ(t) = e −iL0(t−t0) ρ(t), Ξ(t, t 0 ) = P(−iL 1 (t))(1 − θ(t) −1 ) and I(t, t 0 ) = P(−iL 1 (t))θ(t) −1V (t, t 0 )Q. 6. Eq. (15) is then expanded with using the relation as θ(t) −1 = ∞ n=0 σ(t) n . This gives where ρ nm (t) is the (n, m) element of the reduced density operator ρ(t), κ is the trap frequency at the 2nd site and we F n (ω 1 , ω 2 , V 12 , t) for n = 1, 2 are defined as with In Eq. (20) The dynamics described by time convolutionless(TCL) type of master equation is compared to the one by hierarchical equations of motion (HEOM) for the spin-boson system in [37] where they find that the second order and fourth order of TCL equation and HEOM shows almost the same dynamics for weak coupling case.
Relation with the traditional Lindblad master equation: The TCL master equation given by Eq.(4) in a specific limit reduces to the typical Lindblad master equation. In such a case we set the upper integral limit t → ∞ giving where all of the coefficients are time-independent. This is the Born-Markov approximation [30,31,36]. For example, the time-dependent coefficients of the two site model ignoring spatial correlation (c = 0) can be approximated as Taking the limit of the correlation time τ c to approach zero with maintaining ∆ 2 τ c finite [31], we obtain Defining lim τc→0 F 2 (ω 1 , ω 2 , Γ 12 , ∞) ≡ γ, our master equation given by Eq.(4) reduces to the Lindblad form [26,27] : where A m = |m m|. We find that Eq. (24) is the same as the master equation obtained in [8] for the Haken-Strobl model. Moreover, when we substitute the stationary value of ρ 12 (t) into the differential equation, we obtain the same differential equation as in [17]. Besides the above, let us note that, by using the time convolutionless type of master equation, we need only to take long-time limit as the approximation procedure, which is much simpler than using the time-nonlocal (time-convolution) type of master equation [30]. Numerical methods: The effect of the fluctuation on the time evolution of the density operator is described with the time dependent coefficient of the second term in the right hand side of the time-convolutionless type of master equation given by Eq. (4). To obtain the time evolution of the density operator, we numerically solved the master equation by evaluating the time dependent coefficient and iterating the equation step by step with the coefficient. In the evaluation, we scaled the time variable and parameters with the strength of fluctuation ∆.

Quantum yield:
The quantum yield is defined as q = n k t,n τ n n k t,n τ n + n k d,n τ n , where k d,n is the the decay rate at the n-th site by recombination, k t,n is the one by trap and τ n = ∞ 0 ρ nn (t)dt in [17]. q indicates the ratio of trapped quantity to the total loss by recombination and trap. Cao and Silbey showed that the decay rate by recombination is necessary to be much smaller than the trapping rate, k d,n k t,n to obtain a high quantum yield [17]. In such situation, the dependence of τ n on k d,n can be neglected to give n k t,n τ n (k d,n = 0) ≈ 1. Thus q ≈ 1 1 + n k d,n τ n (k d,n = 0) where the last form is obtained by setting the values of k d,n to be the same as k d for all of the state n and with t = n τ n , which is called as the average trapping time.
Data availability: The data that support the findings of this study are available from the corresponding author upon reasonable request.

Aknowledgements
We would like to thank Ryuta Tezuka, Ryuji Nakamura, Shun Imazawa, and Dr. Kazunari Hashimoto for valuable discussions. This work was supported in part from the MEXT KAKENHI Grant-in-Aid for Scientific Research on Innovative Areas Science of hybrid quantum systems Grant No.15H05870, JSPS KAKENHI Grant Numbers 15K05200, and the researcher exchange promotion program of the Research Organization of Information and Systems. This project was also made possible through the support of a grant from the John Templeton Foundation. The opinions expressed in this publication are those of the authors and do not necessarily reflect the views of the John Templeton Foundation (JTF #60478).