Efficient quantum simulation of open quantum dynamics at various Hamiltonians and spectral densities

Abstract

Simulation of open quantum dynamics for various Hamiltonians and spectral densities are ubiquitous for studying various quantum systems. On a quantum computer, only log₂ N qubits are required for the simulation of an N -dimensional quantum system, hence simulation in a quantum computer can greatly reduce the computational complexity compared with classical methods. Recently, a quantum simulation approach was proposed for studying photosynthetic light harvesting [npj Quantum Inf. 4, 52 (2018)]. In this paper, we apply the approach to simulate the open quantum dynamics of various photosynthetic systems. We show that for Drude–Lorentz spectral density, the dimerized geometries with strong couplings within the donor and acceptor clusters respectively exhibit significantly improved efficiency. We also demonstrate that the overall energy transfer can be optimized when the energy gap between the donor and acceptor clusters matches the optimum of the spectral density. The effects of different types of baths, e.g., Ohmic, sub-Ohmic, and super-Ohmic spectral densities are also studied. The present investigations demonstrate that the proposed approach is universal for simulating the exact quantum dynamics of photosynthetic systems.

Graphical abstract

Keywords

nuclear magnetic resonance / quantum simulation / open quantum system

1 Introduction

As commonly understood, the efficiency of exciton energy transfer (EET) in natural photosynthesis is close to unity [1–3]. Because of the discovery of EET with coherent features, the role of quantum coherence in EET efficiency has become a research hotspot in the past two decades [4–13]. Pigment-protein complexes in photosynthesis are essentially open quantum systems. Since the couplings between the system and the environment are of the same order of magnitude as the couplings within the system [2, 3], non-Markovian features arise and make simulating the open quantum dynamics therein difficult [14–18]. Even so, the research of open quantum systems plays an important role in many fields [19, 20]. So far, some theoretical methods have been proposed to effectively simulate EET in photosynthesis [3], such as the numerically-exact hierarchical equation of motion (HEOM) [21–28], the quantum jump approach [29, 30], the small-polaron quantum master equation [31], the modified Redfield theory and its coherent generalization [32, 33]. Among these methods simulating the EET, the HEOM yields exact quantum dynamics in the whole parameter regime, e.g., the Förster regime and the Redfield regime [3, 22]. It is helpful to revealing the role of quantum coherence in optimizing the photosynthetic EET [5] and clearly elucidating the design principals of artificial light-harvesting devices [34–40].

However, despite the fact that the HEOM has been widely used in the study of open quantum dynamics, including EET in natural photosynthesis, in the case of large dimensions and complex spectral densities, the numerical overhead becomes very large. Recently, we proposed a novel experimental approach to simulate EET in photosynthesis [41]. We have proved that this approach is exact and efficient in high-temperature regime. We generate a large number of realizations driven by random Hamiltonians, and by averaging over the ensemble we obtain a density matrix whose dynamics is subject to decoherence. As a demonstration, we adopted a prototype in Ref. [42] and compared the results of nuclear magnetic resonance (NMR) simulation and HEOM calculation under Drude– Lorentz noise. We showed that it is valid to efficiently simulate the exact quantum dynamics in the photosynthetic EET by using NMR if the number of random realizations is sufficiently large.

As we know, for systems with large dimension or complex spectral densities, the HEOM requires a huge amount of computation resources. For example, to simulate an N - level system, the computation cost of an \mathcal{N} -layer HEOM scales exponentially in \mathcal{N}(\mathcal{N}\le N) [43]. However, in the quantum simulation [44, 45], because the quantum dynamics of N states can be simulated by using log₂ N qubits, the computation cost is a polynomial of N. Therefore, this quantum simulation can effectively reduce the computational complexity. And we know that qubits can build different quantum gates, and quantum gates are also the basic components of quantum computing and quantum communication [46–48].

In 2013, it was demonstrated that the efficiency of energy transfer can be improved when there is strong coupling within donor and acceptor pairs by studying energy transfer in a linear-tetramer model [42]. In the same year, del Rey et al. proposed a design principle called phonon antenna. By spectrally sampling optimum in their local environmental fluctuations, the coherence between internal pigments can affect and optimize the way excitation flows [49]. And the strong coupling to an under-damped vibrational mode can help the photosynthetic complex to overcome the energy barrier between the donor and acceptor, and thus increase the efficiency [50].

In Ref. [41], only a quantum simulation with a specific Hamiltonian and Drude–Lorentz noise was demonstrated. The clustered geometry was proven to be optimal for a broad parameter regime by the coherent modified Redfield theory [42]. However, the theory is shown to break down when the reorganization energy is much larger than the intra-system coupling [33, 51]. A natural question arises: does the above discovery still hold in a broad range of parameters by a numerically-exact approach? On the other hand, since both modifying the geometry and spectral-sampling in local environmental fluctuations can optimize energy flow in EET respectively, it is quite natural to ask whether light harvesting can be further optimized if both means have been applied? In this paper, we apply this numerically-exact approach to efficiently verify the hypothesis that we can simultaneously optimize the energy transfer by the geometry and environment. Furthermore, for spectral densities beyond Drude–Lorentz type, the computation complexity of the HEOM will be significantly increased, when there are more exponential components in the correlation function of the bath [43, 52]. In this paper, we also show that the quantum simulation approach is applicable in photosynthetic systems with complicated environments, which are not tractable by the conventional approach.

The structure of this paper is organized as follows. In Section 2, we give a brief introduction to the HEOM and the approach for quantum simulation in NMR. It combines the bath-engineering technique [53, 54] and the gradient ascent pulse engineering (GRAPE) algorithm [55, 56], which will be introduced in Appendixes B and C, respectively. In Section 3.1, we first consider the simulation of the EET dynamics in a linear-tetramer model with different geometries. After the geometry is optimized, we also consider how to improve the EET efficiency from the aspect of the bath, which will be discussed further in Appendix E. Through this specific model, we confirm the discoveries in Refs. [42, 49]. Besides, in Section 3.2 we also consider the effects of different types of spectral densities on the EET dynamics, which was not analyzed in Ref. [41]. In Section 3.3, by N = 8, we show how to generalize our quantum simulation approach to a case beyond 2 qubits. In Section 3.4, we compare the computational complexities of NMR simulation and HEOM, and analyze the possible errors in NMR simulations in Section 3.5. Finally, we discuss the prospect and conclusions in Section 4. In Appendixes A and D, we introduce the initialization and tomography in NMR, respectively.

2 Theory for quantum simulation

In this section, we introduce how to simulate the EET dynamics of photosynthesis in NMR systems. In Ref. [41], we simulated the influence of noise by using the bath-engineering technique [53, 54], and realized the evolution of quantum states by using the GRAPE algorithm [55, 56]. We summarize these two techniques in Appendixes B and C. In the following, we shall give a brief introduction to our model system and the HEOM theory for photosynthetic EET and the process of quantum simulation.

2.1 Model photosynthetic system

Following Ref. [41], we use the tetramer model, which contains four chlorophyll molecules, for the quantum simulation. As shown in Fig. 1(a), the left pair of chlorophylls act as donors and the right pair as acceptors. For molecules, their electronic structures are generally studied [57]. Here, we adopt the Frenkel-exciton Hamiltonian [2, 58]

\begin{array}{}\text{(1)}& H_{\text{EET}}={\displaystyle \sum _{i=1}^4}{\epsilon }_i|i\rangle \langle i|+{\displaystyle \sum _{1\le i\ne j\le 4}}{J}_{ij}|i\rangle \langle j|\end{array}

to describe the dynamics of EET in photosynthesis. Here |i〉 (i = 1, 2, 3, 4) represents that only ith molecule is in the excited state while the others are in the ground state. ε_i is the site energy of ith exciton. For simplicity, the electronic interaction between ith and jth excitons is given by the dipole-dipole interaction as

\begin{array}{}\text{(2)}& J_{ij}=\frac{1}{4\pi {\epsilon }_0{r}_{ij}^3}[{\overrightarrow{\mu }}_i\cdot {\overrightarrow{\mu }}_j-3({\overrightarrow{\mu }}_i\cdot {\widehat{r}}_{ij})({\overrightarrow{\mu }}_j\cdot {\widehat{r}}_{ij})]\mathrm{,}\end{array}

where {\overrightarrow{r}}_{ij}=r_{ij}{\widehat{r}}_{ij} is the displacement vector from site i to site j, {\overrightarrow{\mu }}_i the transition dipole of site i, ε₀ the vacuum permittivity. In numerical simulations, we take µ_j = 7.75 D and r∈[6, 14] Å, which are typical in natural photosynthesis.

In photosynthetic complexes, the energy transfer is assisted by the interaction between the system and the bath, which can be described as

\begin{array}{}\text{(3)}& H_{\text{SB}}={\displaystyle \sum _{i\mathrm{,}k}{g}_{ik}{V}_j\left(a_{ik}^{†}+a_{ik}\right)\mathrm{.}}\end{array}

Here, V_j = |i〉〈i|, a_{ik}^{†} (a_ik) is the creation (annihilation) operator of the kth phonon mode of the ith molecule, and g_ik represents the coupling strength. H_SB is the main cause of energy relaxation in photosynthetic systems. All the information about the couplings between system and environment can be given by the spectral density, i.e.,

\begin{array}{}\text{(4)}& G_{\text{EET}}(\omega )={\displaystyle \sum _k{g}_{ik}^2\delta (\omega -{\omega }_k)\mathrm{.}}\end{array}

The spectral density plays an important role in the optimization of EET in photosynthesis. It has been shown that the energy transfer can be improved by adjusting the parameters of photosynthetic system so that the energy gap matches the optimal frequency of the spectral density [49].

2.2 Hierarchical equation of motion method

The exact EET dynamics can be given by the HEOM [22, 59]

\begin{array}{}\text{(5)}& \begin{array}{ll}\frac{\partial }{\partial t}{\sigma }_{\overrightarrow{n}}=\hfill & (-\text{i}{\mathcal{L}}_e+{\displaystyle \sum _j}{n}_j{\gamma }_j){\sigma }_{\overrightarrow{n}}+{\displaystyle \sum _j}{\text{Φ}}_j{\sigma }_{{\overrightarrow{n}}_{j+}}+{\displaystyle \sum _j}{n}_j{\text{Θ}}_j{\sigma }_{{\overrightarrow{n}}_{j-}}\mathrm{,}\hfill \end{array}\end{array}

where {\sigma }_{\overrightarrow{n}} and {\sigma }_{{\overrightarrow{n}}_{j\pm }} are the auxiliary density matrices with \overrightarrow{n} = (n₁, n₂, · · · , n_j, · · ·) and {\overrightarrow{n}}_{j\pm } = (n₁, n₂, · · · , n_j ± 1, · · ·), and n_j ’s non-negative integers, {\sigma }_{\overrightarrow{0}} = ρ the reduced density matrix of photosynthetic system, {\mathcal{L}}_e the Liouville superoperator of H_EET. Besides, Φ_j =\text{i}{V}_j^{\times } and {\text{Θ}}_j=\text{i}(\frac{2{\lambda }_j}{\beta {\hslash }^2}{V}_j^{\times }-\text{i}\frac{{\lambda }_j}{\hslash }{V}_j^{\circ }) with V_j^{\times }{\sigma }_{\overrightarrow{n}}=V_j{\sigma }_{\overrightarrow{n}}-{\sigma }_{\overrightarrow{n}}{V}_j and V_j^{\circ }{\sigma }_{\overrightarrow{n}}=V_j{\sigma }_{\overrightarrow{n}}+{\sigma }_{\overrightarrow{n}}{V}_j , where V_j = |j〉〈j|. For Drude–Lorentz spectral density, i.e.,

\begin{array}{}\text{(6)}& G_{\text{EET}}(\omega )=\frac{2{\lambda }_j{\gamma }_j\omega }{{\omega }^2+{\gamma }_j^2}\end{array}

with λ_j the reorganization energy and γ_j the relaxation rate. For a generic spectral density, it can be decomposed into a summation of Lorentzian form [52, 60].

2.3 The process of quantum simulation

In the quantum simulation of photosynthetic EET with NMR system in Ref. [41], there are three steps: We first perform the theoretical simulations by random Hamiltonians on the ensemble to obtain the optimal parameters for the GRAPE simulation. Then, we perform the GRAPE simulation to obtain optimal pulse sequence for the experimental realization. Finally, one can then perform the experiment based on the above theoretical simulations. If the experimental realizations are ideal, the experiments will faithfully reproduce the results of the GRAPE simulation. Since in Ref. [41], the experimental demonstrations for the Drude–Lorentz spectral density coincide with those by the HEOM, we will focus on the first two steps in the following discussions. In Appendix, we describe these three steps in detail.

For photosynthetic complexes with N chlorophylls, there are only N single-excitation states in the process of energy transfer. As a result, only log₂ N qubits are required for quantum simulation and thus two qubits are needed to simulate the tetramer model. In this case, the single-excitation states in photosynthesis are encoded as two-qubit product states as |1〉 = |00〉, |2〉 = |01〉, |3〉 = |10〉, and |4〉 = |11〉. As shown in Fig. 1(b), we regard H and C nuclear spins as the two qubits. The Hamiltonian H_NMR implemented in NMR simulations is H_NMR = H_EET/C with a scaling factor

\begin{array}{}\text{(7)}& C=3\times {10}^9\mathrm{.}\end{array}

The total Hamiltonian for simulating EET process is

\begin{array}{}\text{(8)}& H(t)=H_{\text{NMR}}+H_{\text{PDN}}\mathrm{.}\end{array}

H_PDN is introduced to mimic the effects of the local baths in photosynthesis. In order to exactly simulate the bath’s effect, we require that χ(Ct) = Re[g(t)] as shown in Appendix B, where the decoherence factor χ(t) is the Fourier transform of the noise’s correlation function in the quantum simulation, g(t) is the lineshape function of photosynthesis. By the bath-engineering technique [53, 54, 61], we will show that it is possible to simulate various types of spectral densities in Section 3.2.

3 Numerical calculation and analysis

In Ref. [41], only the EET dynamics for a specific Hamiltonian under the Drude–Lorentz noise was demonstrated. In this section, we will use the above quantum simulation approach to investigate the dynamics for different geometries and various types of spectral densities. In addition, we will compare the computational complexity of the quantum simulation and the HEOM, and analyze the errors in the quantum simulation.

3.1 Simultaneous optimization by geometry and bath

In Ref. [42], the effects of the geometry on the energy transfer was investigated. In the study, it was revealed that the dimerized structure explores coherent relaxation to promote the energy transfer within the dimer. Since the coupling between two sites is significantly dependent on the distance between two chlorophylls, different geometries correspond to different Hamiltonians. In this sub-section, we analyze the EET dynamics for different Hamiltonians with the quantum simulation approach.

In Fig. 1(a), we show the comparison among the HEOM, the quantum simulation and the GRAPE simulation for r = 13.4 Å. Therein, all four chlorophylls are equally spaced and the Hamiltonian for the NMR quantum simulation H_NMR = H_EET/C is

\begin{array}{}\text{(9)}& \frac{H_{\text{NMR}}}{2\pi }=\left(\begin{array}{cccc}130& 1.2608& 0.1612& 0.0474\\ 1.2608& 129& 1.3190& 0.1612\\ 0.1612& 1.3190& 123& 1.2608\\ 0.0474& 0.1612& 1.2608& 122\end{array}\right)\text{kHz}\mathrm{.}\end{array}

In Fig. 1(b), we set r = 11.3 Å and thus the distance between the two donors (acceptors) is slightly smaller than the distance between the central two sites. In this case, the donors (acceptors) form a dimer. The corresponding Hamiltonian reads

\begin{array}{}\text{(10)}& \frac{H_{\text{NMR}}}{2\pi }=\left(\begin{array}{cccc}130& 2.1025& 0.1283& 0.0474\\ 2.1025& 129& 0.5759& 0.1283\\ 0.1283& 0.5759& 123& 2.1025\\ 0.0474& 0.1283& 2.1025& 122\end{array}\right)\text{kHz}\mathrm{.}\end{array}

In Fig. 1(c), we further reduce the intra-dimer distance to an even smaller value, i.e., r = 8 Å. Thus, the Hamiltonian is

\begin{array}{}\text{(11)}& \frac{H_{\text{NMR}}}{2\pi }=(\begin{array}{cccc}130& 5.9251& 0.0926& 0.0474\\ 5.9251& 129& 0.2195& 0.0926\\ 0.0926& 0.2195& 123& 5.9251\\ 0.0474& 0.0926& 5.9251& 122\end{array})\text{kHz}\mathrm{.}\end{array}

In order to obtain the above Hamiltonians, we adopt the site energies ε₁ = 13 000 cm⁻¹, ε₂ = 12 900 cm⁻¹, ε₃ = 12 300 cm⁻¹, and ε₄ = 12 200 cm⁻¹, which are typical in natural photosynthetic systems, for the corresponding Hamiltonians of photosynthesis. Except that J₁₄ is invariant, different distances r’s correspond to a set of different coupling terms J_ij ’s (i, j = 1, 2, 3, 4) as determined by Eq. (2).

In addition, we assume the spectral density of Drude– Lorentz form as given in Eq. (6), with the optimal frequency γ_j. In the simulations, we assume identical phonon relaxation rates γ_NMR = γ_EET/C = 2π ×50 Hz and identical reorganization energies λ_NMR = λ_EET/C = 2π × 2 Hz for all local baths. In addition, the influence of temperature is also very important [62]. Since the results of quantum simulation coincide with those of the HEOM only at high temperatures, we take T_EET = 3 × 10⁶ K and T_NMR = 10⁻³ K in our numerical calculations. In this way, we compare the results of the quantum simulation, and the GRAPE simulation, and the HEOM in Fig. 1. As shown in Fig. 1, both the results of the random ensemble and GRAPE are in good agreement with the results of the HEOM simulation. Therefore, we show that our approach can be applied to photosynthetic systems with different geometries, beyond the linear one with equal-spacing.

As shown in Fig. 1(b), the energy transfer is the fastest for the case with r = 11.3 Å. The dimerized geometry explores the coherent relaxation within the donors to accelerate the energy transfer. However, over-dimerization in the geometry significantly reduces the energy transfer rate, because it enlarges the energy gap between the donor and acceptor clusters, cf. Fig. 1(c). These numerically-exact simulations are consistent with those discoveries obtained by the approximate theory in Ref. [42]. We notice that there are some small differences between the HEOM and the quantum simulations, which rely on the assumption that the ensemble average is equivalent to the time average. Therefore, we verify the assumption in Fig. 1(d). As the number of random realizations in the ensemble increases, the results of quantum simulation approach closer and closer to those of the HEOM. Besides, we also observe the deviations of the GRAPE simulations from the quantum simulations, the possible causes of which will be discussed at the end of this section. In this regard, both the quantum simulations and GRAPE simulations successfully reproduce the coherent oscillations at the short-time regime and the incoherent relaxation at the long-time regime, and it is valid to exactly simulate the open quantum dynamics with a generic Hamiltonian by the above approach.

In Ref. [49], it was shown that the optimization of the EET can be achieved when the energy gap of the system matches the optimum of the spectral density. Inspired by this discovery, we subsequently consider the effects of the bath on the EET dynamics for a fixed Hamiltonian with r = 11.3 Å, which was shown to be optimal among different geometries in Fig. 1. Since the peak of the Drude–Lorentz spectral density is located at γ_NMR, we fix λ_NMR and simulate the EET dynamics for a broad range of γ_NMR. Here, in Fig. 2, we only demonstrate the population dynamics for three typical parameters, i.e., γ_NMR/(2π) = 0.5 kHz, 2.668 kHz, and 7 kHz. Obviously, the optimal relaxation rate of the bath is {\gamma }_{\text{NMR}}^{\text{opt}} = 2π × 2.668 kHz as it requires the shortest time to achieve equal populations. The physical mechanism can be described by the energy diagram of the system. Due to their strong couplings, sites 1 and 2 form the donor cluster, while sites 3 and 4 form the acceptor cluster. Since the donor cluster is weakly coupled with the acceptor cluster, the EET between the two clusters can be described by the Förster theory. Therefore, the inter-cluster EET is optimized when the energy gap between the lower-energy eigen-state of the donor cluster and the higher-energy eigen-state of the acceptor cluster matches the optimum of the spectral density. The details about this physical mechanism is elucidated in Appendix E.

3.2 Application to complicated spectral densities

Since a general power-law form spectral density can describe an extremely-large number of physical environments [63, 64], we use it to show the universal applicability of the quantum simulation approach. Following Refs. [17, 65], we analyze the EET dynamics for three types of spectral densities, namely Ohmic, sub-Ohmic, and super-Ohmic spectral densities, which are expressed in a unified manner as

\begin{array}{}\text{(12)}& J(\omega )=\frac{\lambda \omega }{\text{Γ}(s)}{\left(\frac{\omega }{{\omega }_c}\right)}^{s-1}{\text{e}}^{-\frac{\omega }{{\omega }_c}}\mathrm{,}\end{array}

where Γ(s) is the Euler Gamma function, and ω_c is an exponential cutoff frequency, λω_c is the coupling strength between the system and the bath. When 0 < s < 1, s = 1, and s > 1, J (ω) denotes sub-Ohmic, Ohmic, and super-Ohmic spectral densities, respectively.

Here, we take s = 1, 0.5, 3 corresponding to respectively the Ohmic, sub-Ohmic, and super-Ohmic spectral density, i.e.,

\begin{array}{}\text{(13a)}& J_{\text{Ohm}}(\omega )=\lambda \omega {\text{e}}^{-\frac{\omega }{{\omega }_c}}\mathrm{,}\end{array}

\begin{array}{}\text{(13b)}& J_{\text{sub}}(\omega )=\frac{\lambda \omega }{{\pi }^{1/2}}{\left(\frac{\omega }{{\omega }_c}\right)}^{-1/2}{\text{e}}^{-\frac{\omega }{{\omega }_c}}\mathrm{,}\end{array}

\begin{array}{}\text{(13c)}& J_{\sup }(\omega )=\frac{\lambda {\omega }^3}{2{\omega }_c^2}{\text{e}}^{-\frac{\omega }{{\omega }_c}}\mathrm{.}\end{array}

Through the analyses in the above subsection, we can see that the results of the quantum simulations are very coincident with those of the HEOM at high temperatures. In the following, we also analyze the EET dynamics for these three types of spectral densities.

Figure 3 presents the results of our quantum simulations for Ohmic, sub-Ohmic and super-Ohmic spectral densities, with different geometries, i.e., r = 13.4 Å, r = 11.3 Å, and r = 8 Å. In the numerical calculations, we take the temperature as T_EET = 3 × 10⁵ K. Interestingly, we can see that the EET dynamics is strongly dependent on the type of spectral density adopted, i.e., the statistics of the system-bath interactions. Therefore, it is crucial a simulation method can accurately simulate all types of spectral densities. Nevertheless, regardless of the form of the spectral densities, the system identically reaches the equilibrium fastest when the intra-pair distance is r = 11.3 Å, cf. Figs. 3(b1), (b2), (b3), compared with the other distances. As r is reduced, the coherent oscillations in the populations of the donors becomes more and more profound, but the EET times do not decrease monotonically. For these three spectral densities, we may arrive at the same conclusion as the Drude–Lorentzian spectral density, that the moderate-dimerized geometry explores the coherent relaxation within the donors to accelerate the energy transfer, and the over-dimerization in the geometry significantly reduces the energy transfer rate, cf. Figs. 3(c1), (c2), (c3), which is also consistent with the conclusion in Ref. [42].

Based on the above analysis, we demonstrate that this quantum simulation approach can be used to investigate the exact quantum dynamics for different Hamiltonians and various types of spectral densities. It has been proven that although theoretically the quantum dynamics with arbitrary form of spectral density can be simulated by the HEOM through spectral decomposition [27, 28], this may be practically unfeasible. As will be shown in the Section 3.4, the more complex the form of spectral density is, the higher the computational complexity of the HEOM will be. However, since the computational time of the present simulation approach is not affected by the complexity of the spectral density, we show the superiority of our method over the HEOM in the high-temperature limit.

3.3 Three-qubit system

In order to show its universality, we apply our approach to the simulation of quantum dynamics in the reaction center of photosystem II, i.e., an 8-level system. In the reaction center, both the charge separation and electron transfer take place [66, 67]. Here, a three-qubit system is employed in the quantum simulation. We use the Hamiltonian in Ref. [68] as the model to simulate the quantum dynamics, i.e.,

\begin{array}{}\text{(14)}& \frac{H_{\text{NMR}}}{2\pi }=\left(\begin{array}{cccccccc}130.32& 0.1859& -0.0652& -0.0032& -0.041& 0.025& -0.037& 0.6834\\ 0.1859& 128.64& 1.2098& 0.0072& -0.4515& -0.0282& 0.0089& -0.627\\ -0.0652& 1.2098& 128.18& 0.0072& -0.6745& 0.1553& 0.0032& -0.4749\\ -0.0032& 0.0072& 0.0072& 128.14& 0.0266& -0.0402& 0.0024& 0.00\\ -0.041& -0.4515& -0.6745& 0.0266& 127.97& 0.6302& -0.0008& 0.066\\ 0.025& -0.0282& 0.1553& -0.0402& 0.6302& 127.52& -0.0024& -0.037\\ -0.037& 0.0089& 0.0032& 0.0024& -0.0008& -0.0024& 127& 0.0233\\ 0.6834& -0.627& -0.4749& 0.00& 0.066& -0.037& 0.0233& 126.02\end{array}\right)\text{kHz}\text{.}\end{array}

In addition, we assume Drude–Lorentz spectral density as given in Eq. (6). In the simulations, we assume identical phonon relaxation rates γ_NMR = γ_EET/C = 2π × 200 Hz and identical reorganization energies λ_NMR = λ_EET/C = 2π × 0.5 Hz for all local baths. And we take the temperature as T_EET = 3 × 10⁵ K and T_NMR = 10⁻⁴ K in our numerical calculations.

In Fig. 4, we compare the simulation results by the HEOM, the random ensemble, and GRAPE algorithm. Different from Fig. 1, there are no coherent oscillations at the short-time regime due to the large detuning. The dynamics of populations monotonically arrive at the steady state after several hundreds of milliseconds. Since the results of random ensemble and GRAPE algorithm are in good agreement with those of the HEOM, we show that our approach can be effectively generalized to an N -level system with N > 4.

3.4 Computational costs of NMR simulation and HEOM

In the above sections, we show the quantum simulations of open quantum dynamics for a few-level system. Consider an N -level open quantum system, where each level of the system is coupled with an independent bath, and the correlation function of each bath contains K exponentials. The total number of auxiliary density operators for the HEOM is [43, 59]

\begin{array}{}\text{(15)}& \mathcal{I}={\displaystyle \sum _{n=0}^{\mathcal{N}}{\mathcal{I}}_n}={\displaystyle \sum _{n=0}^{\mathcal{N}}\frac{(n+KN-1)!}{n!(KN-1)!}=\frac{(\mathcal{N}+KN)!}{\mathcal{N}!(KN)!}\mathrm{,}}\end{array}

where \mathcal{N} is the hierarchy level of truncation of the HEOM. By using [69]

\begin{array}{}\text{(16)}& n!=\sqrt{2\pi }{n}^{n+\frac{1}{2}}{\text{e}}^{-n}{\text{e}}^{r_n}\mathrm{,}\end{array}

where (12n+1)⁻¹ < r_n < (12n)⁻¹, we can obtain Eq. (15) as

\begin{array}{}\text{(17)}& \begin{array}{ll}{\mathcal{I}}_n=\hfill & \sqrt{\frac{\frac{1}{\mathcal{N}}+\frac{1}{KN}}{2\pi }}{\left(1+\frac{KN}{\mathcal{N}}\right)}^{\mathcal{N}}{(1+\frac{\mathcal{N}}{KN})}^{KN}\hfill \\ \hfill & \times \exp (r_{\mathcal{N}+KN}-r_{\mathcal{N}}-r_{KN})\mathrm{.}\hfill \end{array}\end{array}

When the dimension of the system is large and the spectral density is complex, i.e., KN→∞, for a given hierarchy level of truncation, the total number of auxiliary density operators and thus the computation cost of the HEOM is approximated as

\begin{array}{}\text{(18)}& \underset{KN\to \infty }{\lim }\mathcal{I}\simeq \sqrt{\frac{1}{2\pi \mathcal{N}}}{\left(1+\frac{KN}{\mathcal{N}}\right)}^{\mathcal{N}}{\text{e}}^{KN}\mathrm{.}\end{array}

On the other hand, for a log₂ N -qubit quantum system, the GRAPE algorithm requires [56]

\begin{array}{}\text{(19)}& \begin{array}{l}(4M{{\displaystyle \log }}_{2}N+1)4^{{{\displaystyle \log }}_{2}N}=(4M{{\displaystyle \log }}_{2}N+1)N^2\hfill \end{array}\end{array}

measurements in each iteration to estimate the fitness function and its derivative with respect to the pulse amplitude, which the GRAPE aims to optimize. M is the number of the pulse sequence to be divided, which scales polynomially with the number of qubits [56]. Above all, the complexity of the GRAPE algorithm is a polynomial of N.

In practical simulations, the resources required by the HEOM may be intolerable as complicated spectral densities are widely observed in natural photosynthetic complexes. According to Ref. [70], there are 42 chlorophylls in a trimer of LHC II complex and its spectral density can be described by an overdamped Brownian oscillator and 48 high-frequency modes, i.e., N = 42 and K = 49. We remark that because of the specific mapping from N -level photosynthetic energy transfer to log₂ N -qubit quantum simulation, the complexity has been greatly reduced.

3.5 Error analysis

From the above theoretical calculations, we can see that there are some errors in the quantum simulations as compared to the theoretically-exact results predicted by the HEOM, especially when the number of random realizations in the ensemble is small. Theoretically, the average over the ensemble is equal to the average over the time only in the infinite-large ensemble limit [71]. However, according to a large number of simulations, the quantum simulations agree well with the HEOM for n≥500. For a given random Hamiltonian, the corresponding unitary evolution can be decomposed into a sequence of experimentally-feasible pulses by the GRAPE algorithm. Here, the fidelity is limited by the initial guess of parameters, in combination with both the step and the number of repetitions in attaining the global optimum [56].

In the aspect of NMR realization, there are three main sources of errors. First of all, the prepared initial state is a pseudo-pure state rather than a pure state. The error caused by the preparation of the initial state can be estimated by carrying out a full state tomography and examining the fidelity F({\rho }_{\text{T}}^{\text{P}}\mathrm{,}{\rho }_{\text{e}}^{\text{P}}) between the target density matrix {\rho }_{\text{T}}^{\text{P}} = |00〉〈00| and the experimental density matrix {\rho }_{\text{e}}^{\text{P}} [72]. Here, the fidelity is defined as F (ρ₁, ρ₂) = Tr(ρ₁ρ₂)/\sqrt{\text{Tr}({\rho }_1^2)\text{Tr}({\rho }_2^2)} . Then, in the process of quantum-state evolution, although the fidelity of the unitary evolution U_D calculated by the GRAPE algorithm can nearly approach unity, e.g., F = 0.995 by increasing the step of the repetition, there could still remain imperfection in the experimental realization of the pulse sequence. The error in this step can be also estimated by the fidelity F(U_D{\rho }_{\text{T}}^{\text{P}}{U}_D^{†}\mathrm{,}{\rho }_{\text{e}}^{\text{G}}) with {\rho }_{\text{e}}^{\text{G}} obtained by the tomography. Finally, further errors could be introduced in the process of measurement. The error produced in this step can be characterized by the fidelity between the final state obtained from the experiment and the final state of the theory.

4 Conclusion

In this paper, we discuss the recently-developed approach for the exact simulation of EET dynamics in photosynthesis. By applying the approach to the linear-tetramer model, the energy transfer is shown to be optimal for a moderately-clustered geometry. Based on the optimal geomtry, we show that the energy transfer efficiency can be further improved when the energy gap between the donor and acceptor clusters matches the optimum of the bath’s spectral density. In this regard, we demonstrate that the light-harvesting network can be optimized from two aspects, i.e., the geometry [42] and the bath [49].

Beyond the Drude–Lorentz spectral density, we also show that our approach can be utilized to simulate the EET dynamics for various types of spectral densities. We furthermore demonstrate our approach to the simulation of a quantum system with N = 8, i.e., the reaction center of PSII. Therefore, we can generalize our approach to exactly simulate the open quantum dynamics with N dimensions and different types of environments. By comparing our approach to the HEOM, the complexity is a polynomial of the number of states involved and thus our approach is exponentially accelerated when the system is large and the spectral density is complicated. This is often encountered when simulating open quantum dynamics for natural photosynthesis. To conclude, our approach can be applied to the exact and efficient simulation of open quantum dynamics for various of Hamiltonians and spectral densities. It may shed light on the investigations exploring non-Markovianity in open quantum dynamics, e.g., quantum metrology in non-Markovian environments [65].

References

[1] G. R. Fleming and R. Grondelle, The primary steps of photosynthesis, Phys. Today 47(2), 48 (1994)

[2] Y. C. Cheng and G. R. Fleming, Dynamics of light harvesting in photosynthesis, Annu. Rev. Phys. Chem. 60(1), 241 (2009)

[3] M. J. Tao, N. N. Zhang, P. Y. Wen, F. G. Deng, Q. Ai, and G. L. Long, Coherent and incoherent theories for photosynthetic energy transfer, Sci. Bull. (Beijing) 65(4), 318 (2020)

[4] M. J. Tao, M. Hua, N. N. Zhang, W. T. He, Q. Ai, and F. G. Deng, Quantum simulation of clustered photosynthetic light harvesting in a superconducting quantum circuit, Quantum Eng. 2(3), e53 (2020)

[5] N. Lambert, Y. N. Chen, Y. C. Cheng, C. M. Li, G. Y. Chen, and F. Nori, Quantum biology, Nat. Phys. 9(1), 10 (2013)

[6] J. S. Cao, R. J. Cogdell, D. F. Coker, H. G. Duan, J. Hauer, U. Kleinekathöfer, T. L. C. Jansen, T. Mančal, R. J. D. Miller, J. P. Ogilvie, V. I. Prokhorenko, T. Renger, H. S. Tan, R. Tempelaar, M. Thorwart, E. Thyrhaug, S. Westenhoff, and D. Zigmantas, Quantum biology revisited, Sci. Adv. 6(14), eaaz4888 (2020)

[7] G. S. Engel, T. R. Calhoun, E. L. Read, T. K. Ahn, T. Mančal, Y. C. Cheng, R. E. Blankenship, and G. R. Fleming, Evidence for wavelike energy transfer through quantum coherence in photosynthetic systems, Nature 446(7137), 782 (2007)

[8] H. Lee, Y. C. Cheng, and G. R. Fleming, Coherence dynamics in photosynthesis: Protein protection of excitonic coherence, Science 316(5830), 1462 (2007)

[9] P. G. Wolynes, Some quantum weirdness in physiology, Proc. Natl. Acad. Sci. USA 106(41), 17247 (2009)

[10] E. Collini, C. Y. Wong, K. E. Wilk, P. M. G. Curmi, P. Brumer, and G. D. Scholes, Coherently wired lightharvesting in photosynthetic marine algae at ambient temperature, Nature 463(7281), 644 (2010)

[11] R. Hildner, D. Brinks, J. B. Nieder, R. J. Cogdell, and N. F. van Hulst, Quantum coherent energy transfer over varying pathways in single light-harvesting complexes, Science 340(6139), 1448 (2013)

[12] M. J. Tao, Q. Ai, F. G. Deng, and Y. C. Cheng, Proposal for probing energy transfer pathway by single-molecule pump-dump experiment, Sci. Rep. 6(1), 27535 (2016)

[13] L. G. Mourokh and F. Nori, Energy transfer efficiency in the chromophore network strongly coupled to a vibrational mode, Phys. Rev. E 92(5), 052720 (2015)

[14] H. P. Breuer, E. M. Laine, J. Piilo, and B. Vacchini, Non- Markovian dynamics in open quantum systems, Rev. Mod. Phys. 88(2), 021002 (2016)

[15] I. de Vega and D. Alonso, Dynamics of non-Markovian open quantum systems, Rev. Mod. Phys. 89(1), 015001 (2017)

[16] L. Li, M. J. W. Hall, and H. M. Wiseman, Concepts of quantum non-Markovianity: A hierarchy, Phys. Rep. 759, 1 (2018)

[17] H. P. Breuer and F. Petruccione, The Theory of Open Quantum Systems, Oxford University Press, New York, 2007

[18] A. Ishizaki and G. R. Fleming, On the adequacy of the Redfield equation and related approaches to the study of quantum dynamics in electronic energy transfer, J. Chem. Phys. 130(23), 234110 (2009)

[19] G. Watanabe, Heat engines using small quantum systems, AAPPS Bull. 29, 30 (2019)

[20] J. X. Zhao, J. J. Cheng, Y. Q. Chu, Y. X. Wang, F. G. Deng, and Q. Ai, Hyperbolic metamaterial using chiral molecules, Sci. China Phys. Mech. Astron. 63(6), 260311 (2020)

[21] Y. Tanimura, Stochastic Liouville, Langevin, Fokker- Planck, and master equation approaches to quantum dissipative systems, J. Phys. Soc. Jpn. 75(8), 082001 (2006)

[22] A. Ishizaki and G. R. Fleming, Unified treatment of quantum coherent and incoherent hopping dynamics in electronic energy transfer: Reduced hierarchy equation approach, J. Chem. Phys. 130(23), 234111 (2009)

[23] Y. Yan, F. Yan, Y. Liu, and J. Shao, Hierarchical approach based on stochastic decoupling to dissipative systems, Chem. Phys. Lett. 395(4–6), 216 (2004)

[24] Y. Zhou, Y. Yan, and J. Shao, Stochastic simulation of quantum dissipative dynamics, Europhys. Lett. 72(3), 334 (2005)

[25] J. Shao, Decoupling quantum dissipation interaction via stochastic fields, J. Chem. Phys. 120(11), 5053 (2004)

[26] Z. F. Tang, X. L. Ouyang, Z. H. Gong, H. B. Wang, and J. L. Wu, Extended hierarchy equation of motion for the spin-boson model, J. Chem. Phys. 143(22), 224112 (2015)

[27] H. Liu, L. L. Zhu, S. M. Bai, and Q. Shi, Reduced quantum dynamics with arbitrary bath spectral densities: Hierarchical equations of motion based on several different bath decomposition schemes, J. Chem. Phys. 140(13), 134106 (2014)

[28] M. Schröder, M. Schreiber, and U. Kleinekathöfer, Reduced dynamics of coupled harmonic and anharmonic oscillators using higherorder perturbation theory, J. Chem. Phys. 126(11), 114102 (2007)

[29] A. Olaya-Castro, C. F. Lee, F. F. Olsen, and N. F. Johnson, Efficiency of energy transfer in a light-harvesting system under quantum coherence, Phys. Rev. B 78(8), 085115 (2008)

[30] Q. Ai, Y. J. Fan, B. Y. Jin, and Y. C. Cheng, An efficient quantum jump method for coherent energy transfer dynamics in photosynthetic systems under the influence of laser fields, New J. Phys. 16(5), 053033 (2014)

[31] S. Jang, Y. C. Cheng, D. R. Reichman, and J. D. Eaves, Theory of coherent resonance energy transfer, J. Chem. Phys. 129(10), 101104 (2008)

[32] M. Yang and G. R. Fleming, Influence of phonons on exciton transfer dynamics: Comparison of the Redfield, F rster, and modified Redfield equations, Chem. Phys. 282(1), 163 (2002)

[33] Y. H. Hwang-Fu, W. Chen, and Y. C. Cheng, A coherent modified Redfield theory for excitation energy transfer in molecular aggregates, Chem. Phys. 447, 46 (2015)

[34] H. Dong, D. Z. Xu, J. F. Huang, and C. P. Sun, Coherent excitation transfer via the dark-state channel in a bionic system, Light Sci. Appl. 1(3), e2 (2012)

[35] S. Mostarda, F. Levi, D. Prada-Gracia, F. Mintert, and F. Rao, Structure-dynamics relationship in coherent transport through disordered systems, Nat. Commun. 4(1), 2296 (2013)

[36] G. C. Knee, P. Rowe, L. D. Smith, A. Troisi, and A. Datta, Structure-dynamics relation in physically-plausible multi-chromophore systems, J. Phys. Chem. Lett. 8(10), 2328 (2017)

[37] T. Zech, R. Mulet, T. Wellens, and A. Buchleitner, Centrosymmetry enhances quantum transport in disordered molecular networks, New J. Phys. 16(5), 055002 (2014)

[38] L. Xu, Z. R. Gong, M. J. Tao, and Q. Ai, Artificial light harvesting by dimerized Möbius ring, Phys. Rev. E 97(4), 042124 (2018)

[39] Y. H. Lui, B. Zhang, and S. Hu, Rational design of photoelectrodes for photoelectrochemical water splitting and CO₂ reduction, Front. Phys. 14(5), 53402 (2019)

[40] L. Ju, M. Bie, X. Zhang, X. Chen, and L. Kou, Two-dimensional Janus van der Waals heterojunctions: A review of recent research progresses, Front. Phys. 16(1), 13201 (2021)

[41] B. X. Wang, M. J. Tao, Q. Ai, T. Xin, N. Lambert, D. Ruan, Y. C. Cheng, F. Nori, F. G. Deng, and G. L. Long, Efficient quantum simulation of photosynthetic light harvesting, npj Quantum Inf. 4, 52 (2018)

[42] Q. Ai, T. C. Yen, B. Y. Jin, and Y. C. Cheng, Clustered geometries exploiting quantum coherence effects for efficient energy transfer in light harvesting, J. Phys. Chem. Lett. 4(15), 2577 (2013)

[43] Q. Shi, L. Chen, G. Nan, R. X. Xu, and Y. J. Yan, Efficient hierarchical liouville space propagetor to quantum dissipative dynamics, J. Chem. Phys. 130(8), 084105 (2009)

[44] I. Buluta and F. Nori, Quantum simulators, Science 326(5949), 108 (2009)

[45] I. M. Georgescu, S. Ashhab, and F. Nori, Quantum simulation, Rev. Mod. Phys. 86(1), 153 (2014)

[46] J. Xu, S. Li, T. Chen, and Z. Y. Xue, Nonadiabatic geometric quantum computation with optimal control on superconducting circuits, Front. Phys. 15(4), 41503 (2020)

[47] Z. D. Ye, D. Pan, Z. Sun, C. G. Du, L. G. Yin, and G. L. Long, Generic security analysis framework for quantum secure direct communication, Front. Phys. 16(2), 21503 (2021)

[48] Y. F. Yan, L. Zhou, W. Zhong, and Y. B. Sheng, Measurementdevice-independent quantum key distribution of multiple degrees of freedom of a single photon, Front. Phys. 16(1), 11501 (2021)

[49] M. Rey, A. W. Chin, S. F. Huelga, and M. B. Plenio, Exploiting structured environments for efficient energy transfer: The phonon antenna mechanism, J. Phys. Chem. Lett. 4(6), 903 (2013)

[50] D. J. Gorman, B. Hemmerling, E. Megidish, S. A. Moeller, P. Schindler, M. Sarovar, and H. Haeffner, Engineering vibrationally assisted energy transfer in a trapped-ion quantum simulator, Phys. Rev. X 8(1), 011038 (2018)

[51] Y. Chang and Y. C. Cheng, On the accuracy of coherent modified Redfield theory in simulating excitation energy transfer dynamics, J. Chem. Phys. 142(3), 034109 (2015)

[52] C. Meier and D. J. Tannor, Non-Markovian evolution of the density operator in the presence of strong laser fields, J. Chem. Phys. 111(8), 3365 (1999)

[53] A. Soare, H. Ball, D. Hayes, J. Sastrawan, M. C. Jarratt, J. J. McLoughlin, X. Zhen, T. J. Green, and M. J. Biercuk, Experimental noise filtering by quantum control, Nat. Phys. 10(11), 825 (2014)

[54] A. Soare, H. Ball, D. Hayes, X. Zhen, M. C. Jarratt, J. Sastrawan, H. Uys, and M. J. Biercuk, Experimental bath engineering for quantitative studies of quantum control, Phys. Rev. A 89(4), 042329 (2014)

[55] N. Khaneja, T. Reiss, C. Kehlet, T. Schulte-Herbrüggen, and S. J. Glaser, Optimal control of coupled spin dynamics: Design of NMR pulse sequences by gradient ascent algorithms, J. Magn. Reson. 172(2), 296 (2005)

[56] J. Li, X. D. Yang, X. H. Peng, and C. P. Sun, Hybrid quantum-classical approach to quantum optimal control, Phys. Rev. Lett. 118(15), 150503 (2017)

[57] P. Fulde, Wavefunctions for extended electron systems, AAPPS Bull. 29, 50 (2019)

[58] L. Valkunas, D. Abramavicius, and T. Mančal, Molecular Excitation Dynamics and Relaxation: Quantum Theory and Spectroscopy, Wiley-VCH, Weinheim, Germany, 2013

[59] A. Ishizaki, and G. R. Fleming, Theoretical examination of quantum coherence in a photosythetic system at physiological temperature, Proc. Natl. Acad. Sci. USA 106(41), 17255 (2009)

[60] W. Jiang, F. Z. Wu, and G. J. Yang, Non-Markovian entanglement dynamics of open quantum systems with continuous measurement feedback, Phys. Rev. A 98(5), 052134 (2018)

[61] X. L. Zhen, F. H. Zhang, G. Y. Feng, L. Hang, and G. L. Long, Optimal experimental dynamical decoupling of both longitudinal and transverse relaxations, Phys. Rev. A 93(2), 022304 (2016)

[62] Y. H. Ma, H. Dong, H. T. Quan, and C. P. Sun, The uniqueness of the integration factor associated with the exchanged heat in thermodynamics, Fundamental Research 1(1), 6 (2021)

[63] A. J. Leggett, S. Chakravarty, A. Dorsey, M. Fisher, A. Garg, and W. Zwerger, Dynamics of the dissipative twostate system, Rev. Mod. Phys. 59(1), 1 (1987)

[64] U. Weiss, Quantum Dissipative Systems, World Scientific, Singapore, 2008

[65] A. W. Chin, S. F. Huelga, and M. B. Plenio, Quantum metrology in non-Markovian environments, Phys. Rev. Lett. 109(23), 233601 (2012)

[66] H. G. Duan, V. I. Prokhorenko, E. Wientjes, R. Croce, M. Thorwart, and R. J. D. Miller, Primary charge separation in the photosystem II reaction center revealed by a global analysis of the two-dimensional electronic spectra, Sci. Rep. 7(1), 12347 (2017)

[67] K. L. M. Lewis, F. D. Fuller, J. A. Myers, C. F. Yocum, D. Abramavicius, and J. P. Ogilvie, Simulations of the two-dimensional electronic spectroscopy of the photosystem II reaction center, J. Phys. Chem. A 117(1), 34 (2013)

[68] L. Zhang, D. A. Silva, H. D. Zhang, A. Yue, Y. J. Yan, and X. H. Huang, Dynamic protein conformations preferentially drive energy transfer along the active chain of the photosystem II reaction centre, Nat. Commun. 5(1), 4170 (2014)

[69] H. Robbins, A remark on Stirling’s formula, Am. Math. Mon. 62, 26 (1955)

[70] V. I. Novoderezhkin, M. A. Palacios, H. van Amerongen, and R. van Grondelle, Energy-transfer dynamics in the LHCII complex of higher plants: Modified Redfield approach, J. Phys. Chem. B 108(29), 10363 (2004)

[71] J. W. Goodman, Statistical Optics, 2nd Ed., Wiley, Hoboken, NJ, 2015

[72] D. W. Lu, N. Y. Xu, R. X. Xu, H. W. Chen, J. B. Gong, X. H. Peng, and J. F. Du, Simulation of chemical isomerization reaction dynamics on a NMR quantum simulator, Phys. Rev. Lett. 107(2), 020501 (2011)

[73] I. L. Chuang, L. M. K. Vandersypen, X. L. Zhou, D. W. Leung, and S. Lloyd, Experimental realization of a quantum algorithm, Nature 393(6681), 143 (1998)

[74] L. M. K. Vandersypen and I. Chuang, NMR techniques for quantum control and computation, Rev. Mod. Phys. 76(4), 1037 (2005)

[75] E. Knill, I. Chuang, and R. Laflamme, Effective pure states for bulk quantum computation, Phys. Rev. A 57(5), 3348 (1998)

[76] D. G. Cory, M. D. Price, and T. F. Havel, Nuclear magnetic resonance spectroscopy: An experimentally accessible paradigm for quantum computing, Physica D 120(1–2), 82 (1998)

[77] J. S. Lee, The quantum state tomography on an NMR system, Phys. Lett. A 305(6), 349 (2002)

[78] D. W. Lu, T. Xin, N. K. Yu, Z. F. Ji, J. X. Chen, G. L. Long, J. Baugh, X. H. Peng, B. Zeng, and R. Laflamme, Tomography is necessary for universal entanglement detection with single-copy observables, Phys. Rev. Lett. 116(23), 230501 (2016)

[79] T. Xin, D. W. Lu, J. Klassen, N. K. Yu, Z. F. Ji, J. X. Chen, X. Ma, G. L. Long, B. Zeng, and R. Laflamme, Quantum state tomography via reduced density matrices, Phys. Rev. Lett. 118(2), 020401 (2017)

Rights & permissions

Higher Education Press 2021