Chinese Journal of Chemical Physics  2018, Vol. 31 Issue (4): 421-432

#### The article information

Zhi-hao Gong, Zhou-fei Tang, Jian-shu Cao, Jianlan Wu

Optimal Initialization of a Quantum System for an Efficient Coherent Energy Transfer

Chinese Journal of Chemical Physics, 2018, 31(4): 421-432

http://dx.doi.org/10.1063/1674-0068/31/cjcp1804068

### Article history

Accepted on: June 21, 2018
Optimal Initialization of a Quantum System for an Efficient Coherent Energy Transfer
Zhi-hao Gonga, Zhou-fei Tanga, Jian-shu Caob, Jianlan Wua
Dated: Received on April 14, 2018; Accepted on June 21, 2018
a. Department of Physics, Zhejiang University, Hangzhou 310027, China;
b. Department of Chemistry, Massachusetts Institute of Technology, Cambridge, MA, 02139, USA
*Author to whom correspondence should be addressed. Jianlan Wu, E-mail:jianlanwu@zju.edu.cn
Abstract: For an energy transfer network, the irreversible depletion of excited electron energy occurs through either an efficient flow into an outer energy sink or an inefficient decay. With a small decay rate, the energy transfer efficiency is quantitatively reflected by the average life time of excitation energy before being trapped in the sink where the decay process is omitted. In the weak dissipation regime, the trapping time is analyzed within the exciton population subspace based on the secular Redfield equation. The requirement of the noise-enhanced energy transfer is obtained, where the trapping time follows an exact or approximate 1/Γ-scaling of the dissipation strength Γ. On the opposite side, optimal initial system states are conceptually constructed to suppress the 1/Γ-scaling of the trapping time and maximize the coherent transfer efficiency. Our theory is numerically testified in four models, including a biased two-site system, a symmetric three-site branching system, a homogeneous onedimensional chain, and an 8-chromophore FMO protein complex.
Key words: Noise-enhanced energy transfer    Trapping-free subspace    Optimal initialization    Quantum dissipation
Ⅰ. INTRODUCTION

Optimizing system- and environment-related parameters is a fundamental question in quantum dynamics and thermodynamics, appearing in contexts of energy and charge transfer [1], quantum heat engine [2, 3], optimal control [4, 5], and many other problems. For example, the transfer efficiency of electronic excitation energy in biological photosynthetic protein complexes is maximized ($\sim$100%) at intermediate environmental parameters around the physiological condition. This optimization behavior can be observed in terms of the reorganization energy, the temporal-spatial correlation of bath, the static disorder, and many other parameters [6-22].

The environment surrounding a quantum system induces quantum dissipation [23, 24], which plays an interesting role of adjusting quantum transport processes such as energy transfer. In the strong dissipation limit, the conventional hopping kinetics predicts a small transfer rate and a low efficiency. In the opposite limit of a long-lasting quantum coherence, a delocalized eigenstate (exciton) leads to an instantaneous long-range transfer, but energy oscillates within the quantum system before being irreversibly absorbed by an outer energy sink. The transfer efficiency stays at a low level. As the dissipation strength is gradually applied to disrupt coherence of the quantum system, the transfer efficiency is significantly enhanced until reaching a maximum value at an intermediate dissipation strength. The phenomenon that the efficiency increases with the dissipation strength is termed with different names such as the noise-enhanced energy transfer (NEET) and environment-assisted quantum transport (ENAQT) [6-22].

Taking a biased two-site system as a simple example, we can interpret the NEET using the Förster resonance energy transfer (FRET) theory [25], where the energy transfer rate (more accurately the time integration of the rate kernel) is proportional to the spectral overlap between donor emission and acceptor absorption. The environmental noise broadens two lineshapes, which subsequently increases the spectral overlap and the transfer rate. A similar phenomenon is observed in the classical Kramer's theory, where the friction accelerates the diffusion in the energy space and also leads to the increase of the reaction rate in the weak damping regime [26]. For a general multi-site quantum network, the NEET has been interpreted with the concepts of the invariant subspace [9] and the trapping-free subspace [15]. In the eigen basis representation, the trapping-free subspace consists of excitons free of the irreversible trapping process (orthogonal to the trapping operator) and the assistance of the environmental noise can break this orthogonality for the enhancement of efficiency. In our previous study [15], the trapping-free subspace is demonstrated in a highly-symmetric dendrimer system, and a more comprehensive construction is required.

As a contradiction, the time integration of the rate kernel for an unbiased two-site system approaches the infinity in the complete coherent limit and the transfer efficiency is maximized accordingly. Another example is a homogeneous one-dimensional (1D) chain, which is the simplest polymer model [27]. The coherent transfer efficiency can also be maximized by a ballistic quantum motion. It may seem that such systems disobey the NEET behavior. Instead, we take a different angle to think that the maximized coherent transfer efficiency is induced by a finely tuned initial system state. In other words, the NEET is a universal behavior of a quantum network with an irreversible population depletion but the optimization on the system initialization can suppress the NEET behavior, which is also universal. In this paper, we will perform a mathematical analysis to reveal this conceptual point and verify it in various model systems.

In this work, we briefly review the theoretical modelling of a quantum dynamic network with an irreversible energy trapping process. The trapping time is approximated by its partial value within the exciton population subspace. Following a rational polynomial expression, the partial trapping time is analyzed to derive the constraints of the NEET as well as the optimal initial system states to suppress the $1/\Gamma$-scaling. Four model energy transfer networks are also numerically inspected in detail to verify our theory.

Ⅱ. QUANTUM DISSIPATIVE DYNAMICS WITH AN IRREVERSIBLE TRAPPING PROCESS

For an open quantum system with an irreversible population depletion process, e.g. the trapping process, the time evolution of the system reduced density matrix (RDM) $\rho_{\rm{S}}(t)$ is formally given by [11]

 $\begin{eqnarray} \dot{\rho}_{\rm{S}}(t) = -{\cal L}_{\rm{S}}\rho_{\rm{S}}(t) -{\cal L}_{\rm{t}}[\rho_{\rm{S}}(t)]-{\cal L}_{\rm{d}}\left[\rho_{\rm{S}}(t)\right] \label{eq_01} \end{eqnarray}$ (1)

The three superoperators, $\{{\cal L}_{\rm{S}}, {\cal L}_{\rm{t}}, {\cal L}_{\rm{d}}\}$, refer to the dynamic processes of system quantum oscillation, irreversible trapping (yielding an efficient work), and environment-induced dissipation, respectively. The system Liouville superoperator, ${\cal L}_{{\rm{S}}}$=$i[H_{\rm{S}}, \cdots]$, arises from the commutator of the system Hamiltonian $H_{\rm{S}}$. Notice that an extra unit imaginary number is included in ${\cal L}_{{\rm{S}}}$, while the definition without this number is also widely used in literature [28]. Throughout this paper, the reduced Planck constant is set as unity ($\hbar$=1). In a 'local' basis set $\{|n=1, \cdots, N\rangle\}$, the system Hamiltonian is non-diagonal, expanded as $H^{\rm{L}}_{\rm{S}}$=$\displaystyle\sum_{m, n} H_{mn} |m\rangle\langle n|$. For an energy transfer network with a single excitation, which is the main focus of this paper, $|n\rangle$ represents a quantum state of the excitation localized at site $n$ [29]. Next we build the eigen basis set $\{|\varepsilon_i(i=1, \cdots, N)\rangle\}$ to diagonalize the system Hamiltonian as $H^{\rm{E}}_{\rm{S}}$=$\displaystyle\sum_i \varepsilon_i |\varepsilon_i\rangle\langle \varepsilon_i|$, where $\varepsilon_i$ is the $i$-th eigen energy. In this paper, the two basis representations are distinguished explicitly by the two superscripts, 'L' and 'E'.

For simplicity, a non-Hermitian Hamiltonian $H_{\rm{t}}$ is employed to phenomenologically describe the irreversible trapping process, yielding ${\cal L}_{\rm{t}}[\rho(t)]$=$i (H_{\rm{t}}\rho(t)-\rho(t) H^+_{\rm{t}})$. If a trapping rate $k_{{\rm{t}}; n}$ is defined incoherently at each local site $n$, the trapping Hamiltonian is written as $H^{\rm{L}}_{\rm{t}}$=$\displaystyle\sum_n -i(k_{{\rm{t}}; n}/2)|n\rangle\langle n|$. In the local basis representation, the trapping superoperator ${\cal L}_{\rm{t}}$ is simplified to be [11]

 $\begin{eqnarray} {\cal L}^{\rm{L}}_{{\rm{t}}; mn, m' n'} = \delta_{m', m}\delta_{n', n}\frac{k_{{\rm{t}};m}+k_{{\rm{t}}; n}}{2} \label{eq_02} \end{eqnarray}$ (2)

The trapping rates can be similarly assigned to delocalized excitons (eigenstates) if necessary [19].

The dissipation of a quantum system induced by an interaction between the system and the surrounding environment is reflected by the phenomena of population re-distribution and decoherence [23, 24]. In a microscopic description, we introduce the total Hamiltonian, $H_{\rm{tot}}$=$H_{\rm{S}}$+$H_{\rm{B}}$+$H_{{\rm{S}}{\rm{B}}}$, by excluding the non-Hermitian trapping Hamiltonian $H_\rm{t}$. Here $H_{\rm{B}}$ is the Hamiltonian of the bare bath and $H_{{\rm{S}}{\rm{B}}}$ is the system-bath interaction. In the local basis representation, $H_{{\rm{S}}{\rm{B}}}$ is assumed to follow the form of $H^{\rm{L}}_{{\rm{S}}{\rm{B}}}$=$\displaystyle\sum_n |n\rangle\langle n| B_n$ with $B_n$ being the bath operator. The initial condition of a system-bath factorized state, $\rho_{\rm{tot}}(0)$=$\rho_{\rm{S}}(0)\rho^{\rm{eq}}_{\rm{B}}$ with $\rho^{\rm{eq}}_{\rm{B}}$$\propto$$\exp(-\beta H_{\rm{B}})$, is also assumed.

In the simplest case, the environment is modelled as a classical white noise. The quantum dissipation is described by the Haken-Strobl-Reineker (HSR) model as [30, 31]

 $\begin{eqnarray} {\cal L}^{\rm{L}}_{{\rm{d}}; mn, m' n'} = (1-\delta_{m, n})\delta_{m', m}\delta_{n', n}\Gamma \label{eq_03} \end{eqnarray}$ (3)

where $\Gamma$ is a dephasing rate. In a formal manner, the environment-induced dissipation can be expressed exactly by the Nakajima-Zwanzig projection operator technique [32, 33]. With respect to the reference Hamiltonian, $H_0$=$H_{\rm{S}}$+$H_{\rm{B}}$, the system-bath interaction $H_{{\rm{S}}{\rm{B}}}$ is treated as a perturbation term. In the interaction picture of $H_0$, we perform two unitary transformations, $H_{{\rm{S}}{\rm{B}}}(t)$=$\exp(iH_0t)H_{{\rm{S}}{\rm{B}}}\exp(-iH_0t)$ and $\rho^{({\rm{I}})}_{\rm{S}}(t)$=$\exp(i H_0 t) \rho_{\rm{S}}(t) \exp(-i H_0 t)$. The quantum dissipation is formulated as [24]

 \begin{align} & {{\mathcal{L}}_{\text{d}}}[\rho _{\text{S}}^{(\text{I})}(t)]\to -\int_{0}^{t}{\text{d}}\tau {\mathcal{P}}{{\mathcal{L}}_{\text{SB}}}(t){{\mathcal{T}}_{+}}\left[ {{\text{e}}^{-\int_{\tau }^{t}{\text{d}}{\tau }'{\mathcal{Q}}{{\mathcal{L}}_{\text{SB}}}({\tau }')}} \right]\cdot \\ & \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ {\mathcal{Q}}{{\mathcal{L}}_{\text{SB}}}(\tau ){\mathcal{P}}\rho _{\text{tot}}^{(\text{I})}(\tau ) \\ \end{align} (4)

where ${\cal L}_{{\rm{S}}{\rm{B}}}(t)$=$i[H_{{\rm{S}}{\rm{B}}}(t), \cdots]$ is a Liouville superoperator, and ${\cal P}$=$\rho^{\rm{eq}}_{\rm{B}}\}{\rm{Tr}}_{\rm{B}}\{$ and ${\cal Q}$=${\cal I}$$-$${\cal P}$ are two orthogonal projection operators. The practical difficulty of Eq.(4) is caused by the ambiguous definition of ${\cal Q}$. On the second order perturbation of $H_{{\rm{S}}{\rm{B}}}$, Eq.(4) is simplified to be [24]

 \begin{align} & {{\mathcal{L}}_{\text{d}}}[\rho _{\text{S}}^{(\text{I})}(t)]\approx -\int_{0}^{t}{\text{d}}\tau \text{T}{{\text{r}}_{\text{B}}}\{{{\mathcal{L}}_{\text{SB}}}(t){{\mathcal{L}}_{\text{SB}}}(\tau )\rho _{\text{S}}^{(\text{I})}(\tau )\rho _{\text{B}}^{\text{eq}}\} \\ & \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ =\int_{0}^{t}{\text{d}}\tau \text{T}{{\text{r}}_{\text{B}}}\{[{{H}_{\text{SB}}}(t),[{{H}_{\text{SB}}}(\tau ),\rho _{\text{S}}^{(\text{I})}(\tau )\rho _{\text{B}}^{\text{eq}}]]\} \\ \end{align} (5)

which is a good approximation in the limit of weak dissipation. Furthermore, the Born-Markov approximation and the random phase approximation are applied. In the Schrödinger picture and the eigen basis representation, Eq.(5) is simplified to be

 $\begin{eqnarray} {\cal L}_{\rm{d}}[\rho^{\rm{E}}_{\rm{S}}(t)] &\approx& \sum\limits_{i, j} R_{ii, jj}\rho^{\rm{E}}_{{\rm{S}}; jj}(t)+\sum\limits_{i\neq j} R_{ij, ij} \rho^{\rm{E}}_{{\rm{S}}; ij}(t)\quad \label{eq_06} \end{eqnarray}$ (6)

where $R_{ii, jj}$ and $R_{ij, ij}$ are the Redfield tensors. Eq.(6) is named the secular Redfield equation [29, 34], which can be transformed into the form of the Lindblad equation [24, 25]. Under a condition that the environment is a Gaussian noise (e.g., the bosonic bath), the projection operator ${\cal Q}$ can be explicitly expanded over a series of auxiliary dynamic elements [36-54]. A complete basis set, $\boldsymbol \sigma$=$\{\sigma_0(t)$=$\rho_{\rm{S}}(t), \sigma_1(t), \cdots\}$, is constructed, where $\sigma_h(t)$ is the set of dynamic variables on the $h$-th order of the hierarchic expansion. In a shorthand notation, quantum dissipation is described exactly by the hierarchical equations of motion (HEOM) [36-54],

 $\begin{eqnarray} \dot{\boldsymbol\sigma}(t) = - {\cal W} \boldsymbol\sigma(t) \label{eq_07} \end{eqnarray}$ (7)

where the transition rate matrix ${\cal W}$ is block tridiagonal, ${\cal W}_{h, h'}$=${\cal W}_{h, h} \delta_{h', h}$+${\cal W}_{h, h\pm1}\delta_{h', h\pm 1}$. For the factorized initial state, $\rho_{\rm{tot}}(0)$=$\rho_{\rm{S}}(0)\rho^{\rm{eq}}_{\rm{B}}$, all the high order auxiliary dynamic elements vanish initially, $\sigma_{h\ge 1}(0)$=0, so that the block matrix inversion leads to a time-convolution form,

 $\begin{eqnarray} {\cal L}_{\rm{d}}[\rho_{\rm{S}}(t)] =\int_0^t \textrm{d}\tau {\cal L}_{\rm{d}}(t-\tau)\rho_{\rm{S}}(\tau) \label{eq_08} \end{eqnarray}$ (8)

The dissipation kernel ${\cal L}_{\rm{d}}(t)$ is obtained by the inverse Laplace transform (LT$^{-1}$) of a continued fraction expression [55],

 $\begin{eqnarray} {\cal L}_{\rm{d}}(t) = \mathit{\rm{LT}}^{-1}\bigg[{\cal W}_{0, 1}\cdot\nonumber\\ &&\frac{1}{z+{\cal W}_{1, 1}+{\cal W}_{1, 2}\frac{1}{z+{\cal W}_{2, 2}+\cdots}{\cal W}_{2, 1} }{\cal W}_{1, 0} \bigg] \label{eq_09} \end{eqnarray}$ (9)

In general, the equation of motion in Eq.(1) can be formally solved in the Laplace $z$-space, yielding

 $\begin{eqnarray} \tilde{\rho}_{\rm{S}}(z) = \left[z+{\cal L}_{{\rm{s}}}+{\cal L}_{\rm{t}}+\tilde{{\cal L}}_{\rm{d}}(z)\right]^{-1}\rho_{\rm{S}}(0) \label{eq_10} \end{eqnarray}$ (10)

where $\tilde{\rho}_{\rm{S}}(z)$ and $\tilde{{\cal L}}_{\rm{d}}(z)$ are the Laplace transforms of $\rho_{\rm{S}}(t)$ and ${\cal L}_{\rm{d}}(t)$, respectively. A summary of various quantum dissipation methods can be found in a recent review [1].

For an irreversible quantum dynamic network, a key quantity is the trapping time $\langle t\rangle$, which is the sum of the survival time at each site, i.e., $\langle t\rangle$=$\displaystyle\sum_n \int_0^\infty \rho^{\rm{L}}_{{\rm{S}}; nn}(t) \textrm{d}t$, with $\rho^{\rm{L}}_{{\rm{S}}; nn}(t)$ being the population of site $n$ at time $t$ [11]. Eq.(10) allows us to calculate the trapping time using

 \begin{align} & \langle t\rangle =\text{T}{{\text{r}}_{\text{S}}}\{{{{\tilde{\rho }}}_{\text{S}}}(z=0)\} \\ & \ \ \ \ =\text{T}{{\text{r}}_{\text{S}}}\left\{ {{\left[ {{\mathcal{L}}_{\text{S}}}+{{\mathcal{L}}_{\text{t}}}+{{\widetilde{\mathcal{L}}}_{\text{d}}}(z=0) \right]}^{-1}}{{\rho }_{\text{S}}}(0) \right\} \\ \end{align} (11)

Since $\langle t\rangle$ only depends on the dissipation superoperator at $z$=0, the abbreviation, ${\cal L}_{\rm{d}}$$\equiv$$\tilde{{\cal L}}_{\rm{d}}$($z$=0), is introduced. For an energy transfer network, inefficient energy loss as heat or light can be phenomenologically described by an additional decay process characterized by a decay rate $k_{\mathit{\rm{decay}}}$. The energy transfer efficiency is defined as the cumulated population flow to the energy trapping process, i.e., $q$=$\displaystyle\int_0^\infty {\cal L}_{{\rm{t}}} \rho_{\rm{S}}(t) \textrm{d}t$. If the decay rate $k_{\mathit{\rm{decay}}}$ is much smaller than the trapping rate $k_{\rm{t}}$, the energy transfer efficiency is approximated as [11]

 $\begin{eqnarray} q\approx \frac{1}{1+k_\mathit{\rm{decay}} \langle t\rangle} \label{eq_12} \end{eqnarray}$ (12)

In a real system, $k_{\rm{t}}$ is of the order of ps$^{-1}$ while $k_{\mathit{\rm{decay}}}$ is of the order of ns$^{-1}$. Thus, the trapping time $\langle t\rangle$ in Eq.(11) fully determines the transfer efficiency $q$, i.e., a maximum value of $q$ always corresponds to a minimum value of $\langle t\rangle$.

Ⅲ. TRAPPING TIME IN THE WEAK DISSIPATION LIMIT

In the HSR model [30, 31], the dissipation strength from the classical white noise is purely determined by the dephasing rate $\Gamma$. When the bath is an ensemble of harmonic oscillators, the system-bath interaction is described by the spectral density $J(\omega)$, where $\omega$ is the frequency of a harmonic oscillator. The reorganization energy, $\lambda$=$\displaystyle\int_0^\infty \textrm{d}\omega J(\omega)/\omega$, quantitatively represents the average dissipation strength [23]. At high temperatures, the bosonic bath can be qualitatively mapped onto the HSR model with a linear dependence between $\Gamma$ and $\lambda$ [12, 15]. For conciseness, we use the symbol $\Gamma$ in this section to denote the dissipation strength of a general bath.

A. Noise-enhanced energy transfer (NEET)

In the weak dissipation limit ($\Gamma$$\rightarrow0), the dissipation superoperator follows a Taylor expansion, {\cal L}_{\rm{d}}=\ell_{\rm{d}}\Gamma+O(\Gamma^2), where the reduced \Gamma-independent superoperator is approximated as  \begin{eqnarray} \ell^{\rm{E}}_{{\rm{d}}; ij, kl} \approx \ell^{\rm{E}}_{{\rm{d}}; ii, kk}\delta_{i, j}\delta_{k, l}+\ell^{\rm{E}}_{{\rm{d}}; ij(\neq i), ij}\delta_{k, i}\delta_{l, j} \label{eq_13} \end{eqnarray} (13) in the eigen basis representation. Eq.(13) is obtained from the secular Redfield equation in Eq.(6). In the Liouville space, we partition the RDM into the exciton population and coherence subspaces [56], given by \rho^{\rm{E}}_{\rm{P}}=\displaystyle\sum_i \rho^{\rm{E}}_{{\rm{S}}; ii} |\varepsilon_i\rangle\langle \varepsilon_i| and \rho^{\rm{E}}_{\rm{C}}=\displaystyle\sum_{i\neq j}\rho^{\rm{E}}_{{\rm{S}}; ij}|\varepsilon_i\rangle\langle \varepsilon_j|. The two subscripts, 'P' and 'C', denote the subspaces of population and coherence, respectively. In the eigen basis representation, the three superoperators are written in block matrix forms as  \begin{eqnarray} {\cal L}_{\rm{S}}^{\rm{E}} \hspace{-0.12cm}&=&\hspace{-0.12cm} \left( {\begin{array}{*{20}{c}} 0&0\\ 0&{\cal L}_{\rm{S;C}}^{\rm{E}} \end{array}} \right)\nonumber \\ {\cal L}_{\rm{d}}^{\rm{E}} \hspace{-0.12cm}&=&\hspace{-0.12cm}\left( {\begin{array}{*{20}{c}} \ell_{\rm{d;P}}^{\rm{E}}&0\\ 0&\ell_{\rm{d;C}}^{\rm{E}} \end{array}} \right)\\ {\cal L}_{\rm{t}}^{\rm{E}} \hspace{-0.12cm}&=&\hspace{-0.12cm} \left( {\begin{array}{*{20}{c}} {\cal L}_{\rm{t;P}}^{\rm{E}} &{\cal L}_{\rm{t;PC}}^{\rm{E}}\\ {\cal L}_{\rm{t;CP}}^{\rm{E}}&\ell_{\rm{t;C}}^{\rm{E}} \end{array}} \right)\nonumber \end{eqnarray} (14) In general, the trapping process (\simps) is usually slower than a typical quantum oscillation of excited electrons (\simfs). Through the block matrix inversion of \left[{\cal L}_{\rm{S}}^{\rm{E}} +L_{\rm{d}}^{\rm{E}}+L_{\rm{t}}^{\rm{E}}\right]^{-1} in Eq.(11), we can straightforwardly estimate that the average trapping time \langle t\rangle is reliably approximated by its partial value,  \begin{eqnarray} \langle t\rangle^{\rm{E}}_{\rm{P}} = {\rm{Tr}}_{\rm{S}}\left\{\left[\Gamma\ell^{\rm{E}}_{{\rm{d}}; {\rm{P}}} + {\cal L}^{\rm{E}}_{{\rm{t}}; {\rm{P}}} \right]^{-1} \rho^{\rm{E}}_{\rm{P}}(0) \right\} \label{eq_15} \end{eqnarray} (15) where the matrices in the trace {\rm{Tr}}_{\rm{S}} are restricted to the exciton population subspace. Despite the fact that only dissipation and trapping appear in Eq.(15), the system Hamiltonian enters implicitly through the unitary transformation from the local to eigen basis set. Notice that the NEET may appear in a fast trapping process due to a strong influence of hopping kinetics, which is however beyond the scope of this paper. For a finite N-site energy transfer network, \ell^{\rm{E}}_{{\rm{d}}; {\rm{P}}} and {\cal L}^{\rm{E}}_{{\rm{t}}; {\rm{P}}} are two N$$\times$$N matrices. Through a standard matrix inversion for [\Gamma\ell^{\rm{E}}_{{\rm{d}}; {\rm{P}}}+{\cal L}^{\rm{E}}_{{\rm{t}}; {\rm{P}}}]^{-1} [57], the partial trapping time \langle t\rangle^{\rm{E}}_{\rm{P}} in Eq.(15) can be expressed as a rational polynomial,  \begin{eqnarray} \langle t \rangle^{\rm{E}}_{\rm{P}} &=& \frac{\displaystyle\sum\limits_{k=0}^{N-1}a_k \Gamma^{k}}{\displaystyle\sum\limits_{k=0}^{N-1}b_k \Gamma^{k}}= t_0+\sum\limits_{k=1}^{N-1}\frac{t_k}{\Gamma+\Gamma_k} \label{eq_16} \end{eqnarray} (16) where all the four parameter sets, \{a_k, b_k, t_k, \Gamma_k (k=0, \cdots, N$$-$1)$\}$ are functions of $H_{\rm{S}}$, $\ell^E_{{\rm{d}}; {\rm{P}}}$, and ${\cal L}^{\rm{E}}_{{\rm{t}}; {\rm{P}}}$. The two sets of $\{a_k, t_k\}$ vary with the initial exciton population $\rho^{\rm{E}}_{\rm{P}}(0)$, while the other two sets of $\{b_k, \Gamma_k\}$ are independent of $\rho^{\rm{E}}_{\rm{P}}(0)$. The constraint, $b_N$$\propto$$\mathit{\rm{Det}}[\ell^{\rm{E}}_{{\rm{d}}; {\rm{P}}}]$=0, always holds due to the fact that quantum dissipation conserves the total population. Another constraint, $b_0$$\propto$$\mathit{\rm{Det}}[{\cal L}^{\rm{E}}_{{\rm{t}}; {\rm{P}}}]$, can be also obtained. Here the symbol $\mathit{\rm{Det}}[\cdots]$ denotes the matrix determinant. The set of $\{\Gamma_k\}$ is composed of the roots of

 $\begin{eqnarray} \sum\limits_{k=0}^{N-1}b_k \Gamma^{k}=0 \label{eq_17} \end{eqnarray}$ (17)

To include a weak contribution of the exciton coherence, we add a linear $\Gamma$-term for the total trapping time, giving

 $\begin{eqnarray} \langle t\rangle\approx \langle t\rangle^{\rm{E}}_{\rm{P}}+\delta t_0+f_{\mathit{\rm{hop}}}\Gamma \label{eq_18} \end{eqnarray}$ (18)

where the two positive parameters, $\delta t_0$ and $f_{\mathit{\rm{hop}}}$, depend on the three superoperators and the initial system state. The correction term $\delta t_0$ is usually much smaller than $t_0$ and can be neglected. This linear $\Gamma$-term is the reminiscent of incoherent hopping kinetics in the weak dissipation limit [15]. The total trapping time in the weak dissipation limit is expanded as

 $\begin{eqnarray} \langle t \rangle &\approx& \left(t_0+\delta t_0+\sum\limits_{k=1}^{N-1} \frac{t_k}{\Gamma_k}\right) +\nonumber\\ &&\left(f_{\mathit{\rm{hop}}}-\sum\limits_{k=1}^{N-1}\frac{t_k}{\Gamma^2_k}\right) \Gamma +O(\Gamma^2) \label{eq_19} \end{eqnarray}$ (19)

As a result, the condition

 $\begin{eqnarray} \sum\limits_{k=1}^{N-1}t_k/\Gamma^2_k>f_{\mathit{\rm{hop}}} \label{eq_20} \end{eqnarray}$ (20)

is a general requirement of the NEET that $\langle t\rangle$ decreases with $\Gamma$ (i.e. efficiency $q$ increases with $\Gamma$) in the weak dissipation regime. Under the opposite condition, $\langle t\rangle$ increases with $\Gamma$ and the efficiency is maximized by the coherent energy transfer without dissipation.

However, one or more zero roots in Eq.(17) dramatically change the $\Gamma$-dependence of the trapping time. For example, the partial trapping time $\langle t\rangle^{\rm{E}}_{\rm{P}}$ with $\Gamma_1$=0 is changed to

 $\begin{eqnarray} \langle t \rangle^{\rm{E}}_{\rm{P}} = \left(t_0+\sum\limits_{k=2}^{N-1}\frac{t_k}{\Gamma_k}\right)+ \frac{t_1}{\Gamma}+O(\Gamma) \label{eq_21} \end{eqnarray}$ (21)

where $t_1$ is nonnegative for a physical trapping time. The condition necessary for the exact $1/\Gamma$-scaling in Eq.(21) is

 $\begin{eqnarray} b_0 \propto \mathit{\rm{Det}}\left[{\cal L}^{\rm{E}}_{{\rm{t}}; {\rm{P}}}\right] = 0 \label{eq_22} \end{eqnarray}$ (22)

which is our definition of the rigorous trapping-free subspace $\Phi_\perp$ [15]. For the incoherent trapping process defined in Eq.(2), ${\cal L}^{\rm{E}}_{{\rm{t}}; {\rm{P}}}$ is a diagonal matrix satisfying ${\cal L}^{\rm{E}}_{{\rm{t}}; ii, jj}$=${\cal L}^{\rm{E}}_{{\rm{t}}; ii}\delta_{i, j}$. The matrix determinant becomes $\mathit{\rm{Det}}\left[{\cal L}^{\rm{E}}_{{\rm{t}}; {\rm{P}}} \right]$=$\Pi_{i=1}^N {\cal L}^{\rm{E}}_{{\rm{t}}; ii}$ with ${\cal L}^{\rm{E}}_{{\rm{t}}; ii}$=$\displaystyle\sum_{n=1}^N \langle \varepsilon_i|n\rangle^2 k_{{\rm{t}}; n}$. The trapping-free subspace $\Phi_\perp$ composed of $N_{\rm{t}}$ ($<$$N) trapping-free excitons ({\cal L}^{\rm{E}}_{{\rm{t}}; ii}=0) requires the following two conditions to be satisfied simultaneously: (ⅰ) at least one site is not connected to the energy sink, giving k_{{\rm{t}}; n}=0; (ⅱ) each exciton state in \Phi_\perp is a linear combination of these trapping-free sites, giving \langle \varepsilon_i|n'\rangle=0 for |\varepsilon_i\rangle\in\Phi_\perp and k_{{\rm{t}}; n'}$$\neq$0.

The rigorous solution of $b_0$=0 implies a high symmetry in the system Hamiltonian, which is often not the case in a general quantum network. Instead, we consider a weaker constraint: the characteristic dissipation strengths $\{\Gamma_{k}\}$ are well separated by their magnitudes, satisfying 0$\leftarrow$$|\Gamma_1|, \cdots, |\Gamma_{N_{\rm{t}}}|$$\ll$$|\Gamma_{N_{\rm{t}}+1}|, \cdots, |\Gamma_N|. One possible solution for this constraint is  \begin{eqnarray} |b_0|\ll |b_{N-1}| \Gamma_\textrm{c}^{N-1} \label{eq_23} \end{eqnarray} (23) where \Gamma_\textrm{c} is a critical dissipation strength to balance the units of b_0 and b_{N-1}. Roughly speaking, \Gamma_\textrm{c} can be considered as the boundary of weak dissipation. The exciton states with the small trapping rates ({\cal L}^{\rm{E}}_{{\rm{t}}; ii}$$\rightarrow$0) form an approximate trapping-free subspace $\Phi_\perp$. In the weak dissipation regime, the partial trapping time $\langle t\rangle^{\rm{E}}_{\rm{P}}$ becomes

 $\begin{eqnarray} \begin{array}{l} \left\langle t \right\rangle _{\rm{P}}^{\rm{E}} \approx \left[{{t_0} + \displaystyle\sum\limits_{k = {N_{\rm{t}}} + 1}^{N-1} {\frac{{{t_k}}}{{{\Gamma _k}}}} } \right] + \displaystyle\sum\limits_{k = 1}^{{N_{\rm{t}}}} {\frac{{{t_k}}}{{\Gamma + {\Gamma _k}}}} \\ \hspace{3cm}\textrm{for}\ \Gamma\ll|\Gamma_{k(>N_{\rm{t}})}|;\\ \left\langle t \right\rangle _{\rm{P}}^{\rm{E}} \approx \left[{{t_0} + \displaystyle\sum\limits_{k = {N_{\rm{t}}} + 1}^{N-1} {\frac{{{t_k}}}{{{\Gamma _k}}}} } \right] + \displaystyle\frac{{\displaystyle\sum\limits_{k = 1}^{{N_{\rm{t}}}} {{t_k}} }}{\Gamma }\\ \hspace{3cm}\textrm{for}\ |\Gamma_{k(N_{\rm{t}})}| \end{array} \label{eq_24} \end{eqnarray}$ (24)

To achieve the approximate $1/\Gamma$-scaling in Eq.(24), the two sets of coefficients, $\{t_{k(\le N_{\rm{t}})}\}$ and $\{t_{k(>N_{\rm{t}})}\}$, are needed to be comparable in magnitude.

B. The optimal initial system state for the disappearance of the $1/\Gamma$-scaling

In Section Ⅲ.A, we formulate the exact and approximate $1/\Gamma$-scalings of the trapping time in Eqs. (21) and (24), respectively. On the opposite side, the $1/\Gamma$-scaling can disappear by an appropriate initial system state, $\rho^{\rm{E}}_{\rm{P}}(0)$. For a rigorous trapping-free subspace $\Phi_\perp$, a straightforward way is to prepare zero population in $\Phi_\perp$, i.e., $\rho^{\rm{E}}_{ii}(0)$=0 for $|\varepsilon_i\rangle$$\in$$\Phi_\perp$ [15]. For a general network, we apply a stronger requirement,

 $\begin{eqnarray} a_0:a_1:\cdots:a_{N-1} = b_0:b_1:\cdots:b_{N-1} \label{eq_25} \end{eqnarray}$ (25)

from which all the $\Gamma$-dependent terms vanish in Eq.(16). The partial trapping time $\langle t\rangle^{\rm{E}}_{\rm{P}}$ becomes a constant and the total trapping time $\langle t\rangle$ increases linearly with the dissipation strength $\Gamma$, i.e., $\langle t\rangle$$\approx$$(t_0$+$\delta t_0$)+$f_{\mathit{\rm{hop}}} \Gamma$+$O(\Gamma^2)$. The energy transfer efficiency is globally or locally maximized in the coherent regime (\Gamma\rightarrow0). Next we provide a mathematical procedure of optimizing the initial system state. In the subspace of exciton population, the partial reduced dissipation superoperator \ell^{\rm{E}}_{{\rm{d}}; {\rm{P}}} is diagonalized as follows,  \begin{eqnarray} {\cal D}={\cal V}^{-1}\ell^{\rm{E}}_{{\rm{d}}; {\rm{P}}}{\cal V} = \left( {\begin{array}{*{20}{c}} 0&{}&{}\\ {}&{{\Lambda _2}}&{}\\ {}&{}& \ddots \end{array}} \right) \label{eq_26} \end{eqnarray} (26) where \Lambda_i are the eigen dissipation rates reduced by the dissipation strength \Gamma. The zero eigen rate, \Lambda_1=0, appears due to the condition that system population is conserved by dissipation. The matrix \cal V is not necessarily unitary since \ell^{\rm{E}}_{{\rm{d}}; {\rm{P}}} can be non-Hermitian. The same transformation is applied to the trapping superoperator and the initial RDM, leading to {\cal T}={\cal V}^{-1}{\cal L}^{\rm{E}}_{{\rm{t}}; {\rm{P}}}{\cal V} and \varrho(0)={\cal V}^{-1}\rho^{\rm{E}}_{\rm{P}}(0). The partial trapping time \langle t\rangle^{\rm{E}}_{\rm{P}} in Eq.(15) is rewritten as  \begin{eqnarray} \langle t\rangle^{\rm{E}}_{\rm{P}} = \mathit{\rm{Tr}}_{\rm{S}}\left\{\mathcal V(\Gamma{\cal D}+{\cal T})^{-1}\varrho(0)\right\} \label{eq_27} \end{eqnarray} (27) After tedious but straightforward steps, we obtain b_0={\rm{Det}}[{\cal T}] and b_{N-1}=\displaystyle\prod_{i=2}^{N} \Lambda_{i} [{\cal T}]_{11}. Eq.(23) for the approximate 1/\Gamma-scaling of the NEET becomes  \begin{eqnarray} |\mathit{\rm{Det}}\left[{\cal T}\right]| \ll \prod\limits_{i=2}^{N} |\Lambda_{i} [{\cal T}]_{11}| \Gamma_c^{N-1} \label{eq_28} \end{eqnarray} (28) which implies an upper limit of the trapping rates. With respect to the zero eigen rate of \Lambda_1=0, the optimal condition in Eq.(25) to maximize the coherent energy transfer efficiency is explicitly given by  \begin{align} & {{[\varrho (0)]}_{1}}:{{[\varrho (0)]}_{2}}:\cdots :{{[\varrho (0)]}_{N}}= \\ & \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ {{[T]}_{11}}:{{[T]}_{21}}:\cdots :{{[T]}_{N1}} \\ \end{align} (29) Since the initial exciton coherence varies freely, an infinite possibility of the initial system states (either pure or mixed) can satisfy Eq.(29). If the reduced dissipation superoperator \ell^{\rm{E}}_{{\rm{d}}; {\rm{P}}} is Hermitian as in the HSR model, the eigenvector associated with \Lambda_1=0 is the identity operator, {\cal I}=\displaystyle\sum_{i=1}^N |\varepsilon_i\rangle\langle \varepsilon_i|. The trapping time is equivalent to the survival time along this special eigenvector. The matrix \cal T can be block diagonal as required by the symmetry of the system Hamiltonian H_{\rm{S}}. The calculation of the trapping time is reduced inside an M (\leN)-dimensional subspace including the eigenvector $\cal I$. Consequently, Eqs. (28) and (29) are modified by replacing $\{N, \varrho(0), {\cal T}\}$ with $\{M, \varrho_M(0), {\cal T}_M\}$, where $\varrho_M(0)$ and ${\cal T}_M$ are defined in the $M$-dimensional subspace. The optimal requirement of $\rho_{\rm{S}}(0)$ to maximize the energy transfer efficiency becomes more flexible.

Ⅳ. NEET AND OPTIMAL SYSTEM INITIALIZATION IN MODEL SYSTEMS

In Section Ⅲ, we have derived the requirements of the rigorous and approximate $1/\Gamma$-scaling for the NEET. In the same framework, the initial system state to suppress the $1/\Gamma$-scaling is theoretically obtained, which optimizes the coherent energy transfer efficiency in the weak dissipation regime. In this section, we verify our theory in four model systems of energy transfer, as shown in FIG. 1.

 FIG. 1 The four energy transfer networks studied in this work. (a) a biased two-site system, (b) a symmetric three-site branching system, (c) a homogeneous 1D chain, and (d) an 8-chromophore FMO monomer. The trapping site of each system is highlighted in red color
A. Biased two-site system

The first example is a biased two-site system (see FIG. 1(a)). The system Hamiltonian is defined in the local basis as $H^{\rm{L}}_{\rm{S}}$=$\Delta |1\rangle\langle 1|$+$J(|1\rangle\langle 2|$+$|2\rangle\langle 1|)$, where $\Delta$ is the site energy detuning and $J$ is the site-site coupling. The HSR model in Eq.(3) is applied to model the quantum dissipation. An incoherent trapping process is assumed with the trapping rate $k_{\rm{t}}$ at site 2.

The analytical expression of the trapping time for this model was shown previously [11]. Here we apply the theoretical procedure developed in Section Ⅲ to analyze the NEET and determine the optimal initial system states. In the exciton population subspace, the two relevant superoperators are explicitly written as

 $\begin{eqnarray} \begin{array}{l} \ell^{\rm{E}}_{{\rm{d}}; {\rm{P}}} = \displaystyle\frac{\sin^22\theta }{2}\left( {\begin{array}{*{20}{c}} 1&{{\rm{ - }}1}\\ {{\rm{ - }}1}&1 \end{array}} \right)\\ {\cal L}^{\rm{E}}_{{\rm{t}}; {\rm{P}}} = k_{\rm{t}}\left( {\begin{array}{*{20}{c}} {\sin {\theta ^2}}&0\\ \theta &{\cos {\theta ^2}} \end{array}} \right)\\ \end{array} \label{eq_30} \end{eqnarray}$ (30)

with $\theta$=$-\arctan(2J/\Delta)/2$. Following Eq.(28), we determine a condition of $k_{\rm{t}}$$\ll$$2 \Gamma_\textrm{c}$ for the NEET. For a weak site-site coupling, an intermediate dephasing rate around the energy detuning, $\Gamma_\textrm{c}$$\approx$$|\Delta|$, separates the weak and strong dissipation regimes. To achieve an apparent $1/\Gamma$-scaling, we propose a scenario of a slow trapping process $k_{\rm{t}}$$\ll$$|\Delta|$ accompanied with a large energy mismatch $|J|$$\ll$$|\Delta|$, which is consistent with the result in Ref.[11]. The approximate trapping-free subspace $\Phi_\perp$ is composed of a single exciton state, $|\varepsilon_1\rangle$=$\cos\theta|1\rangle$+$\sin\theta|2\rangle$ with $\theta$$\rightarrow0. The total trapping time is explicitly given by  \begin{eqnarray} \langle t\rangle = t_0+ \frac{t_1}{\Gamma+k_{\rm{t}}/2} +\delta t_0 + f_{\mathit{\rm{hop}}} \Gamma \label{eq_31} \end{eqnarray} (31) where the parameters, \{t_0, \delta t_0, t_1, f_{\mathit{\rm{hop}}}\}, depend on the initial RDM. The first two terms on the right hand side of Eq.(31) arise from the partial trapping time \langle t\rangle^{\rm{E}}_{\rm{P}}. As shown in FIG. 2(a), the partial trapping time \langle t\rangle^{\rm{E}}_{\rm{P}} excellently describes the total trapping time \langle t\rangle in the weak dissipation regime for an example inital RDM, \rho^{\rm{L}}_{{\rm{S}}}(0)=|1\rangle\langle 1|. With the increase of energy mismatch, the approximate 1/\Gamma-scaling of \langle t\rangle^{\rm{E}}_{\rm{P}} gradually becomes exact.  FIG. 2 The trapping time vs. the dephasing rate for a biased two-site sytem (J=1 and k_{\rm{t}}=0.1). (a) The results of the initial system state at \rho^{\rm{L}}_{\rm{S}}(0)=|1\rangle\langle 1|. Subtracted by t_0=2/k_{\rm{t}}, the circles are the total trapping time \langle t\rangle defined in Eq.(11); the solid lines are the partial trapping time \langle t\rangle^{\rm{E}}_{\rm{P}} defined in Eq.(15); the dashed lines are the 1/\Gamma-scaling of (\Delta^2/2J^2)\Gamma^{-1}. The data in black, red, and blue colors refer to \Delta=2, 5, 10, respectively. (b) With \Delta=2, the solid and dashed lines are the results of the initial system states at \phi^{\rm{L}}_{{\rm{S}}}(0)=(|1\rangle$$-$$|2\rangle)/\sqrt{2} and \rho^{\rm{L}}_{{\rm{S}}}(0)=(|1\rangle\langle1|+|2\rangle\langle2|)/2, respectively To suppress the 1/\Gamma-scaling in the trapping time and achieve an optimized coherent energy transfer, we solve Eq.(29) and determine an requirement,  \begin{eqnarray} \rho^{\rm{L}}_{{\rm{S}}; 11}(0)+\frac{2J}{\Delta}\mathit{\rm{Re}} \{\rho^{\rm{L}}_{{\rm{S}}; 12}(0)\}=0 \label{eq_32} \end{eqnarray} (32) in the local basis representation. A trivial solution of Eq.(32) is \rho^{\rm{L}}_{\rm{S}}(0)=|2\rangle\langle 2|, which means that the system is prepared initially at the trapping site 2. Instead, an infinite number of nontrivial optimal initial system states exist, satisfying Eq.(32). One example is a quantum pure state, \phi^{\rm{L}}_{\rm{S}}(0)$$\propto$$2J|1\rangle$$-$$\Delta|2\rangle, and the calculation of \langle t\rangle for this initial RDM is shown in FIG. 2(b). In the coherent limit (\Gamma$$\rightarrow$0), the trapping time from $\phi^{\rm{L}}_{{\rm{S}}}(0)$ is very close to the global minimum result from $\rho^{\rm{L}}(0)$=$|2\rangle\langle 2|$. As a comparison, we construct a quantum mixed state, $\rho^{\rm{L}}_{\rm{S}}(0)$$\propto$$4|J|^2|1\rangle\langle 1|$+$|\Delta|^2|2\rangle\langle 2|$. By violating the requirement of the site-site coherence in Eq.(32), the coherent trapping time $\langle t\rangle_{\Gamma=0}$ is increased significantly and $\langle t\rangle$ is minimized at an intermediate $\Gamma$ (see FIG. 2(b)).

B. Symmetric three-site branching system

The second example is a symmetric three-site branching system (see FIG. 1(b)). The system Hamiltonian is given by $H^{\rm{L}}_{\rm{S}}$=$J(|1\rangle\langle 2|$+$|2\rangle\langle 1|$+$|2\rangle\langle 3|$+$|3\rangle\langle 2|)$, in the local basis representation. The dissipation is approximated by the HSR model, while the trapping rate $k_{\rm{t}}$ is assigned to the middle site 2.

The explicit expression of the trapping time was also provided previously [11]. Following the procedure in Section Ⅲ, we obtain dissipation and trapping superoperators,

 $\begin{eqnarray} \begin{array}{l} \ell^{\rm{E}}_{{\rm{d}}; {\rm{P}}}= \displaystyle\frac{1}{8} \left( {\begin{array}{*{20}{c}} 4&{ - 2}&{ - 2}\\ { - 2}&5&{ - 3}\\ { - 2}&{ - 3}&5 \end{array}} \right)\\ [2ex] {\cal L}^{\rm{E}}_{{\rm{t}}; {\rm{P}}} = \displaystyle\frac{k_{\rm{t}}}{2}\left( {\begin{array}{*{20}{c}} 0&0&0\\ 0&1&0\\ 0&0&1 \end{array}} \right)\\ \end{array} \label{eq_33} \end{eqnarray}$ (33)

in the exciton population subspace. The first exciton state, $|\varepsilon_1\rangle$=$(|1\rangle -|3\rangle)/\sqrt{2}$, is orthogonal to the trapping site 2 and composes a rigorous one-element trapping-free subspace $\Phi_\perp$ with ${\cal L}^{\rm{E}}_{{\rm{t}}; 11}$=0 [15]. FIG. 3(a) demonstrates that the total trapping time $\langle t\rangle$ follows the exact $1/\Gamma$-scaling for a non-optimal initial system state, $\rho^{\rm{L}}_{{\rm{S}}}(0)$=$(|1\rangle\langle 1|$+$|3\rangle\langle 3|)/2$.

 FIG. 3 The trapping time vs. the dephasing rate for a symmetric three-site model ($J$=1). (a) The results of the initial system state at $\rho^{\rm{L}}_{\rm{S}}(0)$=$(|1\rangle\langle 1|$+$|3\rangle\langle 3|)/2$. The circles and solid lines are the total trapping time $\langle t\rangle$ and the partial trapping time $\langle t\rangle^{\rm{E}}_{\rm{P}}$ (i.e., the exact $1/\Gamma$-scaling), respectively. The data in black and red colors refer to $k_{\rm{t}}$=0.1 and 1.0, respectively. (b) The results of an optimal initial system state at $\phi^{\rm{L}}_{\rm{S}}(0)$=$(|1\rangle$+$|3\rangle)/\sqrt{2}$. The solid and dashed lines refer to $k_{\rm{t}}$=0.1 and 1.0, respectively

As discussed in Section Ⅲ.B, a straightforward way to optimize the coherent energy transfer is to avoid the initial population at the trapping-free subspace, i.e., $\rho^{\rm{E}}_{{\rm{S}}; 11}(0)$=0. Alternatively, we apply the mathematical procedure in Section Ⅲ.B to solve the optimization requirement in Eq.(29), which leads to

 $\begin{eqnarray} \rho^{\rm{L}}_{11}(0)+\rho^{\rm{L}}_{33}(0) = 2 \mathit{\rm{Re}}\{\rho^{\rm{L}}_{13}(0) \} \label{eq_34} \end{eqnarray}$ (34)

in the local basis representation. Eq.(34) is identical to the condition of zero population in the trapping-free subspace. Here we design an optimal initial system state, $\phi^{\rm{L}}_{\rm{S}}(0)$=$(|1\rangle+|3\rangle)/\sqrt{2}$, which shares the same site population as in the above non-optimal state but differs in the site-site coherence. The trapping time is changed to a linearly increasing function of $\Gamma$ (see FIG. 3(b)), which indicates the relevance of an appropriate site-site coherence in quantum optimization. The underlying symmetry argument shows that our analysis is not limited to the simple HSR model, but valid in general environments, as demonstrated in our previous study of a tree-like dendrimer system [15].

C. Homogeneous 1D chain systems

The third example is a homogeneous 1D $N$-site chain (see FIG. 1(c)). In the local basis representation, the system Hamiltonian is written as

 $\begin{eqnarray} H^{\rm{L}}_{\rm{S}} = \sum\limits_{n=1}^{N-1} J (|n\rangle\langle n+1|+|n+1\rangle\langle n|) \label{eq_35} \end{eqnarray}$ (35)

with the nearest neighboring interaction. The dissipation is simulated by the HSR model, while the trapping process is defined by an irreversible rate $k_{\rm{t}}$ at the end site $N$.

Following the procedure in Section Ⅲ.B, we numerically calculate Eq.(28) and determine an upper trapping rate limit, $k_{\rm{t}}$$\ll$$N \Gamma_\textrm{c}(N)$, where the critical dephasing rate roughly follows $\Gamma_\textrm{c}(N)$$\sim$$N^{-1}J$. Together, we reach a slow trapping rate of $k_{\rm{t}}$$\ll$$J$, which can lead to the NEET. The participation coefficient of the trapping site in each exciton state decreases with the increased chain size $N$ so that the trapping rate of each exicton decreases accordingly, i.e. ${\cal L}^{\rm{E}}_{{\rm{t}}; ii}$($i$=1, $\cdots, N)$$\rightarrow0 for N$$\rightarrow$$\infty. The NEET is thus expected to be more pronounced with a long spatial extension (N$$\gg$1).

To extract the approximate $1/\Gamma$-scaling, we prepare the initial population at the third site, $\rho^{\rm{L}}_{\rm{S}}(0)$=$|3\rangle\langle 3|$. The usual condition, $\rho^{\rm{L}}_{\rm{S}}(0)$=$|1\rangle\langle 1|$, is actually an optimal initial state and will be discussed later in this subsection. With a small trapping rate, $k_{\rm{t}}$=$0.1J$, we numerically calculate the total trapping time $\langle t\rangle$ and its partial value $\langle t\rangle^{\rm{E}}_{\rm{P}}$ for the chain sizes of 5$\le$$N$$\le$40. FIG.~4(a) shows that $\langle t\rangle^{\rm{E}}_{\rm{P}}$ agrees excellently with $\langle t\rangle$ in the weak dissipation limit. Furthermore, the partial trapping time $\langle t\rangle^{\rm{E}}_{\rm{P}}$ can be fitted by

 $\begin{eqnarray} \langle t\rangle^{\rm{E}}_{\rm{P}} \approx \frac{N}{k_\textrm{t}}+\frac{t_1}{\Gamma+\Gamma_1}+\frac{t_2}{\Gamma+\Gamma_2} \label{eq_36} \end{eqnarray}$ (36)

over a broad range of $\Gamma$. As shown in FIG. 4(a), the difference between the exact and fitting results of $\langle t\rangle^{\rm{E}}_{\rm{P}}$ cannot be seen with naked eyes. The two characteristic dephasing rates, $\Gamma_1$ and $\Gamma_2$, monotonically decrease with the chain size $N$, as shown in FIG. 4(b). In the condition of $\max(\Gamma_1, \Gamma_2)$$\ll$$\Gamma$$\ll$$|\Gamma_{k(>2)}|$, we reach the exact $1/\Gamma$-scaling,

 $\begin{eqnarray} \langle t\rangle^{\rm{E}}_{\rm{P}} \approx \frac{N}{k_{\rm{t}}}+\frac{t_1+t_2}{\Gamma} \label{eq_37} \end{eqnarray}$ (37)
 FIG. 4 (a) The trapping time subtracted by $N/k_{\rm{t}}$ vs. the dephasing rate for the 1D $N$-site chain. The initial system state is at $\rho^{\rm{L}}_{\rm{S}}(0)$=$|3\rangle\langle 3|$. The circles and solid lines are the calculated results of $\langle t\rangle$ and $\langle t\rangle^{\rm{E}}_{\rm{P}}$. The dashed lines are the fitting results of Eq.(36). The dotted-dashed line is the exact $1/\Gamma$-scaling for $N$=40. The data in black, red, blue and green colors refer to $N$=5, 10, 20, and 40, respectively. (b) The fitting parameters of $\Gamma_1$ and $\Gamma_2$ in Eq.(36) vs. the chain size $N$. The solid lines with circles and diamonds are the results of $\Gamma_1$ and $\Gamma_2$, respectively. Here the parameters are $J$=1 and $k_{\rm{t}}$=0.1

which is confirmed for the result of $N$=40 in FIG. 4(a).

Next we calculate the optimal initial system state of the coherent energy transfer. Without tedious details, Eq.(29) is transformed into an $M$-equation array ($k$=1, 2, $\cdots$, $M$),

 $\begin{eqnarray} \sum\limits_{i=2}^{N-1}x^k_{i}\rho^{\rm{L}}_{i}(0)+\sum\limits_{i=1}^{N-2}\sum\limits_{j=1}^{(N-i)/2} 2 x^k_{i, i+2j}\mathit{\rm{Re}}\rho^{\rm{L}}_{i, i+2j}(0)=0\quad \label{eq_38} \end{eqnarray}$ (38)

with $M$=$(N$$-1)/2 for an odd N and M=N/2$$-$1 for an even $N$. For conciseness, we will not provide the explicit forms of the coefficients, $\{x^k_{i}, x^k_{i, i+2j}\}$. However, there are several important properties to be noticed: $\{x^k_{i}, x^k_{i, i+2j}\}$ are independent of the site-site coupling $J$ and the trapping rate $k_{\rm{t}}$; the coefficients obey a mirror symmetry, $x^k_{i}$=$x^k_{N+1-i}$ and $x^k_{i, i+2j}$=$x^k_{N+1-i-2j, N+1-i}$, together with another constraint, $x^k_{i}$=$x^k_{i-1, i+1}$. Interestingly, a unique initial state satisfying Eq.(38) for an arbitrary chain size $N$ is $\rho^{\rm{L}}_{\rm{S}}(0)$=$|1\rangle\langle1|$, where all the population is initially prepared at the starting site. This preference of the coherent energy transfer is consistent with the ballistic quantum diffusion in the infinite homogeneous 1D chain [58, 59]. On the other hand, Eq.(38) allows an infinite number of $N$-dependent solutions (quantum pure and mixed states). Here we present an example,

 $\begin{eqnarray} \rho^{\rm{L}}_{\rm{S}}(0) &=& \frac{1}{2(N+1)}\bigg[3|1\rangle\langle 1|+3|N\rangle\langle N|+2\sum\limits_{i=2}^{N-1}|i\rangle\langle i| -\nonumber\\ &&\sum\limits_{i=1}^{N-2}\left(|i\rangle\langle i+2|+|i+2\rangle\langle i| \right) \bigg] \label{eq_39} \end{eqnarray}$ (39)

As shown in FIG. 5, the total trapping time $\langle t\rangle$ follows the linear $\Gamma$-function exactly for both optimal initial system states.

 FIG. 5 The trapping time vs. the dephasing rate for the 1D $N$-site chain, with the two optimal initial system states satisfying Eq.(38). The solid and dashed lines are the results of $\rho^{\rm{L}}_{\rm{S}}(0)$=$|1\rangle\langle 1|$ and the initial state defined in Eq.(39). The data in black and blue colors refer to $N$=20 and 40, respectively. Here the parameters are $J$=1 and $k_{\rm{t}}$=0.1

By comparing the results of the non-optimal and optimal initial system states in FIGs. 4 and 5, we observe that the trapping time of $\rho^{\rm{L}}_{\rm{S}}(0)$=$|3\rangle\langle 3|$ in the 20-site chain is larger than that of $\rho^{\rm{L}}_{\rm{S}}(0)$=$|1\rangle\langle 1|$ in the 40-site chain for $\Gamma$$<$$10^{-3}$. In contradiction to a classical diffusion picture, the energy can be transferred faster in the quantum coherent limit even when the distance is doubled. The basic mechanism of this phenomenon lies on the the mirror symmetry in the initial system state compatible with dissipation and trapping.

D. Eight-site FMO model

Our final example is the 8-chromophore Fenna-Matthews-Olson (FMO) protein complex (see FIG. 1(d)), which is an important light-harvesting system in green sulfur bacteria [60, 61]. The effective Hamiltonian of an FMO monomer is taken from Refs. [13, 62]. The influence of the bosonic bath is simulated by a Debye spectral density,

 $\begin{eqnarray} J(\omega) = \frac{2\lambda}{\pi} \frac{\omega\omega_\textrm{D}}{\omega^2+\omega^2_\textrm{D}} \label{eq_40} \end{eqnarray}$ (40)

where $\lambda$ is the reorganization energy and $\omega_\textrm{D}$ is the Debye frequency. The trap process is assigned to bacteriochlorophyll (BChl) 3 with a typical trapping rate $k_{\rm{t}}$=1 ps$^{-1}$. To be consistent with the physiological condition, we set the Debye frequency at $\omega^{-1}_\textrm{D}$=50 fs and temperature at $T$=300 K.

For the non-symmetric FMO system, the upper limit of the trapping rate for the NEET is estimated by Eq.(28) as $k_{\rm{t}}$$\ll$$4000~\mathit{\rm{ps}}^{-1}$, where the critical dissipation strength is approximated as $\lambda_\textrm{c}$$\approx50 cm^{-1}. As the realistic trapping rate (\sim1 ps^{-1}) is much smaller than this upper limit, we expect a strong NEET behavior. The intrinsic mechanism is that the excitons (eigen states) in the FMO system are highly localized due to both the large energy mismatch and the long spatial extension. The eigen states almost orthogonal to the trapping site, BChl 3, compose the approximate trapping-free subspace. As a demonstration, we consider a natural initial condition of \rho^{\rm{L}}_{\rm{S}}(0)=|8\rangle\langle 8| [13, 62] and apply the HEOM [36-54] to calculate the trapping time. FIG. 6(a) demonstrates that the total trapping time \langle t\rangle decreases significantly by five orders of magnitude as the reorganization energy increases from 0 cm^{-1} to 1 cm^{-1}. In the weak dissipation regime, the partial trapping time \langle t\rangle^{\rm{E}}_{\rm{P}} from the secular Redfield equation provides an accurate estimation of the total trapping time \langle t\rangle from the HEOM. The approximate 1/\Gamma-scaling,  \begin{eqnarray} \langle t\rangle^{\rm{E}}_{\rm{P}} \approx t_0 + \frac{t_1}{\lambda+\lambda_1}+\frac{t_2}{\lambda+\lambda_2}+\frac{t_3}{\lambda+\lambda_3} \label{eq_41} \end{eqnarray} (41)  FIG. 6 The trapping time vs. the reorganization energy in the 8-site FMO (k_t=1 ps^{-1}). (a) The initial system state is at \rho^{\rm{L}}_{\rm{S}}(0)=|8\rangle\langle 8|. The circles and solid line are the results of \langle t\rangle from the HEOM and \langle t\rangle^{\rm{E}}_{\rm{P}} from the secular Redfield equation. The dashed and dotted-dashed lines are the results of Eqs. (41) and (42). (b) The solid and dashed lines are the results of the initial system state in Eq.(43) with and without site-site coherence can reliably describe the partial trapping time over a broad range of the dissipation strength (\lambda$$ <$100 cm$^{-1}$). Here the three characteristic reorganization energies are $\lambda_1$=7.55$\times$10$^{-5}$ cm$^{-1}$, $\lambda_2$=5.23$\times$10$^{-3}$ cm$^{-1}$, and $\lambda_3$=0.858 cm$^{-1}$. Eq.(41) is consistent with the fact that four sites form a major energy transfer pathway, BChls 8$\rightarrow$(1, 2)$\rightarrow$3 [13, 55, 62-64]. With $\lambda_1, \lambda_2$$\ll$$\lambda_3$, the exact $1/\Gamma$-scaling,

 $\begin{eqnarray} \langle t\rangle^{\rm{E}}_{\rm{P}} \sim \left(t_0+\frac{t_3}{\lambda_3}\right) +\frac{t_1+t_2}{\lambda} \label{eq_42} \end{eqnarray}$ (42)

is extracted in the weak dissipation regime of 0.01 cm$^{-1}$$<$$\lambda$$< 0.2 cm^{-1}. On the opposite side, the optimal initial system state in Eq.(29) is determined by the secular Redfield equation, from which the efficiency of the coherent energy transfer (\lambda$$\rightarrow$0) is optimized. For simplicity, we only present one example,

 $\begin{eqnarray*} \rho^{\rm{L}}_{\rm{S}}(0) = \left( {\begin{array}{*{20}{c}} 0.0081 & 0.0146 &-0.0444 &-0.0127 &-0.0020 &-0.0013 &-0.0015 &-0.0010 \\ 0.0146 & 0.0303 &-0.1252 &-0.0394 &-0.0058 &-0.0042 &-0.0059 &-0.0018 \\ -0.0444 &-0.1252 & 0.8378 & 0.2675 & 0.0392 & 0.0259 & 0.0420 & 0.0045 \\ -0.0127 &-0.0394 & 0.2675 & 0.1065 & 0.0184 & 0.0061 & 0.0242 & 0.0016 \\ -0.0020 &-0.0058 & 0.0392 & 0.0184 & 0.0043 &-0.0009 & 0.0055 & 0.0002 \\ -0.0013 &-0.0042 & 0.0259 & 0.0061 &-0.0009 & 0.0037 &-0.0011 & 0.0002 \\ -0.0015 &-0.0059 & 0.0420 & 0.0242 & 0.0055 &-0.0011 & 0.0091 & 0.0003 \\ -0.0010 &-0.0018 & 0.0045 & 0.0016 & 0.0002 & 0.0002 & 0.0003 & 0.0001 \end{array}} \right) \end{eqnarray*}$ (43)

The trapping time $\langle t\rangle$ under this artificially designed optimal initial state $\rho^{\rm{L}}_{\rm{S}}(0)$ is then calculated using the HEOM. As shown in FIG. 6(b), $\langle t\rangle$ in the coherent limit ($\lambda$$\rightarrow0) is minimized to \langle t\rangle=2.64 ps, significantly smaller than that from \rho^{\rm{L}}_{\rm{S}}(0)=|8\rangle\langle 8|. As the reorganization energy \lambda is increased from 0 to 100 cm^{-1}, the trapping time is increased very weakly to \langle t\rangle=2.77 ps, resisting the influence of dynamic disorders. To illustrate the relevance of an appropriate site-site coherence in the optimal state, we remove the off-diagonal elements in Eq.(43) and \langle t\rangle is changed to a decreasing function of \lambda (see FIG. 6(b)). In fact, the NEET behavior can be observed even with all the initial population localized at the trapping site, BChl 3. Thus, the optimized coherent energy transfer requires an appropriate initial site-site coherence constrained by Eq.(29). Ⅴ. SUMMARY In this paper, we extend our previous studies of efficiency optimization [12, 14, 15] to probing the mechanism of the noise-enhanced energy transfer (NEET) and engineering the initial system state to maximize the efficiency of coherent energy transfer from a conceptual point of view. In the weak dissipation limit, the trapping time \langle t\rangle of an energy transfer network with a slow trapping process can be reliably approximated by the survival time \langle t\rangle^{\rm{E}}_{\rm{P}} of an exciton spanning only inside its population subspace. Following a simple but accurate description of quantum dissipation from the secular Redfield equation, we express \langle t\rangle^{\rm{E}}_{\rm{P}} in a rational polynomial of the dissipation strength \Gamma and formulate the general requirement of the NEET. If the determinant of the trapping superoperator is zero in the exciton population subspace (\mathit{\rm{Det}}[{\cal L}^{\rm{E}}_{{\rm{t}}; {\rm{P}}}]=0), a rigorous trapping-free subspace \Phi_\perp is formed and the trapping time follows the exact 1/\Gamma-scaling, \langle t\rangle$$\sim$$1/\Gamma, in the limit of \Gamma$$\rightarrow$0. Under a weaker condition that $\mathit{\rm{Det}}[{\cal L}^{\rm{E}}_{{\rm{t}}; {\rm{P}}}]$ is nonzero but sufficiently small, the trapping time can be still approximated as $\langle t\rangle$$\sim$$A$+$B/\Gamma$ over a certain range of dissipation strengths. In contrast, the initial system state can be tuned finely to suppress the $1/\Gamma$-scaling of the NEET, globally or locally maximizing the transfer efficiency in the coherent limit ($\Gamma$$\rightarrow$0). A mathematical procedure is proposed by us to fulfill the optimization constraint in Eq.(29), showing an infinite possibility of solutions.

Our theoretical predictions of the $1/\Gamma$-scaling and the optimal initial system states are verified in the four examples of the energy transfer networks: a biased two-site system, a symmetric three-site branching system, a homogeneous 1D chain, and an 8-chromophore FMO monomer. We observe that the NEET can always appear approximately due to a slow trapping process or exactly due to a symmetric system Hamiltonian compatible with trapping. The approximate $1/\Gamma$-scaling induced by the former condition becomes more pronounced with a larger energy mismatch and a longer spatial extension. For each model system, the optimal initial system states to maximize the coherent transfer efficiency are successfully obtained following our procedure.

The studies in this paper confirm that the NEET is a universal behavior in the energy transfer network, which can be applied to other irreversible quantum dynamic systems such as a quantum heat engine. Here we focus on the dephasing rate and reorganization energy, but the concept of the dissipation strength can be extended to other system- and bath-related parameters such as temperature. Our calculation of the optimal initial system states shows that quantum (site-site) coherence to accelerate energy transfer must be tuned finely to be compatible with dissipation and trapping. Although such a fine tuning is difficult to be achieved in natural systems, it is experimentally accessible in precisely-controlled artificial quantum devices, which will be an interesting problem to be explored in the future [65].

Ⅵ. ACKNOWLEDGEMENTS

The work reported was supported by the National Natural Science Foundation of China (No.21573195) and the Ministry of Science and Technology of China (MOST-2014CB921203).

Reference
 [1] A. Chenu, and G. D. Scholes, Annu. Rev. Phys. Chem. 66 , 69 (2015). DOI:10.1146/annurev-physchem-040214-121713 [2] R. Kosloff, and A. Levy, Annu. Rev. Phys. Chem. 65 , 365 (2014). DOI:10.1146/annurev-physchem-040513-103724 [3] D. Z. Xu, and J. S. Cao, Front. Phys. 11 , 110308 (2016). DOI:10.1007/s11467-016-0540-2 [4] P. Brumer, and M. Shapiro, Ann. Rev. Phys. Chem. 43 , 257 (1992). DOI:10.1146/annurev.pc.43.100192.001353 [5] S. J. Glaser, U. Boscain, T. Calarco, C. P. Koch, W. Köckenberger, R. Kosloff, I. Kuprov, B. Luy, S. Schirmer, T. Schulte-Herbrüggen, D. Sugny, and F. K. Wilhelm, Eur. Phys. J. D 69 , 279 (2015). DOI:10.1140/epjd/e2015-60464-1 [6] K. M. Gaab, and C. J. Bardeen, J. Chem. Phys. 121 , 7813 (2004). DOI:10.1063/1.1786922 [7] A. Olaya-Castro, C. F. Lee, F. F. Olsen, and N. F. Johnson, Phys. Rev. B 78 , 085115 (2008). DOI:10.1103/PhysRevB.78.085115 [8] B. M. Plenio, and F. S. Huelga, New J. Phys. 10 , 113019 (2008). DOI:10.1088/1367-2630/10/11/113019 [9] F. Caruso, A. W. Chin, A. Datta, S. F. Huelga, and M. B. Plenio, J. Chem. Phys. 131 , 105106 (2009). DOI:10.1063/1.3223548 [10] F. Caruso, New J. Phys. 16 , 055015 (2014). DOI:10.1088/1367-2630/16/5/055015 [11] J. S. Cao, and R. J. Silbey, J. Phys. Chem. A 113 , 13825 (2009). [12] J. L. Wu, F. Liu, Y. Shen, J. S. Cao, and R. J. Silbey, New J. Phys. 12 , 105012 (2010). DOI:10.1088/1367-2630/12/10/105012 [13] J. Moix, J. L. Wu, P. F. Huo, D. Coker, and J. S. Cao, J. Phys. Chem. Lett. 2 , 3045 (2011). DOI:10.1021/jz201259v [14] J. L. Wu, F. Liu, J. Ma, R. J. Silbey, and J. S. Cao, J. Chem. Phys. 137 , 174111 (2012). DOI:10.1063/1.4762839 [15] J. L. Wu, R. J. Silbey, and J. S. Cao, Phys. Rev. Lett. 110 , 200402 (2013). DOI:10.1103/PhysRevLett.110.200402 [16] M. Mohseni, P. Rebentrost, S. Lloyd, and A. AspuruGuzik, J. Chem. Phys. 129 , 174106 (2008). DOI:10.1063/1.3002335 [17] P. Rebentrost, M. Mohseni, I. Kassal, S. Lloyd, and A. Aspuru-Guzik, New J. Phys. 11 , 033003 (2009). DOI:10.1088/1367-2630/11/3/033003 [18] I. Kassal, and A. Aspuru-Guzik, New J. Phys. 14 , 053041 (2012). DOI:10.1088/1367-2630/14/5/053041 [19] R. de J. Léon-Montiel, I. Kassal, and J. P. Torres, J. Phys. Chem. B 118 , 10588 (2014). DOI:10.1021/jp505179h [20] M. Mohseni, A. Shabani, S. Lloyd, and H. Rabitz, J. Phys. Chem. 140 , 035102 (2014). DOI:10.1063/1.4856795 [21] L. Mühlbacher, and U. Kleinekathöfer, J. Phys. Chem. B 116 , 3900 (2012). DOI:10.1021/jp301444q [22] E. K. Irish, R. Gomez-Bombarelli, and B. W. Lovett, Phys. Rev. A 90 , 012510 (2014). DOI:10.1103/PhysRevA.90.012510 [23] U. Weiss, Quantum Dissipative Systems, New Jersey: World Scientific Publishing Company, (2008). [24] H. Breuer and F. Petruccione, The Theory of Open Quantum Systems, New York: Oxford University Press, (2002). [25] T. Förster, Ann. Phys. (Leipzig) 437 , 55 (1948). DOI:10.1002/(ISSN)1521-3889 [26] P. Hanggi, P. Talkner, and M. Borkovec, Rev. Mod. Phys. 62 , 251 (1990). DOI:10.1103/RevModPhys.62.251 [27] B. J. Schwartz, Annu. Rev. Phys. Chem. 54 , 141 (2003). DOI:10.1146/annurev.physchem.54.011002.103811 [28] S. Mukamel, Principles of Nonlinear Optical Spectroscopy, New York: Oxford University Press, (1995). [29] V. May and K. Oliver, Charge and Energy Transfer Dynamics in Molecular Systems, Weinheim: Wiley-VCH, (2004). [30] H. Haken, and P. Reineker, Z. Phys. 249 , 253 (1972). DOI:10.1007/BF01400230 [31] H. Haken, and G. Strobl, Z. Phys. 262 , 135 (1973). DOI:10.1007/BF01399723 [32] S. Nakajima, Prog. Theor. Phys. 20 , 948 (1958). DOI:10.1143/PTP.20.948 [33] R. Zwanzig, J. Chem. Phys. 33 , 1338 (1960). DOI:10.1063/1.1731409 [34] A. G. Redfield, IBM J. Res. Dev. 19 , 1 (1957). [35] G. Lindblad, Commun. Math. Phys. 48 , 119 (1976). DOI:10.1007/BF01608499 [36] Y. Tanimura, and R. Kubo, J. Phys. Soc. Jpn. 58 , 101 (1989). DOI:10.1143/JPSJ.58.101 [37] A. Ishizaki, and Y. Tanimura, J. Phys. Soc. Jpn. 74 , 3131 (2005). DOI:10.1143/JPSJ.74.3131 [38] Y. Tanimura, J. Chem. Phys. 141 , 044114 (2014). DOI:10.1063/1.4890441 [39] Y. Tanimura, J. Chem. Phys. 142 , 144110 (2015). DOI:10.1063/1.4916647 [40] Y. A. Yan, F. Yang, Y. Liu, and J. S. Shao, Chem. Phys. Lett. 395 , 216 (2004). DOI:10.1016/j.cplett.2004.07.036 [41] J. S. Shao, Chem. Phys. 322 , 187 (2006). DOI:10.1016/j.chemphys.2005.08.007 [42] Y. Zhou, and J. S. Shao, J. Chem. Phys. 128 , 034106 (2008). DOI:10.1063/1.2818095 [43] R. X. Xu, P. Cui, X. Q. Li, Y. Mo, and Y. J. Yan, J. Chem. Phys. 122 , 041103 (2005). DOI:10.1063/1.1850899 [44] J. Hu, M. Luo, F. Jiang, R. X. Xu, and Y. J. Yan, J. Chem. Phys. 134 , 244106 (2011). DOI:10.1063/1.3602466 [45] Y. J. Yan, J. Chem. Phys. 140 , 054105 (2014). DOI:10.1063/1.4863379 [46] R. X. Xu, Y. Liu, H. D. Zhang, and Y. J. Yan, Chin. J. Chem. Phys. 30 , 395 (2017). DOI:10.1063/1674-0068/30/cjcp1706123 [47] J. M. Moix, and J. S. Cao, J. Chem. Phys. 139 , 134106 (2013). DOI:10.1063/1.4822043 [48] C. Y. Hsieh, and J. S. Cao, J. Chem. Phys. 148 , 014103 (2018). DOI:10.1063/1.5018725 [49] C. Y. Hsieh, and J. S. Cao, J. Chem. Phys. 148 , 014104 (2018). DOI:10.1063/1.5018726 [50] Q. Shi, L. P. Chen, G. J. Nan, R. X. Xu, and Y. J. Yan, J. Chem. Phys. 130 , 084105 (2009). DOI:10.1063/1.3077918 [51] H. Liu, L. L. Zhu, S. M. Bai, and Q. Shi, J. Chem. Phys. 140 , 134106 (2014). DOI:10.1063/1.4870035 [52] Z. F. Tang, X. L. Ouyang, Z. H. Gong, H. B. Wang, and J. L. Wu, J. Chem. Phys. 143 , 224112 (2015). DOI:10.1063/1.4936924 [53] C. R. Duan, Z. F. Tang, J. S. Cao, and J. L. Wu, Phys. Rev. B 95 , 214308 (2017). DOI:10.1103/PhysRevB.95.214308 [54] C. R. Duan, Q. L. Wang, Z. F. Tang, and J. L. Wu, J. Chem. Phys. 147 , 164112 (2017). DOI:10.1063/1.4997669 [55] J. L. Wu, Z. F. Tang, Z. H. Gong, J. S. Cao, and S. Mukamel, J. Phys. Chem. Lett. 6 , 1240 (2015). DOI:10.1021/acs.jpclett.5b00227 [56] The terms of 'coherent' and 'incoherent' usually refer to the local basis representation. The evolutions in the exciton population and coherence subspaces correspond to the coherent and incoherent motions of local sites, respectively. [57] J. L. Wu, and J. S. Cao, Adv. Chem. Phys. 146 , 329 (2011). [58] P. W. Anderson, Phys. Rev. 109 , 1492 (1958). DOI:10.1103/PhysRev.109.1492 [59] J. M. Moix, M. Khasin, and J. S. Cao, New J. Phys. 15 , 085010 (2013). DOI:10.1088/1367-2630/15/8/085010 [60] R. E. Fenna, and B. W. Matthews, Nature (London) 258 , 573 (1975). DOI:10.1038/258573a0 [61] D. E. Tronrud, J. Z. Wen, L. Gay, and R. E. Blankenship, Photosynth. Res. 100 , 79 (2009). DOI:10.1007/s11120-009-9430-6 [62] M. Schmidt am Busch, F. Müh, M. E. A. Madjet, and T. Renger, J. Phys. Chem. Lett. 2 , 93 (2011). DOI:10.1021/jz101541b [63] J. Adolphs, and T. Renger, Biophys. J. 91 , 2778 (2006). DOI:10.1529/biophysj.105.079483 [64] A. Ishizaki, and G. R. Fleming, Proc. Natl. Acad. Sci. USA 106 , 17255 (2009). DOI:10.1073/pnas.0908989106 [65] A. Potočnik, A. Bargerbos, F. A. Y. N. Schröder, S. A. Khan, M. C. Collodo, S. Gasparinetti, Y. Salathé, C. Creatore, C. Eichler, H. E. Türeci, A. W. Chin, and A. Wallraff, Nat. Commun. 9 , 904 (2018). DOI:10.1038/s41467-018-03312-x

a. 浙江大学物理系, 杭州 310027;
b. 麻省理工学院化学系, 马萨诸塞州, 剑桥 02139