Quantum simulation of photosynthetic energy transfer

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 N- chromophore photosynthetic complex, with arbitrary structure and bath spectral density, can be effectively simulated by a system with log2 N qubits. The computational cost of simulating such a system with a theoretical tool, like the hierarchical equation of motion, which is exponential in N, can be potentially reduced to requiring a just polynomial number of qubits N 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.

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 Nchromophore photosynthetic complex, with arbitrary structure and bath spectral density, can be effectively simulated by a system with log 2 N qubits. The computational cost of simulating such a system with a theoretical tool, like the hierarchical equation of motion, which is exponential in N , can be potentially reduced to requiring a just polynomial number of qubits N 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 1,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 2 . Much effort has been made to reveal the effects of quantum coherence on efficient energy transfer 2,5-9, 19 . In order to try to mimic EET, a Frenkel-exciton Hamiltonian is required and this can be studied with quantum chemistry approaches 1 , 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 1,12,19 , the hierarchical equation of motion (HEOM) yields a numerically exact solution at the cost of considerable computing time 7,19 . 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 7,19 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 N qubits can powerfully mimic the quantum dynamics of 2 N 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 newlydeveloped 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 4 and NMR 5 . 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 Donor Acceptor (a) (b) Figure 1: Photosynthetic chromophore arrangement and physical system for NMR simulation. (a) Linear geometry with four chromophores for photosynthetic energy transfer; (b) Chemical structure for a 13 C-labeled chloroform molecule, where the H and C nuclear spins are chosen as the two qubits.
transfer. The EET in photosynthesis is described by the Frenkel-exciton Hamiltonian where = 1 4) is the state with a single excitation at site and other sites at the ground state, is the site energy of ij is the excitonic interaction between sites and . In our quantum simulations, we adopt the site energies 12 = 13000 cm = 12900 cm = 12300 cm , and = 12200 cm , and the couplings 12 34 = 126 cm 13 24 = 16 cm 23 = 132 cm , and 14 = 5 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 log 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. 00 01 10 , and 11 . By a straightforward calculation, the Frenkel-exciton Hamiltonian can be mapped to the NMR Hamiltonian as (2) where = 1 2, x, y, z) is the Pauli operator for qubit , numerically πε 10 and ij πJ ij 10, 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 1 cm π/10 kHz) = 3 10 /π The excitonic coupling ij between sites and makes the exciton energy hop in both directions, i.e. Apart from this, the system-bath coup the energy flow towards an energy tra tured photon energy is converted into The interaction between the system a synthetic complexes can be described where ik is the creation operator fo mode of chlorophyll with coupling will induce pure-dephasing on the -th it is in the excited state. Generally couplings are given by the spectral de EET ) = ik which we assume identical for all chrom ical photosynthetic complexes, the syst are of the same order as the intra-sys order to mimic the effects of noise, w engineering technique which has been mented in ion traps and NMR. The im pure-dephasing Hamiltonian where Ω is the constant amplitude o with driving frequency is a ran jω with and Jω being base cies respectively, and is a global sc power-spectral density dτ is the Fourier transform of the autoc of ), By appropriately choosing ) and t bitrary power-spectral density of the ized, e.g. white, Ohmic, or Debye spe the details, please refer to the Supple tion.
As illustrated in Fig. 1(b), the nu carbon atom and hydrogen atom in a ch are chosen to encode the two qubits wi written as CHCl πω πω Figure 1: Photosynthetic chromophore arrangement and physical system for NMR simulation. (a) Linear geometry with four chromophores for photosynthetic energy transfer; (b) Chemical structure for a 13 C-labeled chloroform molecule, where the H and C nuclear spins are chosen as the two qubits.
in photosynthesis is described by the Frenkel-exciton Hamiltonian where |i (i = 1, 2, 3, 4) is the state with a single excitation at site i and other sites at the ground state, ε i is the site energy of |i , J ij is the excitonic interaction between sites i and j. In our quantum simulations, we adopt the site energies 12 ε 1 = 13000 cm −1 , ε 2 = 12900 cm −1 , ε 3 = 12300 cm −1 , and ε 4 = 12200 cm −1 , and the couplings J 12 = J 34 = 126 cm −1 , J 13 = J 24 = 16 cm −1 , J 23 = 132 cm −1 , and J 14 = 5 cm −1 , which are typical parameters in photosynthetic systems. For a photosynthetic complex with N chlorophylls, there are only N single-excitation states involved in the energy transfer. Therefore, only log 2 N 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 |i is encoded as a two-qubit product state, i.e. |00 , |01 , |10 , and |11 . By a straightforward calculation, the Frenkel-exciton Hamiltonian can be mapped to the NMR Hamiltonian as where σ u j (j = 1, 2, u = x, y, z) is the Pauli operator for qubit j, numerically ε ′ j = πε j /10 and J ′ ij = πJ ij /10, but the dimension cm −1 should be replaced by kHz. In other words, all realistic parameters have been scaled down in energy by a factor of 1 cm −1 /(π/10 kHz) = 3 × 10 8 /π.
The excitonic coupling J ij between sites i and j makes the exciton energy hop in both directions, i.e. |i ⇆ |j . 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 where a † ik is the creation operator for the kth phonon mode of chlorophyll i with coupling strength g ik . H SB will induce pure-dephasing on the i-th chromophore when it is in the excited state. Generally, the system-bath couplings are given by the spectral density 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 relies on generating stochastic errors by performing phase modulations on a constant-amplitude carrier, i.e.
where Ω 0 is the constant amplitude of a magnetic field with driving frequency ω µ , ψ j is a random number, ω j = jω 0 with ω 0 and ω J = Jω 0 being base and cutoff frequencies respectively, and α is a global scaling factor. The powerspectral density S z (ω) ≡ dτ φ N (t + τ )φ N (t) e iωτ is the Fourier transform of the autocorrelation function ofφ N (t), By appropriately choosing F (j) and the cutoff J, 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 where ω 1 = 3206.5 Hz, ω 2 = 7787.9 Hz are the chemical shifts of the two spins, and J = 215.1 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 14 , which has been successfully applied to a number of quantum simulations in NMR 19,22,23 , to mimic the quantum dynamics of light harvesting. In Fig. 2(a), the quantum coherent energy transfer for the above Hamiltonian H NMR + H PDN , 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 M 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 M 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 7,19 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 7,19 . Our NMR simulations do not have such shortcomings, as it scales polynomially with respect to the system size 21,22 . 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 ∼ 100 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 ∼ 300 sites) can be faithfully simulated by a 9-qubit NMR quantum computer, which is clearly attainable by the NMR technique 22 . 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. 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 H NMR is the same as the photosynthetic Hamiltonian H EET but the unit cm −1 divided by 3 × 10 4 /π. The diagonal terms are H 11 NMR = 2π × 650 kHz, H 22 NMR = 2π × 645 kHz, H 33 NMR = 2π × 615 kHz, and H 44 NMR = 2π × 610 kHz. The inter-level couplings are the nearest neighboring couplings H 12 NMR = H 34 NMR = 2π × 6.3040 kHz, H 23 NMR = 2π × 6.5950 kHz, the next-nearest-neighboring couplings H 13 NMR = H 24 NMR = 2π × 0.8059 kHz, and the coupling between the two ends H 14 NMR = 2π × 0.2370 kHz. The temperature of the photosynthetic energy transfer and NMR experiment are respectively T EET = 3 × 10 4 K and T NMR = 5 × 10 −5 K. The reorganization energy of the bath is λ EET = 0.2 cm −1 . The cutoff frequency of the bath is γ EET = 900 cm −1 . Both the Hamiltonian and bath parameters for NMR experiment are those for the photosynthetic energy transfer which are scaled down by a factor of 3 × 10 8 /π. Initialization. Starting from the thermal equilibrium state, we prepare a pseudo-pure state 28,29 using the spatialaverage technique 29 . Measurement. The goal is to acquire probability distributions of four states |00 , |01 , |10 , and |11 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 {σ i }(i = 0, 1, 2, 3). Here, σ 0 = I, σ 1 = σ x , σ 2 = σ y , σ 3 = σ z . For the two-qubit system, the observable set is {σ i ⊗ σ j }(i, j = 0, 1, 2, 3). 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 exp(−iπσ y /4) ⊗ I and I ⊗ exp(−iπσ y /4) are respectively implemented on the final states of 1 H and 13 C and the corresponding qubits that need to be observed are respectively 1 H and 13 C. The points in Figs. 2 is obtained by averaging over M = 150 random realizations.
Supplementary information for "Quantum simulation of photosynthetic energy transfer"

RAMSEY-LIKE DYNAMICS IN PHOTOSYNTHESIS
In this section, we provide a detailed derivation for the Ramsey-like dynamics in photosynthesis 1,2 . We assume the total Hamiltonian to be 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 |1 and followed by a π/2 pulse, i.e., And the bath is in thermal equilibrium, i.e., where the partition function of kth bath mode is Then, the system evolves under the Hamiltonian (10) for a time interval t and thus results in Finally, after applying a reverse π/2 pulse, we measure the population of |1 in the final state, i.e., Thus, the populations of |1 reads The off-diagonal element can be calculated as where 8 Hereafter, we shall explicitly give the expression of I where the displacement operator is By using the identity is simplified as where in the last line we have used the Baker-Hausdorff formula 6 Then, we apply the identity to the above equation, we obtain I where the lineshape function reads The population of |1 reads Since the spectral density is defined as with ρ(ω k ) being density of states of bath, the lineshape function can explicitly given as where we assumed a Debye-form spectral density with λ and Λ being the reorganization energy and cutoff frequency respectively. By using a Matsubara expansion 18 , the lineshape function is explicitly calculated as where ν n = 2πn β .
In Sec. , we will demonstrate that in order to simulate the photosynthetic dynamics in NMR, the following relations should be fulfilled

HIERARCHICAL EQUATIONS OF MOTION (HEOM)
The hierarchical equations of motion (HEOM) formalism has become an important method for studying quantum open systems [7][8][9] . In this section, we describe the application of the HEOM method for studying the excitation energy transfer (EET) in photosynthetic systems 8,9 .
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 1,10 , studying EET dynamics, consists of three parts, where In the above, |j represents the state where only the jth pigment is in its electronic excited state and all others are in their electronic ground state. Moreover, this is the so-called site energy of the jth pigment, where ε 0 j is the excited electronic energy of the jth pigment in the absence of phonons and λ j is the reorganization energy of the jth pigment. Furthermore, J jk is the electronic coupling between pigments i and j. Also, ω m , p jm and q jm are the frequency, dimensionless coordinate, and conjugate momentum of the mth phonon mode, respectively. Here, with c jm being the coupling constant between the jth pigment and mth phonon mode. For simplicity, we assume that the phonon modes associated with different pigments are uncorrelated. The reduced density operator of the system with ρ tot being the density operator for the total system can adequately describe the EET dynamics. At the initial time t = 0, we assume that the total system is in the factorized product state of the form In accordance to the vertical Franck-Condon transition 8,9 , the initial condition (47) is appropriate in electronic excitation processes. In this work, we adopt the spectral density of the overdamped Brownian oscillator model, to describe the coupling between the jth pigment and the environmental phonons. For this modeling, the timescale of the phonon relaxation is simply, According to the reorganization dynamics, one can determine the reorganization energy λ j . For high temperatures βγ j < 1, the following hierarchically coupled equations of motion for the reduced density operator with the overdamped Brownian oscillator model is given by where n, n j± are three sets of nonnegative integers, i.e., n = (n 1 , n 2 , n 3 , n 4 ) , n j± = (n 1 , · · · , n j ± 1, · · · , n 4 ) .
The phonon-induced relaxation operators are written by where are the hyper-operator notations. In addition, and the other σ(n = 0, t) are auxiliary operators considering the fluctuation and dissipation. The Liouvillian operator ℓ e corresponds to the electronic Hamiltonian H e . We terminate Eq. (51), when the integers n j 's satisfy where ω e is a characteristic frequency of the system dynamics ℓ e 8 . The required number of auxiliary density operators σ(n, t) is given by

General Case
In this section, inspired by Ref. 3 , we provide a detailed calculation for the dynamics in the classical pure-dephasing noise. The total Hamiltonian is divided into two parts, i.e. the control Hamiltonian and the noise Hamiltonian where In the rotating frame with respect to the noise Hamiltonian reads And the propagator in this frame is correspondingly Therefore, transformed back to the Schrödinger picture, the propagator is written as Let us now consider where and in the last line of Eq. (71) we have used the relation Hereafter, we shall use the compact definition where When according to the Magnus expansion 15 , we have By using the identity the propagator (76) can be rewritten as where with To calculate the fidelity of the operation described by we use the Hilbert-Schmidt inner product to measure the fidelity as where a is the modulus of the vector a(τ ), and · · · is averaged over all possible noise trajectories. In the above equation, the lowest order term is Thus, the fidelity can be expanded as where we have used the relation By introducing the Fourier transform of the cross-power spectrum we have where we have defined

Ramsey Fringes
In the following, we shall consider a special case where [H 0 (t), H c (t)] = 0. At the end of this subsection, we will provide the deviation for Ramsey-interferometer experiment. In this case, we assume the total Hamiltonian as Thus, the control Hamiltonian is and the noise Hamiltonian is In the rotating frame with respect to the noise Hamiltonian reads because And the propagator in this frame is correspondingly Therefore, transformed back to the Schrödinger picture, the propagator is written as In the experiment, the system is initialized to |0 and followed by a π/2 pulse, i.e.
Then, the system evolves under the Hamiltonian (97) for a time interval t and thus results in where Finally, after applying a reverse π/2 pulse, we measure the population of |0 in the final state, i.e. 16 Thus, the population of |0 reads P 0 (t) = cos 2 φ(t) where · · · is averaged over all possible random realizations. If we further assume a Gaussian noise, for any positive integer n, P 0 (t) can be simplified as For the lowest nontrivial order n = 1, we have where we have introduced the Fourier transform of β z (τ 1 )β z (τ 2 ) as with β z (τ 1 )β z (τ 2 ) being only dependent on the time interval τ 2 − τ 1 . For the order with n = 2, we have For the order with arbitrary integer n, we have where there are (2n)!/(2 n n!) terms in the second line according to Isserlis' theorem if it is a Gaussian noise 16 . To conclude, This predicts that before decaying to the steady value 1/2 in the long run, P 0 (t) will experience oscillations with frequency ω L . We assume that where the ψ j 's are random numbers. According to Ref. 16 , the ensemble average is equivalent to the time average for a wide-sense-stationary random process. In this case, the two-time correlation function reads which does not depend on t but τ . The power spectral density is the Fourier transform of the correlation function, i.e. where The power spectral density is a set of equally-spaced peaks with distance ω 0 and height α we have and thus S zz (ω) is an Ohmic spectral density of step-function form with cutoff frequency If we set we have and thus S zz (ω) is the Debye-Drude spectral density of the step-function form with cutoff frequency ω J . The transverse relaxation time T 2 is defined by For photosynthesis, the decoherence is determined by the real part of the lineshape function with spectral density where Therefore, in order to simulate photosynthetic dynamics in NMR, we should relate the following two quantities and cutoff frequencies in two spectra are equal, i.e.

TECHNIQUE OF ARTIFICIALLY INJECTING NOISE
Here we introduce a method of artificially injecting noises in NMR and ion trap systems, including dephasing noise and amplitude noise 4,5 .

Dephasing Noise
Dephasing noise comes from the inhomogeneous and non-static magnetic field in NMR systems. The corresponding Hamiltonian can be written as β z (t)σ z where α i (i = x, y, z) is the noise amplitude and φ j is a random phase. N ω 0 determines the high frequency cutoff and ω 0 is the base frequency with ω j = jω 0 . The types of noise rely on the function F (ω j ). The two-time correlation function for β z (t) is then written as which does not depend on t but on τ . Applying the Wiener-Khintchine theorem 17 , we then obtain the power spectral density which can describe the energy distribution of the stochastic signal in the frequency domain by Fourier transform Hence, we can use the model of power spectral density to reverse the noise distribution in the time domain. For instance, if we want to simulate the power spectral density for S(ω) ∼ ω p , then the modulation function F (ω j ) = (ω j ) p/2−1 . Taking Eq. (132) into above, The initial state |ψ(0) = α|0 + β|1 with the dephasing noise of the Hamiltonian β z (t)σ z in the rotating frame after time τ will become where △θ τ is the integral of β z (t) . Hence we just rotate the angle △θ τ along the z-axis at the desired point to realize the evolution of the quantum system in the dephasing environment.

Amplitude Noise
Similarly, we can obtain the β x (t) of the amplitude noise as a result of the amplitude fluctuation of the control field, and its power spectral density is

Parameters for Dephasing Noise of Debye-Drude Form
For the Debye spectrum of the dephasing noise, taking Re[g(t)] equal to χ(t), we obtain In short, as long as we know the power spectral density of the noise, we can then obtain the time-varying β(t) and χ(t). Table I shows F (ω j ) for distinct types of dephasing and amplitude noises.

Parameters for Dephasing Noise of Arbitrary Form
For the general spectrum J(ω) of the dephasing noise, we make Re[g(t)] and χ(t) be equal according to Eqs. (33), (132), (137) After simplifying the above formula, we can obtain

EXPERIMENTAL DETAILS AND RESULTS
Experiments are carried out at room temperature using a Bruker Avance III 400 MHz spectrometer. The sample is chloroform dissolved in d6-acetone as a two-qubit NMR quantum processor where H is the first qubit and C is second qubit. The internal Hamiltonian of the two-qubit system can be described as where is the J-coupling strength between two spins. The experimental process is divided into three steps, as shown in Fig. 3.

Preparation of the Pseudo-pure State
The thermal equilibrium state for the two qubit system is where I is a 4 × 4 identity matrix, describes the polarization, and γ H and γ C represent the gyromagnetic ratios of the 1 H and 13 C nuclei, respectively. The spatial average technique 12 is used to obtain a pseudo-pure state and the related pulse sequence is depicted in Fig. 4. Thus we only focus on the part |00 as the entire system behaves since the identity part does not influence the unitary operations or measurements in NMR experiments.

Evolution of the Hamiltonian with Debye noise
The total Hamiltonian for simulating photosynthetic EET in the NMR system is where H S is the system Hamiltonian, and β i (t) (i = 1, 2) are the time-dependent Debye noises. In experiments, the evolution can have L discretized steps, and the evolution time is t = L∆t with where H i is the time-independent Hamiltonian at point t i = i∆t. Note U (t) is calculated by the gradient ascent pulse engineering (GRAPE) method with 5 ms of each pulse to reduce the accumulated pulse errors in experiments.

Measure the Probability Distribution of Four States
Our goal now is to acquire probability distributions of four states |00 , |01 , |10 , |11 ; namely, the four diagonal values of the final density matrix. The density matrices of the output states are reconstructed completely via quantum state tomography (QST) 13 . In the QST theory, 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 {σ i } (i = 0, 1, 2, 3), where σ 0 = I is the identity, σ 1 = X, σ 2 = Y , σ 3 = Z are the Pauli matrices. The NMR signal is which is oscillating at the frequency ω and X and Y are obtained in practice by Fourier transforming S(t) and integrating the real and imaginary spectra, respectively. The signal becomes after applying exp[−iπY /4]. The density operator of one-qubit can be estimated by For the two-qubit system, the observable set is {σ i ⊗ σ j }(i, j = 0, 1, 2, 3). 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 exp(−iπY /4) ⊗ I and I ⊗ exp(−iπY /4) are respectively implemented on the final states of 1 H and 13 C and the corresponding qubits that need to be observed are respectively 1 H and 13 C. Then the probability distribution of four states are obtained by the results of measuring ZI , IZ , ZZ . Then the probability distribution of four states are respectively

GRADIENT ASCENT PULSE ENGINEERING (GRAPE) ALGORITHM
Here we describe the GRAPE technique proposed by Glaser et al. 14 which has been frequently used in NMR experiments. For an n-qubit NMR system, the total Hamiltonian contains the internal term and the radio frequency (RF) term where B k and φ k are the amplitude and phase of the control field on the kth nuclear spin. The goal of the GRAPE technique is to find the optimal parameters B k and φ k of the RF field by iteration to control the designed evolution U T very close to desired target evolution U D . Assuming that the total time of RF field is T , which is divided into L discrete segments. The time of each segment is ∆t = T /L and the time propagator of the jth segment can be expressed as Thus, the total evolution is 24 The fidelity to the target evolution U T can be expressed as which is also called the fitness function. The GRAPE algorithm considers the fidelity F as the extreme value optimization of the multi-function. We calculate the gradient function to first order, The fitness functions can be increased in the gradient iteration, where ε is a suitable and small step size. The GRAPE procedure starts from an initial guess input and evaluates the corresponding gradients g k x,y (j) and then keeps iterating until the fitness function reaches the desired value.  In Sec. and , we provide the deviation for the Ramsey experiments in photosynthetic light harvesting and NMR respectively. In this section, we will compare the data by these two different approaches, i.e. the theoretical Ramsey fringes by HEOM and the experimental Ramsey fringes by NMR.
Ramsey experiments are often applied to observe the coherence of the qubits under the impact of noise. In our experiments, we perform the Ramsey spectroscopy for a single qubit subject to the Debye-Drude noise. In particular, the experiment is divided into three steps: • (1) Preparing the |0 state and rotating to the xy plane by a π/2 pulse along the x-axis.
• (3) Applying a reverse π/2 pulse and measuring the population of the |0 state.
We observe that the decay time constant of the population of the |0 state is significantly reduced in the presence of the Debye noise, as shown in Fig. 5. Ramsey fringes are a fit to a cosine with a simple exponential decay envelope, see Eqs. (31) and (118). The result by the HEOM theory (blue curve in Fig. 5) is consistent with that in experiment (red dots in Fig. 5). In other words, we demonstrate the equality of χ(t) and Re[g(t)].

RESULT OF NUMERICAL SIMULATION
Before the experimental demonstration, we shall numerically demonstrate the photosynthetic light harvesting can be exactly mimicked by the NMR quantum simulation using the GRAPE algorithm.
In our numerical simulation, we used the following parameter, i.e., γ NMR = 2π × 45 kHz, λ NMR = 2π × 0.01 kHz, T EET = 3 × 10 4 K, T NMR = 5 × 10 −5 K. And we assume the other physical constants are one, that is, = 1.055 × 10 −34 J · s, k B = 1.381 × 10 −23 J/K. The Hamiltonian for four-pigments in the single-excitation subspace is  Figure S4. The dots show the GRAPE numerical results and the lines show the H time of (b) interval is ms.

COMPUTATIONAL COSTS OF
When the open quantum system consists of levels and baths, exponentials, there are [S19 S24

COMPUTATIONAL COSTS OF NMR AND HEOM
When the open quantum system consists of N levels or sites, each coupled to an independent bath (so N baths), and the correlation function of each bath contains K exponentials, there are 19,24 On the other hand, the computational cost of GRAPE is 21,22 where N is the number of energy levels involved in the energy transfer. Because the N -level photosynthetic light harvesting is simulated by log 2 N -qubit NMR, the computational cost has the potential to be effectively reduced from exponential in N by the HEOM to polynomial in N by GRAPE.

EFFECT OF NUMBER OF RANDOM REALIZATIONS AND ERROR ANALYSIS
Errors are small and mainly caused by imperfections in the initial-ground-state preparation and GRAPE pulses, which can be estimated by numerical simulations. The remaining errors may originate from, e.g., imperfections in the experimental quantum control, the static magnetic field, and the spectral integrals.
As shown in Sec. , in the derivation of quantum dynamics under the influence of noise, we have assumed that the average over random realizations is equivalent to the average over time.
The assumption is valid only if the number of random realizations M is in the infinite-M limit. In order to verify this assumption, we experimentally investigate the effect of number of random realizations on the NMR simulation, as shown in Figs. 7,8. As M increases from M = 50 to M = 150, the experimental simulation approaches closer and closer to the numerical simulation by the HEOM. In Fig. 9, we further compare the numerical results by the GRAPE and the HEOM. As M increases from M = 50 to M = 10 4 , the difference between the results by the GRAPE and the HEOM reduces. When M ≥ 500, the difference is hardly noticeable. Therefore, it is justified to mimic the photosynthetic energy transfer by the NMR with random realizations.
results by the GRAPE and the HEOM. As increases from = 50 to = 10 , the difference between the results by the GRAPE and the HEOM reduces. When 500, the difference is hardly noticeable. Therefore, it is justified to mimic the photosynthetic energy transfer by the NMR with random realizations.  difference is hardly noticeable. Therefore, it is justified to mimic the photos realizations.