# Simulation of Chemical Reaction Dynamics on an NMR Quantum Computer

Quantum simulation can beat current classical computers with minimally a few tens of qubits and will likely become the first practical use of a quantum computer. One promising application of quantum simulation is to attack challenging quantum chemistry problems. Here we report an experimental demonstration that a small nuclear-magnetic-resonance (NMR) quantum computer is already able to simulate the dynamics of a prototype chemical reaction. The experimental results agree well with classical simulations. We conclude that the quantum simulation of chemical reaction dynamics not computable on current classical computers is feasible in the near future.

##### Introduction.

In addition to offering general-purpose quantum algorithms with substantial speed-ups over classical algorithms [1] [e.g., Shor’s quantum factorizing algorithm [2]], a quantum computer can be used to simulate specific quantum systems with high efficiency [3]. This quantum simulation idea was first conceived by Feynman [4]. Lloyd proved that with quantum computation architecture, the required resource for quantum simulation scales polynomially with the size of the simulated system [5], as compared with the exponential scaling on classical computers. During the past years several quantum simulation algorithms designed for individual problems were proposed [6, 7, 8, 9, 10] and a part of them have been realized using physical systems such as NMR [11, 12, 13] or trapped-ions [14]. For quantum chemistry problems, Aspuru-Guzik et al. and Kassal et al. proposed quantum simulation algorithms to calculate stationary molecular properties [15] as well as chemical reaction rates [16], with the quantum simulation of the former experimentally implemented on both NMR [17] and photonic quantum computers [18]. In this work we aim at the quantum simulation of the more challenging side of quantum chemistry problems – chemical reaction dynamics, presenting an experimental NMR implementation for the first time.

Theoretical calculations of chemical reaction dynamics play an important role in understanding reaction mechanisms and in guiding the control of chemical reactions [19, 20]. On classical computers the computational cost for propagating the Schrödinger equation increases exponentially with the system size. Indeed, standard methods in studies of chemical reaction dynamics so far have dealt with up to 9 degrees of freedom (DOF) [21]. Some highly sophisticated approaches, such as the multi-configurational time-dependent Hartree (MCTDH) method [22], can treat dozens of DOF but various approximations are necessary. So generally speaking, classical computers are unable to perform dynamical simulations for large molecules. For example, for a 10–DOF system and if only 8 grid points are needed for the coordinate representation of each DOF, classical computation will have to store and operate 8^{10} data points, already a formidable task for current classical computers. By contrast, such a system size is manageable by a quantum computer with only 30 qubits. Furthermore, the whole set of data can be processed in parallel in quantum simulation.

In this report we demonstrate that the quantum dynamics of a laser-driven hydrogen transfer model reaction can be captured by a small NMR quantum simulator. Given the limited number of qubits, the potential energy curve is modeled by 8 grid points. The continuous reactant-to-product transformation observed in our quantum simulator is in remarkable agreement with a classical computation based also upon an 8-dimensional Hilbert space. To our knowledge, this is the first explicit implementation of the quantum simulation of a chemical reaction process. Theoretical methods and general experimental techniques described in this work should motivate next-generation simulations of chemical reaction dynamics with a larger number of qubits.

##### Theory.

Previously we were able to simulate the ground-state energy of Hydrogen molecule [17]. Here, to simulate chemical reaction dynamics, we consider a one-dimensional model of a laser-driven chemical reaction [23], namely, the isomerization reaction of nonsymmetric substituted malonaldehydes, depicted in (Fig. 1A). The system Hamiltonian in the presence of an external laser field is given by

H(t)=T+V+E(t){\quad\rm with\quad}E(t)=-\mu\varepsilon(t). | (1) |

In (Eq. 1), E(t) is the laser-molecule interaction Hamiltonian, \mu=eq is the dipole moment operator, \varepsilon(t) represents the driving electric field, {T}={p}^{2}/2m is the kinetic energy operator, and

{V}=\frac{\Delta}{2q_{0}}(q-q_{0})+\frac{V^{\ddagger}-\Delta/2}{q_{0}^{4}}(q-q% _{0})^{2}(q+q_{0})^{2} | (2) |

is a double-well potential of the system along the reaction coordinate. In (Eq. 2) V^{\ddagger} is the barrier height, \Delta gives the asymmetry of the two wells, and \pm q_{0} give the locations of the potential well minima. See the figure caption of (Fig. 1B) for more details of this model.

We first employ the split-operator method [24, 25] to obtain the propagator {U}(t+\delta t,t) associated with the time interval from t to t+\delta t. We then have

\displaystyle{U}(t+\delta t,t)\approx | \displaystyle\,e^{-\frac{i}{\hbar}{V}\delta t/2}e^{-\frac{i}{\hbar}{E}(t+% \delta t/2)\delta t/2}e^{-\frac{i}{\hbar}{T}\delta t} | |||

\displaystyle\times e^{-\frac{i}{\hbar}{E}(t+\delta t/2)\delta t/2}e^{-\frac{i% }{\hbar}{V}\delta t/2}. | (3) |

The unitary operator e^{-i{T}\delta t/\hbar} in (Eq. 3) is diagonal in the momentum representation whereas all the other operators are unitary and diagonal in the coordinate representation. Such {U}(t+\delta t,t) can be simulated in a rather simple fashion if we work with both representations and make transformations between them by quantum Fourier transform (QFT) operations. To take snapshots of the dynamics we divide the reaction process into 25 small time steps, with \delta t=1.5\,fs and the total duration t_{f}=37.5\,fs. The electric field of an ultrashort strong laser pulse is chosen as

\varepsilon(t)=\left\{\begin{array}[]{cc}\varepsilon_{0}\sin^{2}(\frac{\pi t}{% 2s_{1}});&\qquad 0\leq t\leq s_{1}\\ \varepsilon_{0};&\qquad s_{1}<t<s_{2}\\ \varepsilon_{0}\sin^{2}[\frac{\pi(t_{f}-t)}{2(t_{f}-s_{2})}];&\qquad s_{2}\leq t% \leq t_{f}\end{array}\right. | (4) |

with s_{1}=5\,fs and s_{2}=32.5\,fs. More details, including an error analysis of the split-operator technique, are given in the supplementary material. The reactant state at t=0 is assumed to be the ground-state \left|\phi_{0}\right\rangle of the bare Hamiltonian T+V, which is mainly localized in the left potential well. The wavefunction of the reacting system at later times is denoted by \left|\psi({t})\right\rangle. The product state of the reaction is taken as the first excited state \left|\phi_{1}\right\rangle of T+V, which is mainly localized in the right potential well.

With the system Hamiltonian, the initial reactant state, the product state, and the propagation method outlined above, the next step is
to encode the time-evolving wavefunction \left|\psi({t})\right\rangle and the T, V, E(t) operators by *n* qubits. To that end
we first obtain the expressions of these operators in representation of a set of N=2^{n} discretized position basis states.
The evolving state can then be encoded as

\displaystyle\left|\psi(t)\right\rangle= | \displaystyle\sum_{q=0}^{2^{n}-1}m_{q}(t)\left|q\right\rangle | |||

\displaystyle= | \displaystyle m_{0}(t)\left|0\cdots 00\right\rangle+...+m_{2^{n}-1}(t)\left|1% \cdots 11\right\rangle, | (5) |

and as a result the system operators become

\displaystyle{T} | \displaystyle\!\!= | \displaystyle\!\!\sum_{k_{1},\cdots,k_{n}=z,i}\alpha_{k_{1}\cdots k_{n}}\sigma% _{k_{1}}^{1}\sigma_{k_{2}}^{2}\cdots\sigma_{k_{n}}^{n}, | (6) | ||

\displaystyle{V} | \displaystyle\!\!= | \displaystyle\!\!\sum_{k_{1},\cdots,k_{n}=z,i}\beta_{k_{1}\cdots k_{n}}\sigma_% {k_{1}}^{1}\sigma_{k_{2}}^{2}\cdots\sigma_{k_{n}}^{n}, | (7) | ||

\displaystyle q | \displaystyle\!= | \displaystyle\!\sum_{j=1}^{n}\gamma_{j}\sigma_{z}^{j}, | (8) |

where \sigma_{z}^{j} (j=1,2,\cdots,n) is the Pauli matrix and \sigma_{i}^{j} is the *N*-dimensional identity matrix.
Because our current quantum computing platform can only offer a limited number of qubits
and the focus of this work is on an
implementation of the necessary gate operations under the above encoding, we have employed
a rather aggressive 8-point discretization using n=3 qubits. The associated diagonal forms of the
T, V, and q matrices are
given in the supplementary material. In particular,
the end grid points are at q=\pm 0.8 Å and the locations of other 6 grid points are shown in (Fig. 1B).
The eigenvalues of the ground and first excited states of the bare Hamiltonian
treated in the 8-dimensional encoding Hilbert space are close to the exact answers. The
associated eigenfunctions are somewhat deformed from exact calculations using, e.g., 64 grid points.
Nonetheless, their unbalanced probability distribution in the two potential wells is maintained.
For example, the probability for the first excited state being found in the right potential well is about 80\%.

##### Experiment.

In our experiment qubits 1,2, and 3 are realized by the {}^{19}F, {}^{13}C, and {}^{1}H nuclear spins of Diethyl-fluoromalonate. The structure of Diethyl-fluoromalonate is shown in (Fig. 2A), where the three nuclei used as qubits are marked by oval. The internal Hamiltonian of this system is given by

\displaystyle\mathcal{H}_{int}= | \displaystyle\sum\limits_{j=1}^{3}{2\pi\nu_{j}}I_{z}^{j}+\sum\limits_{j<k,=1}^% {3}{2\pi}J_{jk}I_{z}^{j}I_{z}^{k}, | (9) |

where \nu_{j} is the resonance frequency of the *j*th spin and
\emph{J}_{jk} is the scalar coupling strength between spins *j* and
*k*, with \emph{J}_{12}=-194.4 Hz, \emph{J}_{13}=47.6 Hz, and \emph{J}_{23}=160.7 Hz.
The relaxation time T_{1} and dephasing time T_{2} for each of the three nuclear spins are tabulated in (Fig. 2B). The experiment is conducted on a Bruker Avance 400 MHz spectrometer at room temperature.

The experiment consists of three parts: (A) Initial state preparation. In this part we prepare the ground state \left|\phi_{0}\right\rangle of the bare Hamiltonian T+V as the reactant state; (B) Dynamical evolution, that is, the explicit implementation of the system evolution such that the continuous chemical reaction dynamics can be simulated; (C) Measurement. In this third part the probabilities of the reactant and product states associated with each of the 25 snapshots of the dynamical evolution are recorded. For the jth snapshot at t_{j}\equiv j\delta t, we measure the overlaps C(\left|\psi(t_{j})\right\rangle,\left|\phi_{0}\right\rangle)=|\langle\phi_{0}% |\psi(t_{j})\rangle|^{2} and C(\left|\psi(t_{j})\right\rangle,\left|\phi_{1}\right\rangle)=|\langle\phi_{1}% |\psi(t_{j})\rangle|^{2}, through which the continuous reactant-to-product transformation can be displayed. The main experimental details are as follows. Readers may again refer to the supplementary material for more technical explanations.

(A) *Initial State Preparation*. Starting from the thermal equilibrium state, firstly we create the pseudo-pure state (PPS) \rho_{000}=(1-\epsilon)\mathbb{{I}}/8+\epsilon\left|000\right\rangle\left%
\langle 000\right| using the spatial average technique [26], where \epsilon\approx 10^{-5} represents the polarization of the system and {\mathbb{{I}}} is the 8\times 8 identity matrix.
The initial state \left|\phi_{0}\right\rangle was prepared from \rho_{000} by applying one shaped radio-frequency
(RF) pulse characterized by 1000 frequency
segments and determined by the GRadient Ascent Pulse Engineering (GRAPE) algorithm [27, 28, 29]. The preparation pulse thus obtained
is shown in (Fig. 2C) with the pulse width chosen as 10 ms and a theoretical fidelity 0.995. Because the central resonance frequencies of the nuclear spins are different, (Fig. 2C) shows the RF field amplitudes vs time in three panels.
To confirm the successful preparation of the state |\phi_{0}\rangle, we carry out
a full state tomography and examine the fidelity between the target density matrix \rho_{0}=|\phi_{0}\rangle\langle\phi_{0}| and the experimental one \rho_{exp}(0). Using the
fidelity definition F(\rho_{1},\rho_{2})\equiv\texttt{Tr}(\rho_{1}{\rho_{2}})/\sqrt{(\texttt{Tr}(%
\rho_{1}^{2})\texttt{Tr}(\rho_{2}^{2})},
we obtain F[\rho_{0},\rho_{exp}(0)]=0.950. Indeed, their real parts shown in (Fig. 4A) are seen to be in agreement.

(B) *Dynamical Evolution*. The reaction process was divided into M=25 discrete time intervals of the same duration \delta t.
Associated with the mth time interval, the unitary evolution operator is given by

U_{m}\approx V_{\delta t/2}{E}_{\delta t/2}(t_{m})U_{QFT}T_{\delta t}U_{QFT}^{% \dagger}{E}_{\delta t/2}(t_{m}){V}_{\delta t/2}, | (10) |

where U_{QFT} represents a QFT operation, and other operators are defined by {V}_{\delta t/2}\equiv e^{-\frac{i}{\hbar}{V}\frac{\delta t}{2}}, {T}_{\delta t}\equiv e^{-\frac{i}{\hbar}{T}\delta t}, and {E}_{\delta t/2}(t_{m})\equiv e^{\frac{i}{\hbar}\varepsilon(t_{m-1}+\delta t/2% )eq\frac{\delta t}{2}}, with V, T, and q all in their diagonal representations. Such a loop of operations is m-dependent because the simulated system is subject to a time-dependent laser field. The numerical values of the diagonal operators T_{\delta t}, V_{\delta t/2} and E_{\delta t/2} are elaborated in the supplementary material. A circuit to realize U_{QFT} and a computational network to realize the U_{m} operator are shown in (Fig. 3).

Each individual operation in the U_{m} loop can be implemented by a particular RF pulse sequence applied to our system. However, in the experiment such a direct decomposition of U_{m} requires a very long gate operation time and highly complicated RF pulse sequences. This bottom-up approach hence accumulates considerable experimental errors and also invites serious decoherence effects. To circumvent this technical problem we find a better experimental approach, which further exploits the GRAPE technique to synthesize U_{m} or their products with one single engineered RF pulse only. That is, the quantum evolution operator U(t_{j},0), which is simulated by \prod_{m=1}^{j}U_{m}, is implemented by one GRAPE coherent control pulse altogether, with a preset fidelity and a typical pulse length ranging from 10 ms to 15 ms. For the 25 snapshots of the dynamics, totally 25 GRAPE pulses are worked out, with their fidelities always set to be larger than 0.99. As a result, the technical complexity of the experiment decreases dramatically but the fidelity is maintained at a high level. The task of finding a GRAPE pulse itself may be fulfilled via feedback learning control [19] that can exploit the quantum evolution of our NMR system itself. However, this quantum procedure is not essential or necessary in our experiment because here the GRAPE pulses on a 3-qubit system can be found rather easily.

(C) *Measurement*. To take the snapshots of the reaction process at t_{j}=j\delta t we need to measure the overlaps of C(\left|\psi(t_{j})\right\rangle,\left|\phi_{0}\right\rangle) and C(\left|\psi({t_{j}})\right\rangle,\left|\phi_{1}\right\rangle). A full state tomography at t_{j} will do, but this will produce much more information than needed. Indeed, assisted by a simple diagonalization technique, sole population measurements already suffice to observe the reactant-to-product transformation.

Specifically, in order to obtain C(\left|\psi({t_{j}})\right\rangle,\left|\phi_{0}\right\rangle)=\texttt{Tr}[%
\rho(t_{j})\rho_{0}] with \rho(t_{j})=\left|\psi({t_{j}})\right\rangle\left\langle\psi(t_{j})\right|, we first find a transformation matrix *R* to diagonalize \rho_{0}, that is, \rho_{0}^{\prime}=R\rho_{0}R^{\dagger}, where \rho_{0}^{\prime} is a diagonal density matrix. Letting \rho^{\prime}(t_{j})=R\rho(t_{j})R^{\dagger} and using the identity \texttt{Tr}[\rho(t_{j})\rho_{0}]=\texttt{Tr}[\rho^{\prime}(t_{j})\rho_{0}^{%
\prime}], it becomes clear
that only the diagonal elements or the populations of \rho^{\prime}(t_{j}) are required to measure \texttt{Tr}[\rho^{\prime}(t_{j})\rho_{0}^{\prime}], namely, the overlap C(\left|\psi(t_{j})\right\rangle,\left|\phi_{0}\right\rangle). To obtain \rho^{\prime}(t_{j}) from \rho(t_{j}), we simply add the extra R operation to the quantum gate network. The actual implementation of the R operation can be again mingled with all other gate operations using one GRAPE pulse. A similar procedure is used to measure C(\left|\psi(t_{j})\right\rangle,\left|\phi_{1}\right\rangle).

The populations of \rho{{}^{\prime}}(t_{j}) can be measured by applying [\pi/2]_{y} pulses to the three qubits and then read the ensuing free induction decay signal. In our sample of natural abundance, only \sim 1\% of all the molecules contain a {}^{13}C nuclear spin. The signals from the {}^{1}H and {}^{19}F nuclear spins are hence dominated by those molecules with the {}^{12}C isotope. To overcome this problem we apply SWAP gates to transmit the information of the {}^{1}H and {}^{19}F channels to the {}^{13}C channel and then measure the {}^{13}C qubit.

To assess the difference between theory and experiment, we carry out one full state tomography for the final state density matrix at t=t_{f}. Because the GRAPE pulse is made to reach a fidelity larger than 0.995, the experimental density matrix \rho_{exp}(t_{f}) is indeed very close to the theoretical density matrix \rho_{theory}(t_{f}) obtained in an 8-dimensional Hilbert space, with a fidelity F[\rho_{theory}(t_{f}),\rho_{exp}(t_{f})]=0.957. The experimental density matrix elements of the final state shown in (Fig. 4B) match the theoretical results to a high degree. With confidence in the experimental results on the full density matrix level, we can now examine the simulated reaction dynamics, reporting only the probabilities of the reactant and product states. (Fig. 4C) shows the time-dependence of the probabilities of both the reactant and product states obtained from our quantum simulator. It is seen that the product-to-reactant ratio increases continuously with time, with the probability of the product state reaching 77% at the end of the simulated reaction. At all times, the experimental observations of the reaction process are in impressive agreement with the smooth curves calculated theoretically on a classical computer. Further, the experimental results are also in qualitative agreement with the exact classical calculation using 64 grid points (see Fig. 1B). A prototype laser-driven reaction is thus successfully simulated by our 3-qubit system. We emphasize that due to the use of GRAPE pulses in synthesizing the gate operations, our simulation experiment lasts about 30 ms only, which is much shorter than the spin decoherence time of our system. The slight difference between theory and experiment can be attributed to imperfect GRAPE pulses, as well as inhomogeneity in RF pulses and in the static magnetic field.

##### Conclusion.

Quantum simulation with only tens of qubits can already exceed the capacity of a classical computer. Before realizing general-purpose quantum algorithms that typically require thousands of qubits, a quantum simulator attacking problems not solvable on current classical computers will be one conceivable milestone in the very near future. The realization of quantum simulations will tremendously change the way we explore quantum chemistry in both stationary and dynamical problems [15, 16]. Our work reported here establishes the first experimental study of the quantum simulation of a prototype laser-driven chemical reaction. The feasibility of simulating chemical reaction processes on a rather small quantum computer is hence demonstrated. Our proof-of-principle experiment also realizes a promising map from laser-driven chemical reactions to the dynamics of interacting spin systems under shaped RF fields. This map itself is of significance because it bridges up two research subjects whose characteristic time scales differ by many orders of magnitude.

References and Notes

- 1 M. A. Nielson and I. L. Chuang, Quantum Computation and Quantum Information (Cambrige Univ. Press, Cambridge, U. K., 2000).
- 2 P. Shor, in Proceddings of the 35th Annual Symposium on Foundations of Computer Science (IEEE Computer Society Press, New York, Santa Fe, NM, 1994), p. 124.
- 3 I. Buluta and F. Nori, Science 326, 108 (2009).
- 4 R. P. Feynman, Int. J. Theor. Phys. 21, 467 (1982).
- 5 S. Lloyd, Science 273, 1073 (1996).
- 6 C. Zalka, in ITP Conference on Quantum Coherence and Decoherence (Royal Soc. London, Santa Barbara, California, 1996), pp. 313-322.
- 7 D. S. Abrams and S. Lloyd, Phys. Rev. Lett. 79, 2586 (1997).
- 8 L. A. Wu, M. S. Byrd, and D. A. Lidar, Phys. Rev. Lett. 89, 057904 (2002).
- 9 A. Y. Smirnov, S. Savel¡¯ev, L. G. Mourokh, and F. Nori, Europhys. Lett. 80, 67008 (2007).
- 10 D. A. Lidar and H. Wang, Phys. Rev. E 59, 2429 (1999).
- 11 X. H. Peng, J. F. Du, and D. Suter, Phys. Rev. A 71, 012307 (2005).
- 12 S. Somaroo, C. H. Tseng, T. F. Havel, R. Laflamme, and D. G. Cory, Phys. Rev. Lett. 82, 5381 (1999).
- 13 C. Negrevergne, R. Somma, G. Ortiz, E. Knill, and R. Laflamme, Phys. Rev. A 71, 032344 (2005).
- 14 A. Friedenauer, H. Schmitz, J. T. Glueckert, D. Porras, and T. Schaetz, Nature Phys. 4, 757 (2008).
- 15 A. Aspuru-Guzik, A. D. Dutoi, P. J. Love, and M. Head-Gordon, Science 309, 1704 (2005).
- 16 I. Kassal, S. P. Jordan, P. J. Love, M. Mohseni, and A. Aspuru-Guzik, Proc. Natl. Acad. Sci. USA 105, 18681-18686 (2008).
- 17 J. F. Du, N. Y. Xu, X. H. Peng, P. F. Wang, S. F. Wu, and D. W. Lu, Phys. Rev. Lett. 104, 030502 (2010).
- 18 B. P. Lanyon, J. D. Whitfield, G. G. Gillett, M. E. Goggin, M. P. Almeida, I. Kassal, J. D. Biamonte, M. Mohseni, B. J. Powell, M. Barbieri, A. Aspuru-Guzik, and A. G. White, Nature Chem. 2, 106 (2010).
- 19 H. Rabitz, R. de Vivie-Riedle, M. Motzkus, and K. Kompa, Science 288, 824 (2000); W. S. Warren, H. Rabitz, and M. Dahleh, Science 259, 1581 (1993).
- 20 S. A. Rice and M. Zhao, Optical Control of Molecular Dynamics (John Wiley, New York, 2000); M. Shapiro and P. Brumer, Principles of the Quantum Control of Molecular Processes (John Wiley, New York, 2003).
- 21 D. Wang, J. Chem. Phys. 124, 201105 (2006).
- 22 H. -D. Meyer and G. A. Worth, Theor. Chem. Acc. 109, 251 (2003).
- 23 N. Došlić, O. Kühn, J. Manz, and K. Sundermann, J. Phys. Chem. A 102, 9645 (1998).
- 24 M. D. Feit, J. A. Fleck, and A. Steiger, J. Comput. Phys. 47, 412 (1982).
- 25 J. Z. H. Zhang, Theory and Application of Quantum Molecular Dynamics (World Scientific, Singapore, 1999).
- 26 D. G. Cory, A. F. Fahmy, and T. F. Havel, Proc. Natl. Acad. Sci. USA. 94, 1634 (1997).
- 27 N. Khaneja, T. Reiss, C. Kehlet, T. S. Herbr\ddot{u}ggen, and S. J. Glaser, J. Magn. Reson. 172, 296 (2005).
- 28 J. Baugh, J. Chamilliard, C. M. Chandrashekar, M. Ditty, A. Hubbard, R. Laflamme, M. Laforest, D. Maslov, O. Moussa, C. Negrevergne, M. Silva, S. Simmons, C. A. Ryan, D. G. Cory, J. S. Hodges, and C. Ramanathan, Phys. in Can. 63, No.4 (2007).
- 29 C. A. Ryan, C. Negrevergne, M. Laforest, E. Knill, and R. Laflamme, Phys. Rev. A 78, 012328 (2008).
- 30 Helpful discussions with J. L. Yang are gratefully acknowledged. This work was supported by National Nature Science Foundation of China, the CAS, and the National Fundamental Research Program 2007CB925200.

Fig. 1. Prototype chemical reaction and potential energy curve. (A) Isomerization reaction of nonsymmetric substituted malonaldehydes. (B) Upper panel: Potential energy curve, together with the eigenfunctions of the ground (red) and the first excited (blue) states. The main system parameters [(Eq. 2)] are taken from Ref. [23], with V^{\ddagger}=0.00625\ E_{\rm h}, \Delta=0.000257\ E_{\rm h}, and q_{0}=1\ a_{0}. As a modification, the potential values for q approaching the left and right ends are increased sharply to ensure rapid decay of the wavefunction amplitudes. In particular, this procedure increases the V value at q=\pm 0.8 Å by a factor of 30. The six discrete squares shown on the potential curve and the two end points at q=\pm 0.8 Å constitute the 8 grid points for our 3-qubit encoding. Lower panel: Numerically exact time-dependence of populations of the ground state (reactant state, denoted P{}_{0}) and the first excited state (product state, denoted P{}_{1}).

Fig. 2. (A) Molecular structure of Diethyl-fluoromalonate. The {}^{1}H, {}^{13}C and {}^{19}F nuclear spins marked by oval are used as three qubits. (B) System parameters and important time scales of Diethyl-fluoromalonate. Diagonal elements are the Larmor frequencies (Hz) and off-diagonal elements are scalar coupling strength (Hz) between two nuclear spins. Relaxation and dephasing time scales (second) T_{1} and T_{2} for each nuclear spin are listed on the right. (C) The GRAPE pulse that realizes the initial state \left|\phi_{0}\right\rangle from the PPS \left|000\right\rangle, with a pulse width 10 ms and a fidelity over 0.995. The (blue) solid line represents the pulse power in x-direction, and the (red) dotted line represents the pulse power in y-direction. The three panels from top to bottom represent the RF features at three central frequencies associated with the {}^{19}F, {}^{13}C and {}^{1}H spins, respectively.

Fig. 3. Upper panel: The network of quantum operations to simulate the chemical reaction dynamics, with the reactant state \left|\phi_{0}\right\rangle. The whole process is divided into 25 loops. The operators T_{\delta t}, V_{\delta t} and E_{\delta t/2} are assumed to be in their diagonal representations. Lower panel: H is the Hadamard gate and S, T are phase gates as specified on the right. Vertical lines ending with a solid dot represent controlled phase gates and the vertical line between two crosses represents a SWAP gate.

Fig. 4. Experimental tomography results and the reaction dynamics obtained both theoretically and experimentally. (A)-(B) Real part of the density matrix of the initial and final states of the simulated reaction. Upper panels show the theoretical results based on an 8-dimensional Hilbert space, and lower panels show the experimental results. (C) The measured probabilities of the reactant and product states to give 25 snapshots of the reaction dynamics. The (red) plus symbols represent measured results of C(\left|\psi(t_{j})\right\rangle,\left|\phi_{0}\right\rangle) and the (blue) circles represent measured results of C(\left|\psi(t_{j})\right\rangle,\left|\phi_{1}\right\rangle), both in agreement with the theoretical smooth curves. Results here also agree qualitatively with the numerically exact dynamics shown in (Fig. 1B).

Fig. 1:

Fig. 2:

Fig. 3:

Fig. 4:

## Supporting Online Material

### BACKGROUND ON QUANTUM DYNAMICS SIMULATION

Let us start with the Schrödinger equation:

i\hbar\dot{\psi}(t)=H(t)\psi(t). | (11) |

Its formal solution can be written as

\psi(t)=U(t,t_{0})\psi(t_{0}), | (12) |

where the quantum propagator U(t,t_{0}) is a unitary operator and is given by

U(t,t_{0})={\cal T}\exp\left[-\frac{i}{\hbar}\int_{t_{0}}^{t}\!\!d\tau\,H(\tau% )\right], | (13) |

with {\cal T} being the time ordering operator. There are a number of established numerical methods for propagating the Schrödinger equation, such as Feynman’s path integral formalism [31], the split-operator method [32] and the Chebychev polynomial method [33, 34], etc. For our purpose here we adopt the split-operator method.

The propagator U(t,t_{0}) satisfies

U(t,t_{0})=U(t,t_{N-1})U(t_{N-1},t_{N-2})\cdots\cdots U(t_{1},t_{0})\,, | (14) |

where, for example, the intermediate time points can be equally spaced with t_{m}=m\delta t+t_{0}. For one of such a small time interval, e.g., from t_{m-1} to t_{m}=t_{m-1}+\delta t, we have

\displaystyle U(t_{m},t_{m-1}) | \displaystyle= | \displaystyle{\cal T}\exp\left[-\frac{i}{\hbar}\int_{t_{m-1}}^{t_{m}}\!\!d\tau% \,H(\tau)\right] | (15) | ||

\displaystyle\approx | \displaystyle\exp\left[-\frac{i}{\hbar}\int_{t_{m-1}}^{t_{m}}\!\!d\tau\,H(\tau% )\right], |

where terms of the order (\delta t)^{3} or higher are neglected. For sufficiently small \delta t, the integral in the above equation can be further carried out by a midpoint rule, leading to

\displaystyle U(t_{m},t_{m-1})\approx\exp\left[-\frac{i}{\hbar}H(t_{m-1}+% \delta t/2)\delta t\right]. | (16) |

This integration step has an error of the order of (\delta t)^{2}, which is still acceptable if the total evolution time is not large. Next we separate the total Hamiltonian into two parts:

H(t)=H_{0}(t)+H^{\prime}(t). | (17) |

For example, H_{0}(t) is the kinetic energy part of the total Hamiltonian and H^{\prime}(t) represents the potential energy part. In general H_{0}(t) and H^{\prime}(t) do not commute with each other. The split-operator scheme [32] applied to (Eq. 16) then leads to

{U}(t_{m},t_{m-1})\approx e^{-\frac{i}{\hbar}H^{\prime}(t_{m-1}+\delta t/2)% \delta t/2}e^{-\frac{i}{\hbar}H_{0}(t_{m-1}+\delta t/2)\delta t}e^{-\frac{i}{% \hbar}H^{\prime}(t_{m-1}+\delta t/2)\delta t/2}. | (18) |

The small error of this operator splitting step arises from the nonzero commutator between H_{0}(t) and H^{\prime}(t), which is at least of the order of (\delta t)^{3}. The advantage of the split-operator method is that each step represents a unitary evolution and each exponential in (Eq. 18) can take a diagonal form in either the position or the momentum representation. The error of this operator splitting step is in general smaller than that induced by the aforementioned midpoint rule integration in (Eq. 16). Because in our work the total duration of the simulated chemical reaction is short, the above low-order approach already has a great performance. If long-time simulations with preferably larger time steps are needed in a quantum simulation, then one may use even higher-order split-operator techniques for explicitly time-dependent Hamiltonians [35, 36].

### EXPERIMENTAL IMPLEMENTATION

The experiment consists of three steps: (a) Initial state preparation, which is to prepare the ground state \left|\phi_{0}\right\rangle of the bare Hamiltonian T+{V} (representing the reactant state); (b) 25 discrete steps of dynamical evolution to simulate the actual continuous chemical reaction dynamics; (c) Measurement of the overlaps C(\left|\psi(t_{j})\right\rangle,\left|\phi_{0}\right\rangle)=|\langle\phi_{0}% |\psi(t_{j})\rangle|^{2} and C(\left|\psi(t_{j})\right\rangle,\left|\phi_{1}\right\rangle)=|\langle\phi_{1}% |\psi(t_{j})\rangle|^{2} at t_{j}=j\delta t, which is to show the transformation between the reactant and product states.

#### A. Initial State Preparation

To prepare the ground state \left|\phi_{0}\right\rangle, firstly we need to create a pseudo-pure state (PPS) from the thermal equilibrium state - a mixed state which is not yet ready for quantum computation purposes. The thermal equilibrium state of our sample can be written as \rho_{ther}=\sum\limits_{i=1}^{3}\gamma_{i}I_{z}^{i}, where \gamma_{i} is the gyromagnetic ratio of the nuclear spins. Typically, \gamma_{\texttt{C}}=1, \gamma_{\texttt{H}}=4 and \gamma_{\texttt{F}}=3.7, with a constant factor ignored. We then use the spatial average technique [37] to prepare the PPS

\rho_{000}=\frac{1-\epsilon}{8}\mathbb{{I}}+\epsilon\left|000\right\rangle% \left\langle 000\right|, | (19) |

where \epsilon\approx 10^{-5} represents the polarization of the system and {\mathbb{{I}}} is the 8\times 8 unity matrix. The unity matrix has no influence on our experimental measurements and hence can be dropped. The pulse sequence to prepare the PPS from the thermal equilibrium state is shown in (Fig. S1A). In particular, the gradient pulses (represented by G{}_{\text{Z}}) destroys the coherence induced by the rotating pulses and free evolutions. After obtaining the PPS, we apply one shaped pulse calculated by the GRadient Ascent Pulse Engineering (GRAPE) algorithm [38, 39, 40] to obtain the initial state \left|\phi_{0}\right\rangle, with the pulse width 10 ms and a fidelity 0.995. In order to assess the accuracy of the experimental preparation of the initial state, a full state tomography [41] is implemented. The fidelity [42] between the target density matrix \rho_{target} and the experimental density matrix \rho_{exp} is found to be

\displaystyle F(\rho_{target},\rho_{exp}) | \displaystyle\equiv | \displaystyle\texttt{Tr}(\rho_{target}\rho_{exp})/\sqrt{(\texttt{Tr}(\rho_{% target}^{2})\texttt{Tr}(\rho_{exp}^{2})} | (20) | ||

\displaystyle\approx | \displaystyle 0.95. |

A detailed comparison between \rho_{target} and \rho_{exp} is displayed in (Fig. S1B).

#### B. Dynamical Evolution

To observe the continuous reactant-to-product transformation, we divide the whole time evolution into 25 discrete steps. For convenience all variables here are expressed in terms of atomic units. For example, in atomic units e=1, \hbar=1, and \delta t=62.02. To exploit the split-operator scheme in (Eq. 18), we let H_{0}=T, and H^{\prime}(t)=V-eq\epsilon(t). The kinetic energy operator T is diagonal in the momentum representation, whereas the V operator and the dipole-field interaction -eq\epsilon(t) operator are both diagonal in the position representation. We then obtain from (Eq. 18)

\displaystyle U(t_{m},t_{m-1})\approx V_{\frac{\delta t}{2}}E_{\frac{\delta t}% {2}}U_{QFT}T_{\delta t}U_{QFT}^{\dagger}E_{\frac{\delta t}{2}}V_{\frac{\delta t% }{2}}, | (21) |

where the operators

\displaystyle V_{\frac{\delta t}{2}} | \displaystyle\equiv | \displaystyle e^{-i{V}\frac{\delta t}{2}}, | (22) | ||

\displaystyle T_{\delta t} | \displaystyle\equiv | \displaystyle e^{-i{T}\delta t}, | (23) | ||

\displaystyle E_{\frac{\delta t}{2}} | \displaystyle\equiv | \displaystyle e^{i{q}\varepsilon(t_{m-1}+\delta t/2)\frac{\delta t}{2}} | (24) |

are assumed to be in their diagonal representations.

To map U(t_{m},t_{m-1}) to our 3-qubit NMR quantum computer, we discretize the potential energy curve using 8 grid points. Upon this discretization, operators V_{\frac{\delta t}{2}}, T_{\delta t}, and q become 8\times 8 diagonal matrices. Numerically, their diagonal elements (denoted V_{diag}, T_{diag}, and q_{diag}, respectively) are found to be

\displaystyle{V}_{diag}= | \displaystyle(293.78,-0.10,1.85,5.41, | ||||

\displaystyle 5.46,2.02,0.18,305.44)\times 10^{-3}; | |||||

\displaystyle{T}_{diag}= | \displaystyle(0,0.91,3.63,8.16, | ||||

\displaystyle 14.51,8.16,3.63,-0.91)\times 10^{-3}; | |||||

\displaystyle{q}_{diag}= | \displaystyle(-1.51,-1.08,-0.65,-0.22, | (25) | |||

\displaystyle 0.22,0.65,1.08,1.51). |

To evaluate the E_{\frac{\delta t}{2}} operator, we also need to discretize the time-dependence of the electric field associated with the ultrashort laser pulse. For 25 snapshots of the reaction dynamics, we discretize the trapezoid-type electric field by 25 points, i.e.,

\varepsilon(t)=[0.05,0.42,0.85,1,...1,0.85,0.42,0.05]\times 10^{-3}. | (26) |

The quantum gate network for the QFT operation that transforms the momentum representation to the coordinate representation is already shown in (Fig. 3). It consists of three Hardmard gates (H), three controlled-phase gates (S and T) and one SWAP gate (vertical line linking crosses). The Hardmard gate H is represented by the Hardmard matrix

\displaystyle H=\frac{1}{\sqrt{2}}\left(\begin{array}[]{cc}1&1\\ 1&-1\\ \end{array}\right), | (27) |

which maps the basis state \left|0\right\rangle to \frac{1}{\sqrt{2}}(\left|0\right\rangle+\left|1\right\rangle) and \left|1\right\rangle to \frac{1}{\sqrt{2}}(\left|0\right\rangle-\left|1\right\rangle). The phase gates S and T are given by

\displaystyle\text{S}=\left(\begin{array}[]{cc}1&0\\ 0&i\\ \end{array}\right) | (28) |

and

\displaystyle\text{T}=\left(\begin{array}[]{cc}1&0\\ 0&e^{i\pi/4}\\ \end{array}\right). | (29) |

The matrix form of the SWAP gate is

\displaystyle\text{SWAP}=\left(\begin{array}[]{cccc}1&0&0&0\\ 0&0&1&0\\ 0&1&0&0\\ 0&0&0&1\\ \end{array}\right), | (30) |

which exchanges the state of the {}^{19}F qubit and with that of the {}^{1}H qubit.

GRAPE Pulses. Since there are hundreds of logical gates in the required network of quantum operations, a direct implementation of the gate operation network will need a large number of single-qubit rotations as well as many free evolutions during the single-qubit operations. This bottom-up approach will then accumulate the errors in every single-qubit operation. Considerable decoherence effects will also emerge during the long process. For example, we have attempted to directly decompose the network into a sequence of RF pulses, finding that the required free evolution time for the 25 loops of evolution is more than 1 s, which is comparable to the T_{2} time of our system. To overcome these problems and to reach a high-fidelity quantum coherent control over the three interacting qubits, the unitary operators used in our experiment are realized by shaped quantum control pulses found by the GRadient Ascent Pulse Engineering (GRAPE) technique [38, 39, 40]. To maximize the fidelity of the experimental propagator as compared with the ideal gate operations, we use a mean gate fidelity by averaging over a weighted distribution of RF field strengths to minimize the inhomogeneity effect of the RF pulses applied to the sample.

For a known or desired unitary operator U_{target}, the goal of the GRAPE algorithm is to find a shaped pulse within a given duration t_{total} to maximize the fidelity

F=|\texttt{Tr}(U_{target}^{\dagger}U_{cal})/2^{n}|^{2}, | (31) |

where U_{cal} is the unitary operator actually realized by the shaped pulse and 2^{n} is the dimension of the Hilbert space. We discretize the evolution time t_{total} into N segments of equal duration \Delta t=t_{total}/N, such that U_{cal}=U_{N}\cdots U_{2}U_{1}, with the evolution operator associated with the *j*th time interval given by

U_{j}=\exp\{-i\Delta t(\mathcal{H}_{int}+\sum_{k=1}^{m}u_{k}(j)\mathcal{H}_{k}% )\}. | (32) |

Here \mathcal{H}_{int} is the three-qubit self-Hamiltonian in the absence of any control field, \mathcal{H}_{k} represents the interaction Hamiltonians due to the applied RF field, and u_{k}(j) is the control vectors associated with \mathcal{H}_{k}. Specifically, in our experiment u_{k}(j) are the time-dependent amplitudes of the RF field along the *x* and *y* directions, for the F-channel, the C-channel and the H-channel. With an initial guess for the pulse shape, we use the GRAPE algorithm to optimize u_{k}(j) iteratively until U_{cal} becomes very close to U_{target}. More details can be found from Ref. [38]. The GRAPE technique dramatically decreases the duration and complexity of our experiment and at the same time increases the quantum control fidelity. In our proof-of-principle demonstration of the feasibility of the quantum simulation of a chemical reaction, the task of searching for the GRAPE pulses is carried out on a classical computer in a rather straightforward manner. It is important to note that this technique can be scaled up for many-qubit systems, because the quantum evolution of the system itself can be exploited in finding the high-fidelity coherent control pulses.

As an example, (Fig. S2) shows the details of one 15 ms GRAPE pulse to realize the quantum evolution from t=0 to t_{7}=7\delta t (also combining the operations for initial state preparation and the extra operation R that is useful for the measurement stage). The shown GRAPE pulse is found by optimizing the frequency spectrum divided into 750 segments. The shown GRAPE pulse has a fidelity over 0.99.

#### C. Measurement

To simulate the process of a chemical reaction, it is necessary to measure the simulated reactant-to-product transformation at different times. To that end we measure the overlaps of C(\left|\psi({t_{m}})\right\rangle,\left|\phi_{0}\right\rangle) and C(\left|\psi(t_{m})\right\rangle,\left|\phi_{1}\right\rangle) at t_{m}=m\delta t. Here we first provide more explanations about how a diagonalization method can reduce the measurement of the overlaps to population measurements.

Without loss of generality, we consider the measurement C(\left|\psi_{7}\right\rangle,\left|\phi_{0}\right\rangle).

C(\left|\psi(t_{7})\right\rangle,\left|\phi_{0}\right\rangle)=|\langle\phi_{0}% |\psi(t_{7})\rangle|^{2}=\texttt{Tr}[\rho(t_{7})\rho_{0}], | (33) |

where \rho(t_{7})=\left|\psi(t_{7}\right\rangle\left\langle\psi(t_{7})\right| and \rho_{0}=\left|\phi_{0}\right\rangle\left\langle\phi_{0}\right|. Let *R* be a transformation matrix which diagonalizes \rho_{0} to a diagonal density matrix \rho_{0}^{\prime}=R\rho_{0}R^{\dagger}. Then

\displaystyle\texttt{Tr}[\rho(t_{7})\rho_{0}]=\texttt{Tr}[R\rho(t_{7})R^{% \dagger}R\rho_{0}R^{\dagger}]=\texttt{Tr}[\rho^{\prime}(t_{7})\rho_{0}^{\prime% }], | (34) |

where \rho^{\prime}(t_{7})=R\rho(t_{7})R^{\dagger}.
Clearly then, only the diagonal terms (populations) of \rho^{\prime}(t_{7}) are relevant when calculating \texttt{Tr}[\rho^{\prime}(t_{7})\rho_{0}^{\prime}],
namely, the overlap C(\left|\psi(t_{7}\right\rangle,\left|\phi_{0}\right\rangle). Hence only population measurement of the density matrix \rho^{\prime}(t_{7}) is needed to obtain the overlap between |\psi(t_{7})\rangle and the initial state.
The GRAPE pulse that combines the operations for initial state preparation, for the quantum evolution, as well as for the extra operation *R* is shown in (Fig. S2).

The three population-readout spectra after applying the GRAPE pulse
are shown in (Fig. S3B-S3D), together with the {}^{13}C spectrum for the PPS |000\rangle.
The populations to be measured are converted to the observable coherent terms by applying a [\pi/2]_{y} pulse to each of the three qubits. For the {}^{13}C spectrum shown in (Fig. S3C), four peaks from left to right are seen, with their respective integration results representing P(5)-P(7), P(6)-P(8), P(1)-P(3), and P(2)-P(4), where P(i) is the *i*th diagonal element of \rho^{\prime}(t_{7}). Experimentally the four integrals associated with the four peaks in (Fig. S3C) are found to be -0.098, -0.482, -0.089 and -0.071, which are close to the theoretical values -0.047, -0.501, -0.114 and -0.041.
Further using other readouts from the {}^{19}F (see Fig. S3B) and {}^{1}H (see Fig. S3D) spectra as well as the normalization condition \sum_{i=1}^{8}{P}({i})=1, we obtain all the 8 populations and hence the overlap C(\left|\psi(t_{7}\right\rangle,\left|\phi_{0}\right\rangle). The theoretical and experimental results for this overlap are 0.535 and 0.529, which are in agreement.
A similar procedure is used to obtain C(\left|\psi(t_{m}\right\rangle,\left|\phi_{1}\right\rangle).

The spectra of the {}^{1}H and {}^{19}F channel are obtained by first transmitting the signals of the {}^{19}F and {}^{1}H qubits to the {}^{13}C qubit using SWAP gates. With this procedure all the spectra shown in (Fig. S3) are exhibited on the {}^{13}C channel. Indeed, because in our sample of natural abundance, only \approx 1\% of all the molecules contain a {}^{13}C nuclear spin, the signals from the {}^{1}H and {}^{19}F nuclear spins without applying SWAP gates would be dominated by those molecules with the {}^{12}C isotope.

References

- 31 R. P. Feynman, Rev. Mod. Phys. 20, 367 (1948).
- 32 M. D. Feit, J. A. Fleck, and A. Steiger, J. Comput. Phys. 47, 412 (1982).
- 33 C. Leforestier, R. H. Bisseling, C. Cerjan, M. D. Feit, R. Friesner, A. Guldberg, A. Hammerich, G. Jolicard, W. Karrlein, H. -D. Meyer, N. Lipkin, O. Roncero, and R. Kosloff, J. Comput. Phys. 94, 59 (1991).
- 34 J. Z. H. Zhang, Theory and Application of Quantum Molecular Dynamics (World Scientific, Singapore, 1999).
- 35 A. D. Bandrauk and H. Chen, Can. J. Chem. 70, 555 (1992).
- 36 W. S. Zhu and X. S. Zhao, J. Chem. Phys. 105, 9536 (1996).
- 37 D. G. Cory, A. F. Fahmy and T. F. Havel, Proc. Natl. Acad. Sci. USA. 94, 1634 (1997).
- 38 N. Khaneja, T. Reiss, C. Kehlet, T. S. Herbruggen, and S. J. Glaser, J. Magn. Reson. 172, 296 (2005).
- 39 J. Baugh, J. Chamilliard, C. M. Chandrashekar, M. Ditty, A. Hubbard, R. Laflamme, M. Laforest, D. Maslov, O. Moussa, C. Negrevergne, M. Silva, S. Simmons, C. A. Ryan, D. G. Cory, J. S. Hodges, and C. Ramanathan, Phys. in Can. 63, No.4 (2007).
- 40 C. A. Ryan, C. Negrevergne, M. Laforest, E. Knill, and R. Laflamme, Phys. Rev. A 78, 012328 (2008).
- 41 J. S. Lee, Phys. Lett. A 305, 349 (2002).
- 42 N. Boulant, E. M. Fortunato, M. A. Pravia, G. Teklemariam, D. G. Cory, and T. F. Havel, Phys. Rev. A 65, 024302 (2002).

Fig. S1. Pulse sequence for the preparation of the PPS and a comparison between experimental and theoretical density matrix elements for the reactant state \left|\phi_{0}\right\rangle. (A) Pulse sequence that implements the PPS preparation, with \theta=0.64\pi and \texttt{X}(\overline{\texttt{X}},\texttt{Y},\overline{\texttt{Y}}) representing rotations around the x(-x, y, -y) direction. G{}_{\text{Z}} represents a gradient pulse to destroy the coherence induced by the rotating pulses and free evolutions. \frac{1}{4\text{J}_{\text{HF}}} and \frac{1}{4\text{J}_{\text{CF}}} represent the free evolution of the system under \mathcal{H}_{int} for 5.252 ms and 1.286 ms, respectively. (B) Comparison between measured density matrix elements of the initial state \left|\phi_{0}\right\rangle and the theoretial target density matrix elements based on the 8-point encoding. Both the real part and the imaginary part of the density matrix elements are shown.

Fig. S2. GRAPE pulse to simulate the quantum evolution of the reacting system from t=0 to t_{7}=7\delta t. The top, middle and bottom panels depict the time-dependence of the RF pulses applied to the F-channel, C-channel and H-channel, respectively. The (blue) solid line represents the pulse power applied in the {x}-direction, and the (red) dotted line represents the pulse power applied in the {y}-direction.

Fig. S3. Measured spectra to extract the populations of the system density matrix before or after applying the GRAPE shown in (Fig. S2). (A) {}^{13}C spectrum of the PPS \left|000\right\rangle as a result of a [\pi/2]_{y} pulse applied to the {}^{13}C qubit. The area (integral) of the absorption peak can be regarded as one benchmark in NMR realizations of quantum computation. (B)-(D) Signals from the {}^{19}F, {}^{13}C and {}^{1}H qubits after applying the GRAPE pulse and a [\pi/2]_{y} pulse to each of the three qubits. All the spectra are exhibited on the {}^{13}C channel through SWAP gates. The integration of each spectral peak gives the difference of two particular diagonal elements of the density matrix \rho^{\prime}(t_{7}).

Fig. S1:

Fig. S2:

Fig. S3: