NMR is an excellent platform for quantum simulation since it is easy to operate and it can have long coherence times.22,23 In this paper, NMR is utilized to simulate the quantum coherent dynamics in photosynthetic light harvesting. As a prototype, a tetramer including four chlorophylls,13 as schematically illustrated in Fig. 1a, is employed in the NMR quantum simulation. In a previous investigation,13 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

$$H_{{\mathrm{EET}}} = \mathop {\sum}\limits_{i = 1}^4 {\kern 1pt} \varepsilon _i\left| i \right\rangle \left\langle i \right| + \mathop {\sum}\limits_{i

e j = 1}^4 {\kern 1pt} J_{ij}\left| i \right\rangle \left\langle j \right|,$$ (1)

where \(\left| i \right\rangle\) (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 \(\left| i \right\rangle\), J ij is the excitonic interaction between sites i and j. In our quantum simulations, we adopt the site energies13 ε 1 = 13,000 cm−1, ε 2 = 12,900 cm−1, ε 3 = 12,300 cm−1, and ε 4 = 12,200 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.

Fig. 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 13C-labeled chloroform molecule, where the H and C nuclear spins are chosen as the two qubits Full size image

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 \(\left| i \right\rangle\) is encoded as a two-qubit product state, i.e., \(\left| {00} \right\rangle\), \(\left| {01} \right\rangle\), \(\left| {10} \right\rangle\), and \(\left| {11} \right\rangle\). By a straightforward calculation, the Frenkel-exciton Hamiltonian can be mapped to the NMR Hamiltonian as

$$\begin{array}{*{20}{l}} {H_{{\mathrm{NMR}}}} \hfill & = \hfill & {\frac{{\varepsilon _1^\prime + \varepsilon _2^\prime - \varepsilon _3^\prime - \varepsilon _4^\prime }}{4}\sigma _1^z + \frac{{\varepsilon _1^\prime - \varepsilon _2^\prime + \varepsilon _3^\prime - \varepsilon _4^\prime }}{4}\sigma _2^z} \hfill \\ {} \hfill & {} \hfill & { + \frac{{J_{13}^\prime + J_{24}^\prime }}{2}\sigma _1^x + \frac{{J_{12}^\prime + J_{34}^\prime }}{2}\sigma _2^x + \frac{{J_{14}^\prime + J_{23}^\prime }}{2}\sigma _1^x\sigma _2^x} \hfill \\ {} \hfill & {} \hfill & { + \frac{{J_{23}^\prime - J_{14}^\prime }}{2}\sigma _1^y\sigma _2^y + \frac{{J_{13}^\prime - J_{24}^\prime }}{2}\sigma _1^x\sigma _2^z} \hfill \\ {} \hfill & {} \hfill & { + \frac{{J_{12}^\prime - J_{34}^\prime }}{2}\sigma _1^z\sigma _2^x,} \hfill \end{array}$$ (2)

where \(\sigma _j^u\) (j = 1, 2, u = x, y, z) is the Pauli operator for qubit j, numerically \(\varepsilon _j^\prime = \pi \varepsilon _j{\mathrm{/}}10\) and \(J_{ij}^\prime = \pi J_{ij}{\mathrm{/}}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 × 108/π.

The excitonic coupling J ij between sites i and j makes the exciton energy hop in both directions, i.e., \(\left| i \right\rangle\) ⇆ \(\left| j \right\rangle\). 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

$$H_{{\mathrm{SB}}} = \mathop {\sum}\limits_{i,k} {\kern 1pt} g_{ik}\left| i \right\rangle \left\langle i \right|\left( {a_{ik}^\dagger + a_{ik}} \right),$$ (3)

where \(a_{ik}^\dagger\) 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

$$G_{{\mathrm{EET}}}(\omega ) = \mathop {\sum}\limits_k {\kern 1pt} g_{ik}^2\delta \left( {\omega - \omega _k} \right),$$ (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

$$H_{{\mathrm{PDN}}} = \vec B(t) \cdot \vec \sigma$$ (5)

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

$$\vec B(t) = {\mathrm{\Omega }}_0{\kern 1pt} {\mathrm{cos}}\left[ {\omega _\mu t + \phi _N(t)} \right]\hat z,$$ (6)

$$\phi _N(t) = \alpha \mathop {\sum}\limits_{j = 1}^J {\kern 1pt} F(j){\mathrm{sin}}\left( {\omega _jt + \psi _j} \right),$$ (7)

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 power-spectral density \(S_z(\omega ) \equiv {\int} {\kern 1pt} d\tau \left\langle {\dot \phi _N(t + \tau )\dot \phi _N(t)} \right\rangle e^{i\omega \tau }\) is the Fourier transform of the autocorrelation function of \(\dot \phi _N(t)\),

$$S_z(\omega ) = \frac{{\pi \alpha ^2\omega _0^2}}{2}\mathop {\sum}\limits_{j = 1}^J {\kern 1pt} j^2F^2(j)\left[ {\delta \left( {\omega - \omega _j} \right) + \delta \left( {\omega + \omega _j} \right)} \right].$$ (8)

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 Drude–Lorentz spectral densities. For the details, please refer to the Supplementary Information.

As illustrated in Fig. 1b, 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

$$H_{{\mathrm{CHCl}}_{\mathrm{3}}} = \pi \omega _1\sigma _1^z + \pi \omega _2\sigma _2^z + \frac{{\pi J}}{2}\sigma _1^z\sigma _2^z,$$ (9)

where ω 1 = 3206.5 Hz, ω 2 = 7787.9 Hz are the chemical shifts of the two spins, and J = 215.1 Hz (ref. 24). 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,25 which has been successfully applied to a number of quantum simulations in NMR,23,26,27 to mimic the quantum dynamics of light harvesting.

Before studying the population dynamics, we shall explore the characteristics of the artificially injected noise. Generally, the noise is characterized by Ramsey interference measurements, as shown in Fig. 2. When there is no noise applied to the system, we would expect coherent evolution without damping. However, as the artificial noise is applied, the oscillations become damped. The envelope of the damped oscillations is fully determined by the noise strength. Interestingly, for the parameters employed in this work for a broad-frequency background bath, the NMR experimental results exactly fit the theoretical prediction by HEOM. This suggests that our technique for injecting artificial noise can, in this limit, faithfully reproduce the effect of the protein environment on natural photosynthetic complexes.

Fig. 2 The decay of the population of the state \(\left| 0 \right\rangle\). The theoretical simulation adopting HEOM is shown by the blue curve, and the experimental data by the red dots. In our experiment, the parameters are ε D − ε A = 2π × 15 kHz, dt = 20 μs, and M = 50. The parameters of the Drude–Lorentz noise are γ NMR = 2π × 45 kHz, λ NMR = 2π × 0.01 kHz, T EET = 300 K, ω 0 = 2π × 5 Hz, and ω J = 2π × 25 kHz Full size image

In Fig. 3a, 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. 3b, 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) and (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. In Fig. 4, we show the effect on the simulation error of the number of realizations in each ensemble. Both results, by NMR and by HEOM, converge as M increases, such that the difference between them can hardly be noticed when M ~ 500. 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) and (7)).

Fig. 3 Simulation of the energy transfer governed by H NMR + H PDN for M = 150 random realizations. a Long-time quantum dynamics; b short-time quantum dynamics. The symbols show the experimental data, and the curves are obtained from the numerical simulation using the HEOM. In all figures, the horizontal coordinates of the curves for EET dynamics have been magnified by 3 × 108/π times Full size image

Fig. 4 Effect of the ensemble size on the simulations. The symbols show the theoretical results and the curves show the HEOM results. The number of realizations in each ensemble are M = 50 a, 500 b, and 5000 c Full size image

It has been suggested in the literature6,28,29 that the long-lasting electronic coherence observed in experiments might be attributed to strong coupling to vibrational modes. Therefore, it is interesting to consider whether it is possible to simulate some features of underdamped vibrational modes in the NMR platform. In photosynthesis, the coupling to a vibrational mode can be described by \(\mathop {\sum}

olimits_i {\kern 1pt} g_{iv}\left| i \right\rangle \left\langle i \right|\left( {a_{iv}^\dagger + a_{iv}} \right)\). In order to simulate the underdamped vibrational mode in our protocol, we can supplement an additional term, i.e., αF(v) sin(ω v t) into Eq. (7). In Fig. 5, we show the numerical simulation for the above 4-level system including such a “semi-classical” vibrational mode. Because we set the vibrational mode to be closely-resonant with the two lower exciton states, there is significantly prolonged coherence in the two lower energy states in the short-time regime, as compared to Fig. 3. Since there is no random phase associated with this modulation, and quantum properties of the vibrational mode are neglected, each realization will evolve uniformly in Fig. 5. Therefore, the long-lasting coherence is driven by the nearly-resonant driving. Going beyond a semi-classical approach, and accurately capturing the quantum features of this mode, requires the use of ancillary degrees of freedom, as discussed in the conclusions.