A multi-site variational master equation approach to dissipative energy transfer

Unitary transformations can allow one to study open quantum systems in situations for which standard, weak-coupling type approximations are not valid. We develop here an extension of the variational (polaron) transformation approach to open system dynamics, which applies to arbitrarily large exciton transport networks with local environments. After deriving a time-local master equation in the transformed frame, we go on to compare the population dynamics predicted using our technique with other established master equations. The variational frame dynamics are found to agree with both weak coupling and full polaron master equations in their respective regions of validity. In parameter regimes considered difficult for these methods, the dynamics predicted by our technique are found to interpolate between the two. The variational method thus gives insight, across a broad range of parameters, into the competition between coherent and incoherent processes in determining the dynamical behaviour of energy transfer networks.


Introduction
The theory of open systems is necessary to describe any quantum system in contact with an uncontrollable and non-negligible environment. In problems of energy transport one is often interested in a regime where the environment, which consists of a very large number of degrees of freedom, is highly influential. We are here concerned with the dynamics of electronic excitations across some discrete network of molecules, in which the environment can play a key role -and this is exemplified by recent observations in various microscopic biological systems. Examples include: Molecular complexes involved in photosynthesis, such as the Fenna-Matthews-Olsen (FMO) complex in green sulphur bacteria [1,2,3,4,5,6,7,8,9,10,11] and the light harvesting photosystems in green plants [12]; the electron transport chain in Respiratory Complex I [13]; the donor-bridgeacceptor model of olfaction [14] and certain models of magnetoreception in birds [15,16]. Similar energy transport models have also been used extensively outside of biologyapplications range from the dynamics of coupled quantum dots [17,18,19,20] in solid state physics, to those of impurities in lattice Bose-Einstein condensates (BECs) [21].
Several techniques have been developed over the years to calculate the explicit time domain dynamics of open quantum systems. Some are numerically exact, meaning that given sufficient computational resources, they will converge to the correct dynamics under some well-controlled approximations. Such techniques include the path integral [22,23,24], hierarchical equations of motion (HEOM) [4,25,26] and density matrix renormalization group (DMRG) [27] methods. Though powerful, these approaches typically place restrictions on the kind of system that can be modelled, and they may also scale badly (in terms of computing resources) with the size and complexity of said system [4,23]. Often -as will be the case in this article -more numerically tractable, though approximate, methods are used. For example, those based on master equations [28,29,30]. This technique provides an equation of motion for the reduced density matrix of the system in question without having to track the full evolution of the environment, though normally involves some kind of perturbative expansion in a small parameter, such as the system-environment coupling strength. A further advantage of the master equation approach is that it can offer insights into the mechanisms underlying the dynamics of a system by relating rates and energy shifts directly to microscopic parameters. However, the obvious drawback of many master equations is that they rely on certain Hamiltonian parameters being small. If this condition is not fulfilled, then the truncation of the perturbative expansion often leads to (potentially unphysical) results which can diverge wildly from the true dynamics [31].
In certain parameter regimes, performing unitary transformations, such as the polaron transformation [32,33,34,18,35,36,37,38,39,40], on the combined systemenvironment Hamiltonian can result in a smaller interaction energy in the transformed frame. The transformed system is then amenable to being modelled using a perturbative master equation. For example, the polaron transformation can work well over a broad range of parameters when the relevant environmental timescales are short compared to those in the system -in fact, the polaron transformation diagonalizes the Hamiltonian we use below when no electronic couplings are present between the sites. It is thus often used when the coupling between system and environment is strong or when internal system couplings are small. Between the weak-coupling and polaron regimes, however, lies a region of parameter space for which neither model is appropriate. In addition, the polaron transformation runs into problems when applied to a system with an Ohmic or sub-Ohmic environment (one for which the environment spectral density scales linearly or sub-linearly, respectively, at low frequencies). In this case, infrared divergences arise which prevent certain master equations (such as the time-local form used in this work) from correctly predicting the dynamics in the transformed frame.
As an extension of the standard polaron approach, the variational (polaron) transformation [41,42,43,44,45,46,47] allows one to attempt to derive a perturbative series which is as valid as possible (given the restricted form of the transformation) in all parameter regimes. This is achieved by performing an optimized, partial polaron displacement on each of the environmental phonon modes relevant to their particular mode frequency, thereby interpolating between the weak-coupling and polaron representations for separate modes, as well as in the final master equation. Here we build on previous work on the variational transformation for two-site systems [43,44,46] which is in turn based on an idea originating with Silbey & Harris [41,42]. The major new contribution of this paper is the generalization of the formalism to any number of sites, allowing for the simulation of large networks across a range of environmental coupling parameters and temperatures. For comparison to other techniques, we have included examples of dynamics for systems in several different regimes.
In section 2 we describe the mathematical model for which the transformation is valid, including some of its limitations. We then go on to discuss the form of the variational polaron transformation and the accompanying optimization procedure in section 3. In section 4 we outline the master equation formalism in the variationally transformed frame, and in section 5 we present some example dynamics, including that of the FMO system. Finally, in section 6, we shall conclude by briefly discussing the various advantages and drawbacks of the method outlined herein.

The transport model
The system (S) considered in this work is that of N coupled two-level systems, known as sites. Between them they carry exactly one excitation -for molecular networks these are electronic excitations which, for charge neutral systems, are called excitons. The latter restriction to a single excitation allows for a great reduction in the size of the system Hilbert space, and is sufficient to describe the behaviour of many physically and biologically relevant systems [18,19,48]. For example, it is thought to be a valid approximation to the in vivo dynamics of the FMO complex studied in section 5 [5]. Each of the N sites is linearly-coupled to its own, independent phonon environment (E). The Hamiltonian for the combined system and environment is given by where b n,k is the annihilation operator for phonon mode k on site n, and |n is the state of S in which only site n is excited (see figure 1 for a cartoon visualization of the above Hamiltonian). The environmental spectrum is usually taken to be a continuum, such that the couplings g n,k can be described in terms of a spectral density J n (ω). The spectral density is defined by J n (ω) = k |g n,k | 2 δ(ω − ω k ); (2) which takes into account the density of states, dispersion relation and interaction mechanism with the environment. In the continuum case, a good measure of the strength of the system-environment coupling at each site is the reorganization energy [25]: This model ignores any spatial correlations between phonon excitations at different sites, meaning that the Hamiltonian in 1 is not relevant for systems with strong, long-range correlations, such as impurities in BECs. For the case of FMO it has been claimed, based on detailed molecular dynamics simulations, that spatial correlations do not play a significant role in the exciton dynamics [49]. Since the environments at each site are independent, one can take the couplings to be real for this model without loss of generality. In the case of a global environment, however, the phases of the couplings to each site can encode the environmental correlations between them. We assume that the environment is initially in a thermal (Gibbs) state at temperature T = 1/(k B β): ρ E (0) = e −βH E / tr(e −βH E ). In addition, the combined system-environment state is assumed to be initially separable, such that there are no system-environment correlations:

The variational polaron transformation
The Hamiltonian in 1 describes the transport problem in a way which is intuitive and transparent, in that H S and H E are the Hamiltonians for the system and environment in isolation and the interaction term H I is simple in form. However, one of the features of quantum mechanics is that the physics of a system -even a composite one -is invariant under unitary transformations. As a consequence, one is not restricted to a single way of distinguishing two subsystems. In the model outlined above, much of the energy associated with the interaction between system and environment is due to the excitation deforming the surrounding molecular structure, and hence affecting the state of the phonon environment.
Applying the polaron transformation allows us to move into a reference frame where this back-action from system to environment is accounted for at the Hamiltonian level. The system is 'dressed' by the environment, and the environmental phonon modes are displaced in phase space conditional on the state of the system. The explicit form of the transformation is leading to a transformed Hamiltonian: The interaction HamiltonianH I now contains two terms. One,H L , is of the same linear form as the interaction in the untransformed Hamiltonian, albeit with modified coupling strength. The other term,H D , contains a new kind of interaction between off-diagonal system operators and the environmental displacement operators Here, the expectation values (B n B m ) of the displacement operators in the B nm have been taken into the system Hamiltonian, and are thus treated as renormalized couplings between the sites. The B n 's are given by for a thermal equilibrium environmental state ρ E . The site energies after transformation are also shifted in comparison to the original frame by a factor R n , defined as Usually when the polaron transformation is discussed in the literature, what is meant is the fully-displaced version of 4 [39,40] where f n,k = g n,k , for all n and k. This results inH L = 0, which leaves only the new displacement interaction term in the Hamiltonian:H I =H D . In the variational case, however, the f n,k are left as free parameters, and minimization over an upper bound on the free energy, as described below, determines their values. The idea is that the optimization inherent to the variational approach allows us to minimize the effect of the interaction Hamiltoniañ H I , given the transformation form. This is done here in order to validate the use of perturbation series in various master equation approaches, which must in practice be truncated at some finite order. In general, the 'smaller'H I , the more accurate the low-order dynamics are likely to be.
Since, most of the time, there is no single parameter inH I which determines exactly how small its effect is, we choose instead to optimize the variational transformation by minimizing the contribution ofH I to the free energy (the average energy of a thermal state of the system). This choice is consistent with earlier variational treatments [41,42] and ensures that the steady state of the resulting dynamics is as accurate as possible -in equilibrium the free energy should be at a minimum. As it is generally impossible to find an exact analytical expression for the free energy, it is the Feynman-Bogoliubov upper bound [50] that we shall minimize. The bound is given by where X H 0 = tr Xe −βH 0 . The true free energy A is related to this bound by the inequality A ≤ A B . Given that we want to end up withH I small, it is reasonable to neglect the higher order terms in 8 as a first approximation. Furthermore, the interaction Hamiltonian in the transformed frame has been constructed such that the second term goes to zero, H I H 0 = 0. Therefore, minimization amounts to maximizing the value of tr(e −βH 0 ). Although, perhaps counter-intuitively,H I now appears to be absent from A B its influence is, in fact, still present implicitly inH 0 .
The transformed system Hamiltonian can be written as a function of the renormalization parameters {R n , B n }, therefore the minimization condition can be written: which, after using the expressions for the renormalization parameters in 7 allows us to write f n,k = F n (ω n,k )g n,k , with In the continuum limit for the environment, the minimization procedure for an N-site system amounts to solving the 2N coupled integral equations given by the definitions of the renormalization parameters:

Master equation formulation in the variational frame
As outlined above, a truncated perturbative expansion in the new interaction HamiltonianH I should, following the optimization procedure, be as accurate as possible, given the polaron form of the transformation and the minimization condition used. The next step, therefore, is to derive a master equation in the variational frame using standard techniques. By utilizing a projection operator P with the following action on the combined system-environment state: Pρ = tr E (ρ) ⊗ ρ R , where ρ R is an arbitrary reference state for the environment, we can separate out the reduced system dynamics.
Here, we choose to derive a time-local or 'time-convolutionless' master equation due to the relative ease with which it can be solved numerically. In addition, we cut off the perturbation series at second order. The resulting master equation has the following general form in the interaction picture [28]: where K 2 and I 2 are superoperators acting on Pρ and (1 − P)ρ respectively, which have been curtailed to second order inH I . In the untransformed frame the separable initial state means that the second, inhomogeneous term in 13 disappears for the choice ρ R = ρ E (0). In the variational frame this is no longer the case and the inhomogeneous term must be taken into account. However, for two-site systems the inhomogeneous term was seen to have only a small, transient effect on the dynamics at finite temperatures [44] for single site initial excitations. Therefore, we shall henceforth neglect it even in the transformed frame. This amounts to assuming that the environment relaxes into its displaced state instantaneously. One would expect this to be a good approximation at finite temperatures for smooth spectral densities and when the typical environment timescales are shorter than the transition timescales in the system [37,38,39,40,44]. The examples we present in section 5 satisfy each of these conditions.
The remaining (homogeneous) term is written explicitly as By writing the interaction Hamiltonian in the formH I = we can rewrite the master equation in terms of system operators S i and two-time environmental correlation func- After moving back into the Schrödinger picture, the master equation takes the form: The interaction Hamiltonian system operators can be split into three distinct groups in the following way: which leads in turn to three varieties of non-zero time correlation function. The first type are due to the linear interaction term,H L , and are therefore of the same form as those that appear in the standard weak coupling master equation: (17) where F n (ω) is the continuum version of the optimized function in 10. The second type come from the displacement operator interaction,H D , and are the only type to appear in the fully displaced polaron master equation: where φ xy and the δ nm are Kronecker deltas. Finally, the third type appear in the more general variational master equation due to an overlap between the two types of interaction: with φ yz The dynamics calculated using 15 will be of the density matrix in the variationally transformed frame. In order to consider quantities in the original frame, one must perform the inverse of the transformation in 4. The site populations (diagonal elements of the density matrix) are unchanged, since the operators |n n| commute with the transformation. However, in general, the inverse transformations of the coherences (offdiagonal elements of the density matrix) are much more difficult to calculate. In the case that the inhomogeneous term in 15 is ignored, one can make the approximation [40] (ρ S (t)) nm = B n B m (ρ S (t)) nm for n = m, whereρ S (s) is the system density matrix in the variational frame and ρ S (s) is the system density matrix in the untransformed frame. This is equivalent to making a Born approximation in the transformed frame. Under such an approximation the transformed frame state isρ(t) ≃ρ S (t) ⊗ ρ R for all times. This state will transform in the same way under the inverse variational transformation as the system Hamiltonian 1 does under the forward transformation, leading to the factors of B n B m for the off-diagonal elements mentioned above.

Three sites
We would like to compare dynamics calculated in the variational frame, using 15, with other techniques. After two sites, the next simplest system with a Hamiltonian of the form of that in 1 has three sites and only nearest neighbour couplings. Figure 2 shows the dynamics for such a system in a variety of parameter regimes, calculated in the variational frame. For comparison, we have also plotted the dynamics calculated in the fully-displaced polaron frame as well as that calculated using the untransformed Hamiltonian (weak-coupling, or Redfield, approximation). The system is characterized by its on-site energies {E n }, inter-site couplings V 12 and V 23 , and spectral densities of the form Column (a) in figure 2 represents a regime where system frequencies (∼ 20cm −1 ) are much smaller than the environment cutoff frequency (200cm −1 ). In this case, the fully-displaced polaron transformation is expected to do well [45], and we would also expect the variational transformation to match it, as it indeed does. The untransformed dynamics fail to reach the correct steady state due to a reasonably large reorganization energy for the environment. We also found that the weak coupling approximation can lead to unphysical results for parameters differing from those in column (a) only by their intersite coupling. This was not the case for the dynamics in the variational or full polaron frames. The second column, (b), shows dynamics in a regime for which neither the weakcoupling nor full-polaron transformation are ideally suited. System frequencies are comparable to environment frequencies, and the coupling to the environment is not small. The variational dynamics appear to interpolate between the two other results. It is clear that it agrees with the weak-coupling dynamics at short times, before settling on a different set of long time populations. The comparatively large system energies prevent the polaron transformation from dealing correctly with the lower frequency parts of the environment, leading to the dynamics in the full-polaron frame overestimating the damping of coherent oscillations. The variational transformation preserves coherence here precisely due to the fact that it optimizes the frequency dependence of the polaron transformation, as opposed to indiscriminately displacing every phonon mode by the full amount.
Finally, in column (c) the reorganization energy is much lower than in column (b), but the system frequencies are still large. A weak coupling approximation is therefore valid in this case, and one would expect the dynamics in the untransformed frame to be more accurate than that in the (full) polaron frame. The red and black curves in this panel almost sit on top of each other, showing that the variational dynamics agree with the weak coupling results, and corroborates the fact that the variational transformation allows us to capture the dynamics across a broad range of coupling strengths. Note that, interestingly, both the weak-coupling and polaron approaches can overestimate the damping of coherence in comparison to the variational method, dependent on the parameter regime (cf. panels (a) and (c)).

The Fenna-Matthews-Olsen complex
We now analyse how the variational master equation performs for a larger system, namely the FMO complex. This is usually assumed to be a seven site system (although recent results suggest there is in fact an additional eighth site [8]) and thus has a much larger parameter space than the three site system in figure 2. Despite this, one can see from figure 3, which compares FMO dynamics across a range of reorganization energies and temperatures, that the variational transformation performs the same kind of interpolation between weak-coupling and full-polaron dynamics as in the three-site case. The FMO system Hamiltonian used in this section -taken from [51] -is: and the spectral density -the smooth part of that from [51] -is of the form: with ω 1 = 0.575cm −1 and ω 2 = 2cm −1 . Figure 3 provides a clear example of the importance of the interplay between coherent and incoherent dynamics in excitonic energy transport. The biological purpose  Figure 3. Variational (solid), weak-coupling (dashed) and full-polaron (dotted) dynamics for the populations of site 1 (black, red, blue) and site 3 (olive, orange, cyan) of the FMO complex. Although all seven sites were modelled, only the input (1) and output (3) site populations are shown for clarity. The system Hamiltonian used is the same as that in [51] and the spectral density is given in 24. Panels (a-c) have η = 1 2 , (d-f) have η = 1, (g-i) have η = 2, and (j-l) have η = 4. The dynamics in the first column of plots were calculated at 5 K, the second a 77 K, and the third at 300 K.
of the FMO complex is to transport excitations from site 1 (or sometimes site 6) to site 3, from which the excitation is then removed [3]. It is therefore beneficial to have the population on site 3 build up as fast as possible. One can see from the figure (most clearly in the second column), that the optimum rate of transfer (panels (j), (h) and (i)) occurs in the variational theory when coupling to the environment is neither too strong nor too weak. That is, phonon-assisted transport is enhanced by the presence of some degree of coherence in the system. These optimal cases appear to lie in the intermediate region of parameter space, outside the remit of weak-coupling or polaron master equations, for which something like the variational approach is required.
The full, seven-site dynamics for the cases where η = 1, at T = 77 K and T = 300 K, are shown in figure 4. These plots correspond to the parameter regimes of the FMO  [51], with the exception that only the smooth, non-peaked part of the spectral density was used. The dynamics in panel (a) were calculated at 300 K whilst those in (b) were calculated at 77 K. Two different initial states: ρ S (0) = |1 1| and ρ S (0) = |6 6| were used for the upper and lower plots, respectively. dynamics in figure 2 of [51]. However, the exact calculations presented in that paper include a significant peak in the spectral density which, in any master equation approach, would ideally be treated separately from the rest of the environment in order to capture its effect on the dynamics non-perturbatively. Moreover, a variational polaron treatment of such a peak would likely cause there to be significant system-environment correlations in the transformed frame, which would not be taken into account without the inclusion of inhomogeneous terms in the master equation. That being said, the qualitative agreement of the variational dynamics shown in the figure with the results in [51] is surprisingly good.

Discussion
There are several advantages to using master equations over other approaches to open quantum systems dynamics. Primarily, these are the efficiency with which one can solve them, and the potential insights into underlying physics which they can give. We can see from the various terms in 15 exactly how the parameters in the Hamiltonian enter into the system dynamics and the relative magnitude of these terms can give us an idea of the parameter regime of a given system -in the sense of which quantities are contributing most to the dynamics.
Moving into the variational frame prior to solving the master equation allows us to calculate sensible dynamics over a larger range of parameters compared to the more standard weak-coupling or full-polaron approaches. As can be seen from section 5, the variational master equation can capture the dynamics in both the weak and strong coupling regimes, when the weak coupling and polaron master equations, respectively, are expected to do well. It is also able to bridge the gap between the two in intermediate regimes.
The variational master equation is expected to work well over a wide range of parameters, and can in principle handle arbitrary spectral densities. However, like the full polaron transformation, it works best in the scaling limit [45], in which important environmental frequencies (ω c ) are large compared to relevant system frequencies (V nm ). The corresponding downside is that moving into the variational frame provides less of an advantage in terms of improving the accuracy of the dynamics when the typical environmental timescales are significantly longer than those of the system, and the coupling between the two is strong. The intuition for this is that the low frequency phonon modes are too 'sluggish' to keep up with the motion of the exciton as it moves through the system and do not, therefore, dress the system in the same way as higher frequency modes. They may still, however, have a profound impact on the system dynamics, which a transfomation of displacement form is unable to capture.
One aspect of our method which might benefit from modification is the specific minimization condition used. Whilst we expect the first term of 8 to be a good metric for the size of the interaction, it does not directly correspond to the quantity which we expand perturbatively in the master equation. In fact, the next highest order term, (O(H 2 I )), is mathematically more similar [52], albeit much more complicated. As one smoothly varies the Hamiltonian parameters, the optimum transformation can jump, as different local minima become global minima. This effect has been studied for the case of two sites [47], and it was found that the variational predictions are less accurate around the discontinuity. For multiple sites, the free energy landscape becomes more complex and more local minima emerge in parameter space, hence there is greater scope for this kind of jumping. Whilst not optimal, the transformations corresponding to such local minima are still likely to lead to more accurate dynamics than those calculated in the untransformed frame.
In summary, we have outlined a variational method for solving open quantum systems dynamics in molecular networks with local environments. The method is valid over a wide range of parameters and is efficient to compute. By moving into a reference frame in which system and environment are less strongly interacting, one is able to use a perturbative master equation to more accurately calculate dynamics. For our Hamiltonian 1, the approach surpasses both weak-coupling and full-polaron master equations in terms of breadth of applicability. It can be used to model interesting biological systems which sit in difficult intermediate coupling regimes, such as FMO, and allows for the systematic study of the effects of certain parameters on the dynamics.
There is still ample room for improvement, and the general concept of redrawing the boundary between system and environment has far greater reach than the implementation presented in this paper. For instance, one could generalize the transformation to a larger class of Hamiltonians, say those whose environments couple to multiple sites, or one could augment the form of the transformation itself, perhaps by including squeezing in addition to displacement. The technique we have developed utilizes one of the most fundamental properties of quantum mechanics, namely the invariance of dynamical laws under unitary transformations, and gives insight into the important physical mechanisms underlying the evolution of open quantum systems. Whilst future master equation approaches may go beyond the polaron transformation, they are likely to benefit from a kind of variational minimization in the spirit of that which we have outlined here.