Exact and efficient quantum simulation of open quantum dynamics
for various of Hamiltonians and spectral densities

Following Ref. Wang18 , 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. We adopt the Frenkel-exciton Hamiltonian Cheng09 ; Valkunas13

to describe the dynamics of EET in photosynthesis. Here $|i\rangle~{}(i=1,2,3,4)$ represents that only $i$th molecule is in the excited state while the others are in the ground state. $\varepsilon_{i}$ is the site energy of $i$th exciton. For simplicity, the electronic interaction between $i$th and $j$th excitons is given by the dipole-dipole interaction as

where $\vec{r}_{ij}=r_{ij}\hat{r}_{ij}$ is the displacement vector from site $i$ to site $j$, $\vec{\mu}_{i}$ the transition dipole of site $i$, $\varepsilon_{0}$ the vacuum permittivity. In numerical simulations, we take $\mu_{j}=7.75~{}\textrm{D}$ and $r\in[6,14]~{}\rm{\mathring{A}}$, which are typical in natural photosynthesis.

For photosynthetic complexes with $N$ chlorophylls, there are only $N$ single-excitation states in the process of energy transfer. As a result, only $\log_{2}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\rangle=|00\rangle$, $|2\rangle=|01\rangle$, $|3\rangle=|10\rangle$, and $|4\rangle=|11\rangle$. As shown in Fig. 1(b), we regard H and C nuclear spins as the two qubits. The Hamiltonian $H_{\rm{NMR}}$ implemented in NMR simulations is $H_{\rm{NMR}}=H_{\rm{EET}}/C$ with a scaling factor

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

Here, $V_{j}=|i\rangle\langle i|$, $a_{ik}^{\dagger}$ ($a_{ik}$) is the creation (annihilation) operator of the $k$th phonon mode of the $i$th molecule, and $g_{ik}$ represents the coupling strength. $H_{\rm{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.,

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

The exact EET dynamics can be given by the HEOM Ishizaki09-2 ; Ishizaki09-3

where $\sigma_{\vec{n}}$ and $\sigma_{\vec{n}_{j\pm}}$ are the auxiliary density matrices with $\vec{n}=(n_{1},n_{2},\cdots,n_{j},\cdots)$ and $\vec{n}_{j\pm}=(n_{1},n_{2},\cdots,n_{j}\pm 1,\cdots)$, and $n_{j}$’s non-negative integers, $\sigma_{\vec{0}}=\rho$ the reduced density matrix of photosynthetic system, $\mathcal{L}_{e}$ the Liouville superoperator of $H_{\textrm{EET}}$. Besides, $\Phi_{j}=iV_{j}^{\times}$ and $\Theta_{j}=i(\frac{2\lambda_{j}}{\beta\hbar^{2}}V_{j}^{\times}-i\frac{\lambda_{j}}{\hbar}V_{j}^{\circ})$ with $V_{j}^{\times}\sigma_{\vec{n}}=V_{j}\sigma_{\vec{n}}-\sigma_{\vec{n}}V_{j}$ and $V_{j}^{\circ}\sigma_{\vec{n}}=V_{j}\sigma_{\vec{n}}+\sigma_{\vec{n}}V_{j}$, where $V_{j}=|j\rangle\langle j|$. For Drude-Lorentz spectral density, i.e.,

with $\lambda_{j}$ the reorganization energy and $\gamma_{j}$ the relaxation rate. For a generic spectral density, it can be decomposed into a summation of Lorentzian form Meier99 ; Jiang18 .

In the quantum simulation of photosynthetic EET with NMR system, there are three steps: First of all, we initialize the system for quantum simulation, i.e., preparation of pseudo-pure states. Afterwards, the system evolves under the evolution operator $U$. Finally, the state of the system is measured. In the following, we shall discuss the details of the three steps.

In Ref. Ai13 , 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. 2(a), we show the comparison among the HEOM, the quantum simulation and the GRAPE simulation for $r=13.4~{}\rm{\mathring{A}}$. Therein, all four chlorophylls are equally spaced and the Hamiltonian for the NMR quantum simulation $H_{\rm{NMR}}=H_{\rm{EET}}/C$ is

In Fig. 2(b), we set $r=11.3~{}\rm{\mathring{A}}$ 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

In Fig. 2(c), we further reduce the intra-dimer distance to an even smaller value, i.e., $r=8~{}\rm{\mathring{A}}$. Thus, the Hamiltonian is

In order to obtain the above Hamiltonians, we adopt the site energies $\varepsilon_{1}\!=\!13000~{}\rm{cm}^{-1}$, $\varepsilon_{2}\!=\!12900~{}\rm{cm}^{-1}$, $\varepsilon_{3}\!=\!12300~{}\rm{cm}^{-1}$, and $\varepsilon_{4}\!=\!12200~{}\rm{cm}^{-1}$, which are typical in natural photosynthetic systems, for the corresponding Hamiltonians of photosynthesis. Except that $J_{14}=J_{41}$ are 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-Lorentzian form as given in Eq. (7), with the optimal frequency $\gamma_{j}$. In the simulations, we assume identical phonon relaxation rates $\gamma_{\textrm{NMR}}=\gamma_{\textrm{EET}}/C=2\pi\times 50~{}\rm{Hz}$ and identical reorganization energies $\lambda_{\textrm{NMR}}=\lambda_{\textrm{EET}}/C=2\pi\times 2~{}\rm{Hz}$ for all local baths. Since the results of quantum simulation coincide with those of the HEOM only at high temperatures, we take $T_{\rm{EET}}=3\times 10^{6}$ K and $T_{\rm{NMR}}=10^{-3}~{}\rm{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. 2.

As shown in Fig. 2(b), the energy transfer is the fastest for the case with $r=11.3~{}\rm{\mathring{A}}$. 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. 2(c). These numerically-exact simulations are consistent with those discoveries obtained by the approximate theory in Ref. Ai13 . 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. 2(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, of which causes 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. Rey13 , it was shown that the optimization of the EET can be achieved when the energy gap of the system matches the optimum in 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~{}\rm{\mathring{A}}$, which was shown to be optimal among different geometries in Fig. 2. Since the peak of the Drude-Lorentz spectral density is located at $\gamma_{\textrm{NMR}}$, we fix $\lambda_{\textrm{NMR}}$ and simulate the EET dynamics for a broad range of $\gamma_{\textrm{NMR}}$. Here, in Fig. 3, we only demonstrate the population dynamics for three typical parameters, i.e., $\gamma_{\textrm{NMR}}/2\pi=0.5~{}\rm{kHz}$, $2.668~{}\rm{kHz}$, and $7~{}\rm{kHz}$. Obviously, the optimal relaxation rate of the bath is $\gamma_{\textrm{NMR}}^{\textrm{opt}}=2\pi\times 2.668~{}\rm{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 B.

Since a general power-law form spectral density can describe an extremely-large number of physical environments Weiss08 ; Leggett87 , we use it to show the universal applicability of the quantum simulation approach. Following Refs. Breuer07 ; Chin12 , 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

where $\Gamma(s)$ is the Euler Gamma function, and $\omega_{c}$ is an exponential cutoff frequency, $\lambda\omega_{c}$ is the coupling strength between the system and the bath. When $0<s<1$, $s=1$, and $s>1$, $J(\omega)$ 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.,

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 4 presents the results of our quantum simulations for Ohmic, sub-Ohmic and super-Ohmic spectral densities, with different geometries, i.e., $r=13.4~{}\rm{\mathring{A}}$, $r=11.3~{}\rm{\mathring{A}}$, and $r=8~{}\rm{\mathring{A}}$. In the numerical calculations, we take the temperature as $T_{\textrm{EET}}=3\times 10^{5}~{}\rm{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~{}\rm{\mathring{A}}$, cf. Figs. 4(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. 4(c1),(c2),(c3), which is also consistent with the conclusion in Ref. Ai13 .

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 Liu14 ; Schroder07 , this may be practically unfeasible. As will be shown in the next subsection, 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.

In the above sections, we show the quantum simulations of open quantum dynamics for a 4-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 Shi09 ; Ishizaki09-3

where $\mathcal{N}$ is the hierarchy level of truncation of the HEOM. By using Robbins1955

where $(12n+1)^{-1}<r_{n}<(12n)^{-1}$, we can obtain Eq. (33) as

When the dimension of the system is large and the spectral density is complex, i.e., $KN\rightarrow\infty$, 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

On the other hand, for a $\log_{2}{N}$-qubit quantum system, the GRAPE algorithm requires Li17

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 Li17 . 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. Novoderezhkin04 , 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_{2}{N}$-qubit quantum simulation, the complexity has been greatly reduced.

From the above theoretical calculations, we can see that there are some errors in the quantum simulations as compared to the theoretically-exact results 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 Goodman15 . However, according to a large number of simulations, the quantum simulations agree well with the HEOM for $n\geq 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 Li17 .

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. Then, in the process of quantum-state evolution, although the fidelity of the unitary evolution $U_{D}$ calculated by the GRAPE algorithm can in principle approach unity, e.g. by increasing the step of the repetition, there could still remain imperfection in the experimental realization of the pulse sequence. Finally, further errors could be introduced in the process of measurement.