# Quantum simulation of photosynthetic energy transfer

^{1}

^{2}

^{3}

Near-unity energy transfer efficiency has been widely observed in natural photosynthetic complexes. This phenomenon has attracted broad interest from different fields, such as physics, biology, chemistry and material science, as it may offer valuable insights into efficient solar-energy harvesting. Recently, quantum coherent effects have been discovered in photosynthetic light harvesting, and their potential role on energy transfer has seen heated debate. Here, we perform an experimental quantum simulation of photosynthetic energy transfer using nuclear magnetic resonance (NMR). We show that an -chromophore photosynthetic complex, with arbitrary structure and bath spectral density, can be effectively simulated by a system with qubits. The computational cost of simulating such a system with a theoretical tool, like the hierarchical equation of motion, which is exponential in , can be potentially reduced to requiring a just polynomial number of qubits using NMR quantum simulation. The benefits of performing such quantum simulation in NMR are even greater when the spectral density is complex, as in natural photosynthetic complexes. These findings may shed light on quantum coherence in energy transfer and help to provide design principles for efficient artificial light harvesting.

Efficient exciton energy transfer (EET) is crucial in photosynthesis and solar cells ^{31, 2},
especially when the systems are large ^{3},
e.g., as in Photosystem I (PSI) and Photosystem II (PSII),
which have hundreds of chromophores.
A comprehensive knowledge of the quantum dynamics of such systems
would be of potential importance to the study on EET ^{32}.
Much effort has been made to reveal the effects of quantum coherence
on efficient energy transfer ^{32, 5, 6, 7, 8, 9, 49}.
In order to try to mimic EET,
a Frenkel-exciton Hamiltonian is required and
this can be studied with quantum chemistry approaches ^{31},
e.g., fitting experimental spectra, or calculations by density functional theory.
Because photosynthetic pigment-protein complexes are intrinsically open quantum systems,
with system-bath couplings comparable to the intra-system couplings,
it is difficult to faithfully mimic the exact quantum dynamics of EET.
Among the methods for describing EET ^{31, 49, 12},
the hierarchical equation of motion (HEOM) yields a numerically exact solution
at the cost of considerable computing time ^{37, 49}. For the widely-studied Fenna-Matthews-Olson (FMO) complex, a reliable result could be produced by 16,170 coupled differential equations at low temperatures, based on the HEOM ^{37, 49} with 4 layers. Moreover, the HEOM encounters difficulties when the system
has non-trivial spectral density,
making it often difficult to verify the results.
On the other hand, due to heterogeneity, e.g. static disorder and conformational change, it is
difficult to experimentally verify the theoretical predictions in natural
photosynthetic systems. Because quantum simulations with qubits can powerfully mimic
the quantum dynamics of states by virtue of quantum mechanics ^{13, 14, 15},
the quantum simulation of photosynthetic energy transfer with an arbitrary Hamiltonian
by nuclear magnetic resonance (NMR) provides an intriguing approach to verify theoretical predictions. Recently, a newly-developed technique which could simulate
the effect of a bath with an arbitrary spectral density
by a set of classical pulses has been successfully realized with ion traps ^{34} and NMR ^{35}. Therefore, photosynthetic light-harvesting with an arbitrary Hamiltonian
and spectral density,
describing a structured environment, can be experimentally simulated by NMR.

NMR is an excellent platform for quantum simulation since it is easy
to operate and it can have long coherence times ^{18, 19}.
In this paper, NMR is utilized to simulate the quantum coherent dynamics in photosynthetic
light harvesting. As a prototype, a tetramer including four chlorophylls ^{12},
as schematically illustrated in Fig. 1(a), is
employed in the NMR quantum simulation. In a previous investigation ^{12},
the tetramer model was exploited to study the clustered geometry utilizing
exciton delocalization and energy matching to accelerate the energy transfer.
The EET in photosynthesis is described
by the Frenkel-exciton Hamiltonian

(1) |

where () is the state with a single excitation at site
and other sites at the ground state, is the site energy of ,
is the excitonic interaction between sites and .
In our quantum simulations, we adopt the site energies ^{12} cm, cm, cm, and cm, and the couplings cm, cm, cm, and cm, which are typical parameters in
photosynthetic systems.

For a photosynthetic complex with chlorophylls, there are only single-excitation states involved in the energy transfer. Therefore, only qubits are required to realize the quantum simulation. To mimic the energy transfer described by the above Hamiltonian (1), two qubits are necessary for the quantum simulation. In this case, the photosynthetic single-excitation state is encoded as a two-qubit product state, i.e. , , , and . By a straightforward calculation, the Frenkel-exciton Hamiltonian can be mapped to the NMR Hamiltonian as

(2) | |||||

where (, ) is the Pauli operator for qubit , numerically and , but the dimension cm should be replaced by kHz. In other words, all realistic parameters have been scaled down in energy by a factor of .

The excitonic coupling between sites and
makes the exciton energy hop in both directions,
i.e. .
Apart from this, the system-bath couplings will
facilitate the energy flow towards an energy trap,
where the captured photon energy is converted into
chemical energy ^{2, 3}.
The interaction between the system and bath in photosynthetic complexes
can be described by

(3) |

where is the creation operator for the th phonon mode of chlorophyll with coupling strength . will induce pure-dephasing on the -th chromophore when it is in the excited state. Generally, the system-bath couplings are given by the spectral density

(4) |

which we assume identical for all chromophores. For typical photosynthetic complexes, the system-bath couplings are of the same order as the intra-system couplings. In order to mimic the effects of noise, we utilize the bath-engineering technique which has been successfully implemented in ion traps and NMR. The implementation of a pure-dephasing Hamiltonian

(5) |

relies on generating stochastic errors by performing phase modulations on a constant-amplitude carrier, i.e.

(6) | |||||

(7) |

where is the constant amplitude of a magnetic field with driving frequency , is a random number, with and being base and cutoff frequencies respectively, and is a global scaling factor. The power-spectral density is the Fourier transform of the autocorrelation function of ,

(8) |

By appropriately choosing and the cutoff , an arbitrary power-spectral density of the bath can be realized, e.g. white, Ohmic, or Debye spectral densities. For the details, please refer to the Supplementary Information.

As illustrated in Fig. 1(b), the nuclear spins of the carbon atom and hydrogen atom in a chloroform molecule are chosen to encode the two qubits with the Hamiltonian written as

(9) |

where Hz, Hz are the chemical shifts of the two spins,
and Hz ^{20}.
Therefore, in order to simulate the quantum dynamics of energy
transfer, the entire evolution is decomposed into repetitive identical cycles.
After each cycle, there is a small difference between the exact evolution
and the simulated one. However, because there are tens of thousands of cycles before
completing the energy transfer, the accumulated error might be significant so
that the simulation becomes unreliable. To avoid this problem,
we utilize the gradient ascent pulse engineering (GRAPE) algorithm ^{44}, which has been successfully applied to a number of quantum simulations in NMR ^{22, 19, 23},
to mimic the quantum dynamics of light harvesting.

In Fig. 2(a), the quantum coherent energy transfer for the above Hamiltonian , using an initial condition where an excitation is localized on the first chromophore, is simulated in NMR. In the short-time regime, cf. Fig. 2(b), there are Rabi-like oscillations of coherent energy transfer between the two levels with the highest energies, because there is a strong coupling between these two levels. Furthermore, the oscillation quickly damps as the energy transfer is irreversible due to the pure-dephasing noise. After an exponential decay process, the populations on all levels reach thermal equilibrium. Noticeably, there are small oscillations for the two lowest-energy levels as a result of their strong coupling. For each point in the quantum simulation, the data is averaged over random realizations. For a given realization, the system undergoes a coherent evolution by applying a time-dependent magnetic field with fluctuating phases, as shown in Eqs. (6,7). However, for an ensemble of random realizations, since each realization experiences a different phase at a given time, the ensemble average manifests itself as a single realization undergoing a pure-dephasing noise. In this regard, the deviation of the NMR simulation from that predicted by the HEOM decreases as the number of realizations in the ensemble increases, cf. the Supplementary Information. This effect would be more remarkable if were increased further, as confirmed by our theoretical simulations in the Supplementary Information. In conclusion, the dephasing Hamiltonian (3) in photosynthesis is effectively mimicked by a classical time-fluctuating magnetic field (6,7).

Note that the proof-of-concept NMR experiment presented in this work
goes beyond reproducing the HEOM results. To faithfully simulate
the EET dynamics on a large photosynthetic system, using the HEOM would
be computationally unaffordable. In addition, the HEOM ^{37, 49}
is known to encounter difficulties when the spectral density is not in a simple
Drude-Lorentz form.
As shown in the Supplementary Information, the computational cost of HEOM scales
exponentially with the system size and the number of exponentials in the bath
correlation function ^{37, 49}.
Our NMR simulations do not have such shortcomings, as it scales polynomially
with respect to the system size ^{51, 52}.
Moreover, the quantum simulation algorithm also scales more favorably
in terms of the number of exponentials in the bath correlation function.
With our current approach, in principle, a photosynthetic system with sites (e.g. the PSI complex)
requires a 7-qubit quantum simulator. This means that the coherent EET
dynamics of a full functional biological unit for photosynthesis
(e.g. the PSII supercomplex that contains sites) can be
faithfully simulated by a 9-qubit NMR quantum computer,
which is clearly attainable by the NMR technique ^{52}. Thus,
we believe that NMR quantum simulation has the potential to
help clarify the mysteries of light harvesting in natural photosynthesis.
In this regard, other approaches, which do not take advantage of this scaling provided by encoding multiple sites into a single qubit ^{26},
currently lack the required size to simulate such large photosynthetic systems, and thus may also benefit from the approach we develop here.

In this paper, the photosynthetic energy transfer is experimentally
simulated in NMR. As a prototype, a two-qubit NMR system is utilized
to demonstrate both the coherent oscillations at the short-time regime
and the steady-state thermalization at the long-time.
By using the GRAPE technique, an arbitrary photosynthetic system can
be faithfully mapped to the NMR system. Besides,
the effect on EET can be effectively
mimicked by a set of well-designed pulses, which act as a classical pure-dephasing noise.
The quantum simulation of photosynthetic energy transfer in NMR would probably facilitate the investigation of quantum coherent effects on the EET
and more clear design principles for artificial
light-harvesting devices together with structured baths ^{12, 27}.

A recent work ^{26} experimentally verified that a structured bath can optimize the energy transfer in EET ^{27}. However, here we showed we can use NMR and bath-engineering
techniques to efficiently mimic photosynthetic light harvesting in large systems with arbitrary system structure and bath spectral density. Future extensions of our approach include encoding quantum properties of the environment in ancillary qubits.

Acknowledgements We thank Xin-Long Zhen, and Ke-Ren Li for discussions. GLL was supported by National Natural Science Foundation of China under Grant Nos. 11175094 and 91221205, and the National Basic Research Program of China under Grant No. 2015CB921001. FGD was supported by the National Natural Science Foundation of China under Grant No. 11474026 and No. 11674033. FN was supported by the MURI Center for Dynamic Magneto-Optics via the AFOSR Award No. FA9550-14-1-0040, the Japan Society for the Promotion of Science (KAKENHI), the IMPACT program of JST, JSPS-RFBR grant No. 17-52-50023, CREST grant No. JPMJCR1676, and RIKEN-AIST Challenge Research Fund. QA was supported by the National Natural Science Foundation of China under Grant No. 11505007.

Author Contributions All work was carried out under the supervision of G.L.L., F.G.D., Y.C.C., and F.N. B.X.W. and T.X. performed the experiments. Y.C.C., B.X.W., Q.A., M.J.T. and N.L. analysed the experimental data. M.J.T., Q.A. and N.L. wrote the HEOM simulation code and performed the numerical simulation. Y.C.C. and Q.A. contributed the theory. All authors contributed to writing the manuscript.

Author Information Reprints and permissions information is available at www.nature.com/reprints. The authors declare no competing financial interests. Readers are welcome to comment on the online version of the paper. Correspondence and requests for materials should be addressed to G.L.L.(gllong@tsinghua.edu.cn) and F.G.D.(fgdeng@bnu.edu.cn).

## Methods

Physical parameters. The Hamiltonian implemented in NMR simulation is the same as the photosynthetic Hamiltonian but the unit cm divided by . The diagonal terms are kHz, kHz, kHz, and kHz. The inter-level couplings are the nearest neighboring couplings kHz, kHz, the next-nearest-neighboring couplings kHz, and the coupling between the two ends kHz. The temperature of the photosynthetic energy transfer and NMR experiment are respectively K and K. The reorganization energy of the bath is cm. The cutoff frequency of the bath is cm. Both the Hamiltonian and bath parameters for NMR experiment are those for the photosynthetic energy transfer which are scaled down by a factor of .

Initialization. Starting from the thermal equilibrium state,
we prepare a pseudo-pure state ^{28, 29}
using the spatial-average technique ^{29}.

Measurement. The goal is to acquire probability
distributions of four states , ,
, and after the pulse sequences that we
designed are implemented on the two-qubit NMR system, namely, four
diagonal values of the final density matrix. The density matrices of the
output states are reconstructed completely via quantum state
tomography ^{30}. Therefore, the density
matrix of the system can be estimated from ensemble averages of a
set of observables. For the one-qubit system, the observable set is
. Here, , ,
, . For the two-qubit system, the observable
set is . In our
experiments, the complete density matrix tomography is not
necessary. All we need is to perform two experiments in which the
reading-out pulses and are
respectively implemented on the final states of H and
C and the corresponding qubits that need to be observed are
respectively H and C. The points in
Figs. 2 is obtained by averaging over random realizations.

Supplementary information for ”Quantum simulation of photosynthetic energy transfer”

## Appendix A Ramsey-like Dynamics in Photosynthesis

In this section, we provide a detailed derivation for the Ramsey-like dynamics in photosynthesis ^{31, 32}. We assume the total Hamiltonian to be

(10) |

where we have assumed that the dimer is subject to two local harmonic-oscillator baths with the same parameters.

In the experiment, the system is initialized to and followed by a pulse, i.e.,

(11) |

And the bath is in thermal equilibrium, i.e.,

(12) |

where the partition function of th bath mode is

(13) |

Then, the system evolves under the Hamiltonian (10) for a time interval and thus results in

(14) | |||||

Finally, after applying a reverse pulse, we measure the population of in the final state, i.e.,

(15) |

Thus, the populations of reads

(16) |

The off-diagonal element can be calculated as

(17) | |||||

where

(18) | |||||

(19) | |||||

(20) | |||||

(21) | |||||

Hereafter, we shall explicitly give the expression of as

(22) | |||||

where the displacement operator is

(23) |

By using the identity

(24) |

is simplified as

(25) | |||||

where in the last line we have used the Baker-Hausdorff formula ^{36}

(26) |

Then, we apply the identity

(27) |

to the above equation, we obtain as

(28) | |||||

By inserting into , we have

(29) | |||||

where the lineshape function reads

(30) |

The population of reads

(31) |

Since the spectral density is defined as

(32) | |||||

with being density of states of bath, the lineshape function can explicitly given as

(33) | |||||

where we assumed a Debye-form spectral density

(34) |

with and being the reorganization energy and cutoff frequency respectively.

By using a Matsubara expansion ^{48}, the lineshape function is explicitly calculated as

(35) |

where

(36) |

In Sec. C, we will demonstrate that in order to simulate the photosynthetic dynamics in NMR, the following relations should be fulfilled

(37) | |||

(38) |

## Appendix B Hierarchical Equations of Motion (HEOM)

The hierarchical equations of motion (HEOM) formalism has become an important method for studying quantum open systems ^{37, 38, 39}. In this section, we describe the application of the HEOM method for studying the excitation energy transfer (EET) in photosynthetic systems ^{38, 39}.

We discuss the EET dynamics in a photosynthetic complex containing four pigments, and each pigment is modeled by a two-level system. The following Frenkel exciton Hamiltonian ^{31, 40}, studying EET dynamics, consists of three parts,

(39) |

where

(40) | |||||

(41) | |||||

(42) | |||||

In the above, represents the state where only the th pigment is in its electronic excited state and all others are in their electronic ground state. Moreover, this

(43) |

is the so-called site energy of the th pigment, where is the excited electronic energy of the th pigment in the absence of phonons and is the reorganization energy of the th pigment. Furthermore, is the electronic coupling between pigments and . Also, , and are the frequency, dimensionless coordinate, and conjugate momentum of the th phonon mode, respectively. Here,

(44) | |||||

(45) |

with being the coupling constant between the th pigment and th phonon mode. For simplicity, we assume that the phonon modes associated with different pigments are uncorrelated.

The reduced density operator of the system

(46) |

with being the density operator for the total system can adequately describe the EET dynamics. At the initial time , we assume that the total system is in the factorized product state of the form

(47) |

with

(48) |

In accordance to the vertical Franck-Condon transition ^{38, 39},
the initial condition (47) is appropriate in electronic excitation processes. In this work,
we adopt the spectral density of the overdamped Brownian oscillator model,

(49) |

to describe the coupling between the th pigment and the environmental phonons. For this modeling, the timescale of the phonon relaxation is simply,

(50) |

According to the reorganization dynamics, one can determine the reorganization energy .

For high temperatures , the following hierarchically coupled equations of motion for the reduced density operator with the overdamped Brownian oscillator model is given by

(51) |

where are three sets of nonnegative integers, i.e.,

(52) | |||||

(53) |

The phonon-induced relaxation operators are written by

(54) | |||||

(55) |

where

(56) | |||||

(57) |

are the hyper-operator notations. In addition,

(58) |

and the other are auxiliary operators considering the fluctuation and dissipation. The Liouvillian operator corresponds to the electronic Hamiltonian .

We terminate Eq. (51), when the integers ’s satisfy

(59) |

where is a characteristic frequency of the system dynamics ^{38}. The required number of auxiliary density operators is given by

(60) |

## Appendix C Dynamics in Classical Pure-dephasing Noise

### c.1 General Case

In this section, inspired by Ref. ^{33}, we provide a detailed calculation for the dynamics in the classical pure-dephasing noise.
The total Hamiltonian

(61) |

is divided into two parts, i.e. the control Hamiltonian

(62) |

and the noise Hamiltonian

(63) |

where

(64) | |||||

(65) | |||||

(66) |

In the rotating frame with respect to

(67) |

the noise Hamiltonian reads

(68) |

And the propagator in this frame is correspondingly

(69) |

Therefore, transformed back to the Schrödinger picture, the propagator is written as

(70) |

Let us now consider

(71) | |||||

where

(72) |

and in the last line of Eq. (71) we have used the relation

(73) |

Hereafter, we shall use the compact definition

(74) |

where

(75) |

When

(76) | |||||

(77) |

according to the Magnus expansion ^{}