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

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–11] and the light harvesting photosystems in green plants [12]; the electron transport chain in respiratory complex I [13]; the donor-bridge-acceptor model of olfaction [14] and certain models of magnetoreception in birds [15, 16]. Similar energy transport models have also been used extensively outside of biology—applications range from the dynamics of coupled quantum dots [17–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–24], hierarchical equations of motion [4, 25, 26] and density matrix renormalization group [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 paper—more numerically tractable, though approximate, methods are used. For example, those based on master equations [28–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 [18, 32–40], on the combined system–environment 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–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 and 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 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

$$\begin{eqnarray} &&H = H_{\mathrm{S}} + H_{\mathrm{E}} + H_{\mathrm{I}}, \quad H_{\mathrm{S}} = \sum_{n}{\varepsilon_n \mbox{$|n\rangle\langle n|$}} + \sum_{n \neq m}{V_{nm} \mbox{$|n\rangle\langle m|$}}, \nonumber \\ &&H_{\mathrm{E}} = \sum_{n,\mathbf{k}}{ \omega_{n,\mathbf{k}} b^\dagger_{n,\mathbf{k}} b_{n,\mathbf{k}}}, \quad H_{\mathrm{I}} = \sum_{n,\mathbf{k}}{\mbox{$|n\rangle\langle n|$}(g_{n,\mathbf{k}} b^\dagger_{n,\mathbf{k}} + g^*_{n,\mathbf{k}} b_{n,\mathbf{k}})}, \end{eqnarray} \tag{ 1 }$$

where b_n,k is the annihilation operator for phonon mode k on site n and $\left | n \right >$ is the state of S in which only site n is excited (see figure 1 for a cartoon visualization of the above Hamiltonian).

Zoom In Zoom Out Reset image size

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

$$\begin{equation} J_n(\omega) = \sum_{\mathbf{k}}{|g_{n,\mathbf{k}}|^2 \delta(\omega-\omega_k),} \end{equation} \tag{ 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]:

$$\begin{equation} \lambda_n = \int_0^\infty{{\rm d}\omega \frac{J_n(\omega)}{\omega}}. \end{equation} \tag{ 3 }$$

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: ρ(0) = ρ_S(0)⊗ρ_E(0).

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

$$\begin{equation} \tilde{H} = {\rm e}^G H {\rm e}^{-G}, \quad {\rm where} \quad G = \sum_{n,\mathbf{k}}{\mbox{$|n\rangle\langle n|$} \omega_{n,\mathbf{k}}^{-1} (f_{n,\mathbf{k}} b^\dagger_{n,\mathbf{k}} - f^*_{n,\mathbf{k}} b_{n,\mathbf{k}})}, \end{equation} \tag{ 4 }$$

leading to a transformed Hamiltonian:

$$\begin{eqnarray} \tilde{H} &=& \tilde{H}_0 + \tilde{H}_{\mathrm{I}}, \nonumber \\ \tilde{H}_0 &=& \tilde{H}_{\mathrm{S}} + \tilde{H}_{\mathrm{E}}, \quad \tilde{H}_{\mathrm{I}} = \tilde{H}_{\mathrm{L}} + \tilde{H}_{\mathrm{D}},\\ \tilde{H}_{\mathrm{S}} &=& \sum_{n}{(\varepsilon_n + R_n) \mbox{$|n\rangle\langle n|$}} + \sum_{n \neq m}{\mathcal{B}_{n}\mathcal{B}_{m} V_{nm} \mbox{$|n\rangle\langle m|$}}, \nonumber \end{eqnarray} \tag{ 5 }$$

$$\begin{eqnarray} \tilde{H}_{\mathrm{L}} &=& \sum_{n,\mathbf{k}}{\mbox{$|n\rangle\langle n|$}[(g_{n,\mathbf{k}}-f_{n,\mathbf{k}}) b^\dagger_{n,\mathbf{k}} + (g_{n,\mathbf{k}}-f_{n,\mathbf{k}})^* b_{n,\mathbf{k}}]},\nonumber\\ \tilde{H}_{\mathrm{D}}&=& \sum_{n \neq m}{ V_{nm} \mbox{$|n\rangle\langle m|$} B_{nm}}, \quad \tilde{H}_{\mathrm{E}} = H_{\mathrm{E}}. \end{eqnarray} \tag{ 6 }$$

The interaction Hamiltonian $\tilde {H}_{\mathrm {I}}$ now contains two terms. One, $\tilde {H}_{\mathrm {L}}$, is of the same linear form as the interaction in the untransformed Hamiltonian, albeit with modified coupling strength. The other term, $\tilde {H}_{\mathrm {D}}$, contains a new kind of interaction between off-diagonal system operators and the environmental displacement operators

$$\begin{eqnarray*} &&\fl B_{nm} = B_{mn}^\dagger = \exp\left[\sum_{\mathbf{k}}{\omega_{n,\mathbf{k}}^{-1} (f_{n,\mathbf{k}} b^\dagger_{n,\mathbf{k}} - f^*_{n,\mathbf{k}} b_{n,\mathbf{k}})} - \sum_{\mathbf{k}}{\omega_{m,\mathbf{k}}^{-1} (f_{m,\mathbf{k}} b^\dagger_{m,\mathbf{k}} - f^*_{m,\mathbf{k}} b_{m,\mathbf{k}})}\right] - \mathcal{B}_{n}\mathcal{B}_{m}. \end{eqnarray*}$$

Here, the expectation values ($\mathcal {B}_{n}\mathcal {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 $\mathcal {B}_n$ are given by

$$\begin{eqnarray*} \mathcal{B}_n &=& {\rm tr}\left\{ \exp\left[\sum_{\mathbf{k}}{\omega_{n,\mathbf{k}}^{-1} (f_{n,\mathbf{k}} b^\dagger_{n,\mathbf{k}} - f^*_{n,\mathbf{k}} b_{n,\mathbf{k}})}\right] \rho_{\mathrm{E}} \right\}\\ &=& \exp\left[-\frac{1}{2} \sum_{\mathbf{k}}{\frac{|f_{n,\mathbf{k}}|^2}{\omega_{n,\mathbf{k}}^{2}} \coth(\beta \omega_{n,\mathbf{k}}/2)}\right] \end{eqnarray*}$$

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

$$\begin{equation} R_n = \sum_{\mathbf{k}}{\omega_{n,\mathbf{k}}^{-1}\left[|f_{n,\mathbf{k}}|^2-2\Re(f_{n,\mathbf{k}} g_{n,\mathbf{k}}^*)\right]}. \end{equation} \tag{ 7 }$$

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 in $\tilde {H}_{\mathrm {L}} = 0$, which leaves only the new displacement interaction term in the Hamiltonian: $\tilde {H}_{\mathrm {I}} = \tilde {H}_{\mathrm {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 Hamiltonian $\tilde {H}_{\mathrm {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' $\tilde {H}_{\mathrm {I}}$, the more accurate the low-order dynamics are likely to be.

Since, most of the time, there is no single parameter in $\tilde {H}_{\mathrm {I}}$ which determines exactly how small its effect is, we choose instead to optimize the variational transformation by minimizing the contribution of $\tilde {H}_{\mathrm {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

$$\begin{equation} A_{\mathrm{B}} = -\frac{1}{\beta} \ln [{\rm tr}({\rm e}^{-\beta \tilde{H}_0})] + \langle \tilde{H}_{\mathrm{I}} \rangle_{\tilde{H}_0} + \mathcal{O}(\langle \tilde{H}_{\mathrm{I}}^2 \rangle_{\tilde{H}_0}), \end{equation} \tag{ 8 }$$

where $\langle X \rangle _{\tilde {H}_0} = {\rm tr}(X {\mathrm { e}}^{-\beta \tilde {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 with $\tilde {H}_{\mathrm {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, $\langle \tilde {H}_{\mathrm {I}} \rangle _{\tilde {H}_0} = 0$. Therefore, minimization amounts to maximizing the value of ${\rm tr}({\mathrm { e}}^{-\beta \tilde {H}_0})$. Although, perhaps counter-intuitively, $\tilde {H}_{\mathrm {I}}$ now appears to be absent from A_B its influence is, in fact, still present implicitly in $\tilde {H}_0$.

The transformed system Hamiltonian can be written as a function of the renormalization parameters $\{ R_n, \mathcal {B}_n \}$, therefore the minimization condition can be written

$$\begin{equation} \frac{\partial A_{\mathrm{B}}}{\partial f_{n,\mathbf{k}}} = \frac{\partial A_{\mathrm{B}}}{\partial R_n}\frac{\partial R_n}{\partial f_{n,\mathbf{k}}} + \frac{\partial A_{\mathrm{B}}}{\partial \mathcal{B}_n}\frac{\partial \mathcal{B}_n}{\partial f_{n,\mathbf{k}}} = 0 \quad \forall \; n, \mathbf{k}, \end{equation} \tag{ 9 }$$

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

$$\begin{equation} F_n(\omega_{n,\mathbf{k}},\{ R_n, \mathcal{B}_n \}) = \frac{2 \omega_{n,\mathbf{k}} \frac{\partial A_{\mathrm{B}}}{\partial R_n}}{2 \omega_{n,\mathbf{k}} \frac{\partial A_{\mathrm{B}}}{\partial R_n}- \mathcal{B}_n \frac{\partial A_{\mathrm{B}}}{\partial \mathcal{B}_n} \coth(\beta \omega_{n,\mathbf{k}} / 2)}. \end{equation} \tag{ 10 }$$

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

$$\begin{equation} \mathcal{B}_n = \exp\left[-\frac{1}{2} \int_0^\infty{{\rm d}\omega \; \frac{F_n(\omega,\{ R_n, \mathcal{B}_n \})^2 J_n(\omega)}{\omega^{2}} \coth(\beta \omega/2)}\right], \end{equation} \tag{ 11 }$$

$$\begin{equation} R_n = \int_0^\infty{{\rm d}\omega \; \frac{J_n(\omega)}{\omega} F_n(\omega,\{ R_n, \mathcal{B}_n \})\left[ F_n(\omega,\{ R_n, \mathcal{B}_n \})-2 \right]}. \end{equation} \tag{ 12 }$$

As outlined above, a truncated perturbative expansion in the new interaction Hamiltonian $\tilde {H}_{\mathrm {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 $\mathcal {P}$ with the following action on the combined system–environment state: $\mathcal {P} \tilde {\rho } = {\rm tr}_{\mathrm {E}}(\tilde {\rho }) \otimes \rho _{\mathrm {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]:

$$\begin{equation} \frac{\partial}{\partial t} {\rm tr}_{\mathrm{E}}\left\{\mathcal{P} \tilde{\rho}(t)\right\} = {\rm tr}_{\mathrm{E}}\left\{\mathcal{K}_2(t) \mathcal{P} \tilde{\rho}(t)\right\} + {\rm tr}_{\mathrm{E}}\left\{\mathcal{I}_2(t) (1-\mathcal{P}) \tilde{\rho}(0)\right\}, \end{equation} \tag{ 13 }$$

where $\mathcal {K}_2$ and $\mathcal {I}_2$ are superoperators acting on $\mathcal {P} \tilde {\rho }$ and $(1-\mathcal {P}) \tilde {\rho }$ respectively, which have been curtailed to second order in $\tilde {H}_{\mathrm {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–40, 44]. The examples we present in section 5 satisfy each of these conditions.

The remaining (homogeneous) term is written explicitly as

$$\begin{equation} {\rm tr}_{\mathrm{E}}\left\{\mathcal{K}_2(t) \mathcal{P} \tilde{\rho}(t)\right\} = - \int_0^t{{\rm d} s \; {\rm tr}_{\mathrm{E}}\{[\tilde{H}_{\mathrm{I}}(t),[\tilde{H}_{\mathrm{I}}(s),\mathcal{P} \tilde{\rho}(t)]]\}}. \end{equation} \tag{ 14 }$$

By writing the interaction Hamiltonian in the form $\tilde {H}_{\mathrm {I}} = \sum _{i=1}^{N^2}S_i \otimes E_i$ (with interaction picture counterpart $\tilde {H}_{\mathrm {I}}(t) = \sum _{i=1}^{N^2}S_i(t) \otimes E_i(t)$) we can rewrite the master equation in terms of system operators S_i and two-time environmental correlation functions $\Lambda _{ij}(t-s) = {\rm tr}_{\mathrm {E}}\left \{E_i(t)E_j(s) \rho _{\mathrm {R}}\right \}$. After moving back into the Schrödinger picture, the master equation takes the form

$$\begin{eqnarray} &&\fl \frac{\partial \tilde{\rho}_{\mathrm{S}}(t)}{\partial t} = - {\rm i} [\tilde{H}_{\mathrm{S}},\tilde{\rho}_{\mathrm{S}}(t)] \nonumber\\ &&- \sum_{i,j} \int_0^t{{\rm d} s \; (\Lambda_{ij}(s)\{S_i S_j(s) \tilde{\rho}_{\mathrm{S}}(t) - S_j(s) \tilde{\rho}_{\mathrm{S}}(t) S_i\} + {\rm h.c.})}. \end{eqnarray} \tag{ 15 }$$

The interaction Hamiltonian system operators can be split into three distinct groups in the following way:

$$\begin{eqnarray} S_i = \left\{\begin{array}{@{}lrclr@{}} \mbox{$|n\rangle\langle n|$} = S^z_n, \quad & 1 \leqslant & {i} & \leqslant N, & \\ \mbox{$|n\rangle\langle m|$} + \mbox{$|m\rangle\langle n|$} = S^x_{nm}, \quad & N < & {i} & \leqslant \frac{1}{2}N(N+1), & \\ i\mbox{$|n\rangle\langle m|$} - {\rm i}\mbox{$|m\rangle\langle n|$} = S^y_{nm}, \quad & \frac{1}{2}N(N+1) < & {i} & \leqslant N^2 & \end{array} \right. \end{eqnarray} \tag{ 16 }$$

which leads in turn to three varieties of non-zero time correlation function. The first type are due to the linear interaction term, $\tilde {H}_{\mathrm {L}}$, and are therefore of the same form as those that appear in the standard weak coupling master equation:

$$\begin{eqnarray} \Lambda^{zz}_{n}(t) &=& \phi^{zz}_n(t) \nonumber \\ &=& \int_0^\infty{{\rm d}\omega \; J_n(\omega)\left[1-F_n(\omega)\right]^2\left[\cos(\omega t) \coth\,(\beta \omega /2) - {\rm i} \sin\,(\omega t)\right]}, \end{eqnarray} \tag{ 17 }$$

where F_n(ω) is the continuum version of the optimized function in (10). The second type come from the displacement operator interaction, $\tilde {H}_{\mathrm {D}}$, and are the only type to appear in the fully displaced polaron master equation:

$$\begin{eqnarray*} &&\fl \Lambda^{xx}_{nmpq}(t) = \frac{1}{2}V_{nm}V_{pq}\mathcal{B}_n \mathcal{B}_m \mathcal{B}_p \mathcal{B}_q \left\{ \exp\left[\delta_{np} \phi^{xy}_n(t) +\delta_{mq} \phi^{xy}_m(t) \right]\right.\\ &&\left.+ \exp \left[-\delta_{np} \phi^{xy}_n(t) -\delta_{mq} \phi^{xy}_m(t) \right]-2\right\}, \nonumber \end{eqnarray*}$$

$$\begin{eqnarray} &&\fl \Lambda^{yy}_{nmpq}(t)= \frac{1}{2}V_{nm}V_{pq}\mathcal{B}_n \mathcal{B}_m \mathcal{B}_p \mathcal{B}_q \left\{ \exp\left[\delta_{np} \phi^{xy}_n(t) +(\delta_{mq}-\delta_{mp}) \phi^{xy}_m(t) \right]\right. \nonumber\\ &&\left.- \exp\left[-\delta_{np}\phi^{xy}_n(t) -(\delta_{mq}-\delta_{mp})\phi^{xy}_m(t) \right]\right\}, \end{eqnarray} \tag{ 18 }$$

where

$$\begin{equation} \phi^{xy}_n(t) = \int_0^\infty{{\rm d}\omega \; \frac{J_n(\omega)}{\omega^2} F_n(\omega)^2 \left[\cos(\omega t) \coth(\beta \omega /2) - {\rm i} \sin(\omega t)\right]} \end{equation} \tag{ 19 }$$

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:

$$\begin{equation} \Lambda^{yz}_{nmp}(t) = \delta_{np} V_{nm}\mathcal{B}_n \mathcal{B}_m \phi^{yz}_n(t) \end{equation} \tag{ 20 }$$

with

$$\begin{equation} \fl \phi^{yz}_n(t) = \int_0^\infty{\rm d}\omega \; \frac{J_n(\omega)}{\omega} F_n(\omega)\left[1-F_n(\omega)\right] \left[\sin\,(\omega t) \coth(\beta \omega /2) + {\rm i} \cos\,(\omega t)\right]. \end{equation} \tag{ 21 }$$

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 (off-diagonal 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] $(\rho _{\mathrm {S}}(t))_{nm} = \mathcal {B}_n\mathcal {B}_m (\tilde {\rho }_{\mathrm {S}}(t))_{nm}$ for n ≠ m, where $\tilde {\rho}_{\rm S}(t)$ is the system density matrix in the variational frame and ρ_S(t) 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 $\tilde {\rho }(t)\simeq \tilde {\rho }_{\mathrm {S}}(t)\otimes \rho _{\mathrm {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 $\mathcal {B}_n\mathcal {B}_m$ for the off-diagonal elements mentioned above.

5.1. 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₁₂ and V₂₃ and spectral densities of the form

$$\begin{equation} J_n(\omega) = \frac{\lambda_n}{2}\frac{\omega^3}{\omega_{\mathrm{c}}^3}{\rm e}^{-\omega/\omega_{\mathrm{c}}}. \end{equation} \tag{ 22 }$$

Column (a) in figure 2 represents a regime where system frequencies (∼20 cm⁻¹) are much smaller than the environment cutoff frequency (200 cm⁻¹). 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.

Zoom In Zoom Out Reset image size

The second column, (b), shows dynamics in a regime for which neither the weak-coupling 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)).

5.2. 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

$$\begin{equation} H_{\mathrm{FMO}}\!=\!\left(\!\kern-2pt\begin{array}{@{}ccccccc@{}} 240 & -87.7 & 5.5 & -5.9 & 6.7 & -13.7 & -9.9 \\ -87.7 & 315 & 30.8 & 8.2 & 0.7 & 11.8 & 4.3 \\ 5.5 & 30.8 & 0 & -53.5 & -2.2 & -9.6 & 6.0 \\ -5.9 & 8.2 & -53.5 & 130 & -70.7 & -17.0 & -63.3 \\ 6.7 & 0.7 & -2.2 & -70.7 & 285 & 81.1 & -1.3\\ -13.7 & 11.8 & -9.6 & -17.0 & 81.1 & 435 & 39.7\\ -9.9 & 4.3 & 6.0 & -63.3 & -1.3 & 39.7 & 245 \end{array}\!\right) {\rm cm}^{-1}\quad \end{equation} \tag{ 23 }$$

and the spectral density—the smooth part of that from [51]—is of the form

$$\begin{equation} J_n(\omega) = 3.053 \times 10^{-5} \times \eta\frac{\omega^5}{\omega_1^4}{\rm e}^{-\sqrt{\omega/\omega_1}}+1.908 \times 10^{-5} \times \eta\frac{\omega^5}{\omega_2^4}{\rm e}^{-\sqrt{\omega/\omega_2}} \end{equation} \tag{ 24 }$$

with ω₁ = 0.575 cm⁻¹ and ω₂ = 2 cm⁻¹.

Zoom In Zoom Out Reset image size

Figure 3 provides a clear example of the importance of the interplay between coherent and incoherent dynamics in excitonic energy transport. The biological purpose 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 and 300 K, are shown in figure 4. These plots correspond to the parameter regimes of the FMO 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.

Zoom In Zoom Out Reset image size

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 transformation 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, ($\mathcal {O}(H_{\mathrm {I}}^2)$), 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.

We thank the Leverhulme Trust and the National Research Foundation and Ministry of Education of Singapore for financial support. BWL thanks the Royal Society for a University Research Fellowship. DPSM thanks the EPSRC, CHIST-ERA project SSQN, and CONICET. AN is supported by Imperial College. In addition, we are grateful to Alex Chin, Agata Checinska and Libby Heaney for useful discussions.