## Abstract

Quantum computers have been proved to be able to mimic quantum systems efficiently in polynomial time. Quantum chemistry problems, such as static molecular energy calculations and dynamical chemical reaction simulations, become very intractable on classical computers with scaling up of the system. Therefore, quantum simulation is a feasible and effective approach to tackle quantum chemistry problems. Proof-of-principle experiments have been implemented on the calculation of the hydrogen molecular energies and one-dimensional chemical isomerization reaction dynamics using nuclear magnetic resonance systems. We conclude that quantum simulation will surpass classical computers for quantum chemistry in the near future.

## 1. Introduction

What type of computer can we use to simulate quantum systems? Certainly the classical computer is not the right choice. Although it can deal with tremendous common problems at present, when the target system is quantum mechanical, whose evolution obeys the Schrödinger equation, the classical computer unfortunately fails. Suppose that the system size scales up and expands rapidly. Then the computational cost grows exponentially and will soon exceed the capacity of any current computers. In 1982, Richard Feynman presented a path to solve this problem [1], p. 468: ‘… there is to be an *exact* simulation, that the computer will do *exactly* the same as nature’. The idea of quantum simulation, using controllable quantum systems to mimic other quantum systems, was suggested for the first time. Lloyd proved that universal quantum simulators can be built with quantum computation architecture, whose required resources scale polynomially with the size of the simulated system [2].

Compared with quantum algorithms, which, in general, require thousands of qubits to surpass classical computing (e.g. Shor's quantum factorization algorithm [3]), a quantum simulator could be realized with merely tens of qubits [4]. Moreover, quantum simulation is less demanding than quantum computing because, typically, it does not require either explicit quantum gates or error corrections. Furthermore, the level of coherent control techniques for quantum systems necessary for the physical realization of quantum simulation is now within reach [5]. Therefore, quantum simulation is attracting plentiful interest in an enormous range of fields, including quantum chemistry, materials science, quantum many-body problems, condensed matter physics, etc. During the past few years, several quantum simulation algorithm designs have been proposed, which were engineered for particular problems [6–10]. Pioneering experiments were first implemented in nuclear magnetic resonance (NMR) systems. The first simulation was on quantum oscillators [11], followed by quantum phase transition [12], quantum many-body problems [13] and a pairing Hamiltonian [14,15]. Later, simulation of quantum phase transitions of a two-spin quantum magnet [16] and simulation of the Dirac equation [17] were performed using trapped ions. In this paper, our discussion is focused on the issue of quantum simulation of quantum chemistry [18], essentially including experiments on static molecular energies and dynamic chemical reactions.

## 2. Simulating molecular ground-state energy

In quantum chemistry, one of the most significant problems is how to calculate molecular energies, which is a fundamental problem in computational quantum chemistry. The resources required on classical computers for a full simulation of the molecular system increase exponentially with increasing number of atoms involved, limiting such full configuration interaction (FCI) calculations of molecular energies to simple diatomic and triatomic molecules [19]. This limitation has led to the development of several approximate methods to predict the chemical properties of large systems, such as density functional theory (DFT) calculations. Although these approximate methods can be superior in solving some problems, there are weaknesses associated with them. For instance, the DFT approach would fail for systems involving long-range charge-transfer excited states [20]. Thus, quantum chemistry desires a universal method—perhaps quantum simulation would be a good choice.

The problem of solving the molecular ground-state energy can equally well be considered as how to calculate the lowest eigenvalue of the Schrödinger equation, which contains a time-independent Hamiltonian. Although classical computers can deal with certain molecules to a fixed accuracy, the desired resources would scale exponentially with the size of molecules and soon exceed the capacities of classical computers. In contrast, theoretically, quantum computers are able to simulate the ground-state energy to the precision required by quantum chemistry with just polynomial time complexity. In 2005, Aspuru-Guzik *et al.* proposed an algorithm to simulate a molecular system and to calculate the energies of molecules [21]. Information about the energy is coded in the phase shift of a quantum register, measured with a quantum phase estimation algorithm (PEA) [7,22]. The main process can be summarized into three steps (figure 1*a*).

—

*Encoding or preparation*. This involves mapping the molecular wave function to the status of qubits. This initializing manipulation is the first step of almost all kinds of quantum simulation tasks, establishing a bridge between the physical system and the simulated system. Undoubtedly, as the qubits needed increase with the size of the target molecule, the method of preparing qubits to the ground state*ϕ*of the molecular Hamiltonian will become more complicated.—

*Evolution.*This involves simulating the time evolution of the molecular Hamiltonian. The goal of this step is to generate a phase shift on the probe qubit, which reflects the information of the eigenvalue of the ground state energy. The PEA indicates that, if we apply a unitary operator*U*=e^{−iHτ}on the ground state |*ψ*〉, a phase will be generated as*U*|*ψ*〉=e^{−iHτ}|*ψ*〉=e^{−iEτ}|*ψ*〉=e^{i2πϕ}|*ψ*〉, where*E*=−2*πϕ*/*τ*is the energy of the ground state. Therefore, measurement of the phase is equivalent to calculating the ground-state energy.—

*Measurement.*This involves measuring the phase shift of the probe qubit, in order to get the value of the ground-state energy. To fulfil this purpose, a four-bit inverse quantum Fourier transform (QFT) is adopted for the relative phase measurement, needing four qubits to obtain one precise bit with a probability of success of .

In the theoretical proposal, it is demonstrated that the number of qubits required scales linearly with the number of basis functions, and the number of gates required grows polynomially with the number of qubits. Additionally, calculations of the water and lithium hydride molecular ground-state energies were carried out on a quantum simulator using a recursive PEA.

Although this proposal is almost perfect in theory, examining such large molecules, i.e. water and lithium hydride, is infeasible at the time of the experiment. When the system expands to a large Hilbert space, it is difficult for the current quantum computing technology to support requirements such as the initial state preparation, the decomposition of the circuit and the readout of the phase shift. As the first step of quantum simulation for quantum chemistry, both Lanyon *et al.* [23] and the present authors [24] selected the simplest but significant scenario: simulating the ground-state energy of the hydrogen molecule. In Lanyon's paper, a key step to the realization of the algorithm was carried out on a photonic system, and the ground-state energy was calculated with 20-bit precision. In our experiment, we chose a widely used STO-3G basis set in FCI [25], and obtained the ground-state energy to 45-bit precision using an NMR quantum simulator.

From the theory, adopting the Born–Oppenheimer approximation, the Hamiltonian of the electron in the hydrogen molecule can be described as
2.1where *T*_{i} is the kinetic energy of the *i*th electron, *V* _{ij} is the Coulomb potential energy between the *i*th electron and the *j*th nucleus, and *O*_{ij} is the Coulomb potential energy between the *i*th electron and *j*th electron. In this model, each atom has a 1s Gaussian-type function in STO-3G basis, consisting of one bonding orbital with *gerade* symmetry and one antibonding orbital with *ungerade* symmetry. This structure forms four spin orbits corresponding to six possible configurations. Furthermore, only two configurations contribute to the ground-state energy calculation: the ground-state configuration |*Ψ*_{0}〉 and the double excitation configuration . So the singlet symmetry and spatial symmetry of the ground state can be simplified. In atomic units (a.u.), the nuclear distance *γ* is 1.4 a.u., and the Hamiltonian of this model is [26]
2.2whose eigenvalue is −1.851 570 929 351 19 a.u.

To perform this simulation, we selected a two-qubit NMR quantum simulator. Each qubit is represented by ^{13}C (system qubit) and ^{1}H (probe qubit) nuclear spins in ^{13}C-labelled chloroform dissolved in *d*_{6}-acetone. The molecular structure is shown in figure 1*b*, where the two qubits are marked. The internal Hamiltonian of this two-qubit system is
2.3where *ω*_{H}/2*π* and *ω*_{C}2*π* are the Larmor frequencies and *J*_{CH} represents the weak coupling constant, typically, *J*_{CH}=214.6 Hz.

The whole experiment can be divided into three sections: (*a*) the adiabatic preparation of the system qubit ^{13}C to the ground state of the Hamiltonian; (*b*) the application of the time evolution of the Hamiltonian on the qubits to generate the phase shift on the probe qubit; and (*c*) measurement of the phase shift on the qubit to extract the energy information. Now, we will discuss these three parts in detail.

### (a) Initial-state preparation

The thermal equilibrium state is a mixed state, which is not suitable for most quantum simulating tasks. Therefore, we first create a pseudo-pure state (PPS) using the spatial average technique [27], with **I** representing the 4×4 unity operator and *ϵ*≈10^{−5} the polarization. Then we prepare the system qubit to the ground state |*Ψ*〉 through the adiabatic state preparation (ASP) process, and the probe qubit to the state by a pseudo-Hadamard gate .

To prepare the ground state |*Ψ*〉, the ASP method based on the quantum adiabatic theorem [28,29] is a simple and accurate approach. Besides, the adiabatic computational model has been proposed [30] and has been proved equivalent to the conventional circuit model [31]. The quantum adiabatic theorem states that a quantum system remains in its instantaneous eigenstate if the system Hamiltonian varies slowly enough and if there is a gap between this eigenvalue and the rest of the Hamiltonian's spectrum [28,29]. Applying this technique, we drove the qubits from the ground state of a simple Hamiltonian *H*_{0}=*σ*_{x} to the target by varying the system Hamiltonian sufficiently slowly. The system qubit is prepared into the ground-state of *H*_{0} by a conjugated pseudo-Hadamard gate . The slowly varying discrete Hamiltonian is interpolated from a linear form *H*_{ad}=(1−*s*)*σ*_{x}+*sH* with *s*=*t*/*T*. The success of ASP is ensured by *T*=5.52 a.u. yielded by the numerical simulation [32,33]. Thus, the unitary operator for each adiabatic step, using Trotter's formula, is
2.4where the duration *δt*=*T*/(*M*+1) and *s*_{m}=*m*/(*M*+1). A simple sequence, including three single-rotating pulses , can implement *U*_{m} directly, and the total evolution can be realized by repeating this sequence consequently. So far we have completed the preparation of the initial state.

### (b) Time evolution

In this step, we will first introduce two techniques, the NMR interferometer [34,35] and the iterative scheme, followed by the implementation of the methods.

In the theoretical proposal mentioned already, a four-bit QFT is essential for evaluating the phase shift. However, in our case, in the NMR system the relative phase was read out by an NMR interferometer, whose principle is very similar to that of an interferometer in optics. This device is a sensitive detector and more precise than the original four-bit QFT gate, as the error bound of the phase shift measured is less than ±5^{°}. As well, the theory ensures that arbitrary precision can be achieved using the iterative process. The initial condition was *U*_{0}=*U* and the iterative form was selected as *U*_{k+1}=[e^{−i2πϕk′}*U*_{k}]^{2n}, where *n* is the number of bits and *ϕ*_{k}′=max{*ϕ*_{k}−*ϕ*_{errbd},0}.

The controlled-rotating gate (represent by the solid dot) functions as follows: when the control qubit ^{1}H is |0〉, the system ^{13}C remains unchanged; and when the control qubit is |1〉, the system undergoes the *U*_{k} evolution. Thus, the operator can be described as
2.5Applying *U*_{0}=e^{−iHτ}, the initial state transformed into
2.6The relative shift can be read out directly in NMR with this interferometer. Here the evolution time is
2.7The pulse sequence for conducting *C*_{Uk} is shown in figure 1*c*.

### (c) Phase measurement

After each iteration, we measured the phase shift and thereby adjusted the operator for the next iteration using the scheme described earlier. The quadrature detection in NMR serves as a phase detector. If the initial state is set as the reference phase, all the relative phases can be read out directly in the NMR spectra. Additionally, the value *ϕ*_{k}′ in each iteration is stored for further use. After measuring all the 15 phases, we use a recursive method to rebuild *ϕ* as the experimental result, which is formulated as
2.8where is merely the intermediate value with no physical meaning. Finally, we get the experimental result, which approaches the theoretical value of *ϕ*. In 15 iterations, the eigenvalue of the molecular Hamiltonian is extracted as −1.851 570 929 351 124, and the intermediate result after each iteration is shown in figure 2.

To conclude, we simulate the hydrogen molecule and calculate its ground-state energy to a precision of 45 bits in our NMR quantum simulator. These two experiments [23,24] are significant progress towards the goal of achieving a rapid quantum chemistry simulation. In terms of scalability, the main issues are the efficient decomposition of the evolution operator [23] and the complexity of the adiabatic passage [36], both of which are at least feasible for medium molecular sizes. It is foreseeable that, with the confronted technical challenges overcome gradually, a quantum simulator is expected to speed up quantum chemistry calculations on a medium-scale quantum computer.

## 3. Simulating chemical isomerization reaction dynamics

We have discussed an implementation of quantum simulation of a static chemical model in the previous section. Now our interests turn to quantum simulation of a dynamic chemical process. Classical computers have done remarkably well at simulating small quantum chemical systems, and, assisted by the Born–Oppenheimer approximation and various techniques [37–39], these methods have reduced the original exponential cost problem to a polynomial one, at given accuracy. However, Kassal *et al.* [40] proposed an algorithm illustrating that quantum chemical reaction simulations can be done universally in polynomial time complexity without approximations. In this algorithm, they introduced ancillary qubits into the procedure, for ease of decomposing operators into elementary quantum gates, when the system scales up to more than a few qubits. Also they iterated the process to improve the accuracy of the quantity of interest. Furthermore, they demonstrated that simulating the entire electronic and nuclear interaction evolving in time is more effective than applying the Born–Oppenheimer approximation for a quantum computer, compared to classical computers, provided that the chemical reactions involve more than four atoms.

Based on their algorithm, we experimentally realized the quantum simulation of a chemical isomerization reaction using an NMR quantum simulator [41]. In particular, we testified a one-dimensional model of the reaction, driven by a laser, of relocating hydrogen on non-symmetric substituted malonaldehydes [42] shown in figure 3*a*. In this model, we used the method called the ‘hydrogen-subway’ with low intensities of laser to avoid the Keldysh limit [43], which is a potential problem in conventional above-barrier pump–dump approaches involving high-intensity laser fields. Both methods achieve the same percentage of product yields (figure 3*b*).

The experiment is carried out on a Bruker Avance 400 MHz spectrometer at room temperature. 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 figure 4*a*, where the three nuclei used as qubits are marked by an oval. The internal Hamiltonian of this system is given by
3.1where *ν*_{j} is the resonance frequency of the *j*th spin and *J*_{jk} is the scalar coupling strength between spins *j* and *k*.

In our experiment, we conducted the process in the following steps: preparing the molecule in the desired initial state, followed by evolving it with time, and then taking measurements during and after the experiment. Starting from the originally time-independent Hamiltonian system operator, we modified it by adding the external laser field into consideration,
3.2where *T* and *V* are, respectively, the kinetic and the potential energy operator and *E*(*t*) is the contributing term from the laser operation. Additionally, *T* and *V* solely depend on positions, and the latter takes the form of a double-well potential of the system. Thus *E*(*t*), the time-varying laser–molecule interaction Hamiltonian, offered the ‘hydrogen-subway’ approach, enabling the reaction to occur while the energy of the molecule remains below the barrier height [42].

Owing to the asymmetry of the two wells, we assumed that the initial reactant state is the ground state of the system and is primarily localized in the left potential well shown in figure 3*b*. The corresponding product state is taken as being chiefly localized in the right potential well, which is expected to be the first excited state of *T*+*V* . For the purpose of preparing our chemical molecule in the initial state, the state under the bare Hamiltonian *T*+*V* , in practice, we created a pseudo-state from the thermal equilibrium state, which was then applied with a pre-set radio-frequency (RF) pulse, determined by the GRadient Ascent Pulse Engineering (GRAPE) algorithm [44–46]. At the end of the procedure, in order to confirm that the molecule was in the reactant state, from a full-state tomography and a fidelity test, we acquired results indicating a strong agreement between the theoretical initial state and the prepared state (figure 5*a*)."

Now we have obtained our initial wave function for the molecule, we would like to evolve it in time. In other words, we applied the propagator of the Hamiltonian operator to the function labelled as *U*_{m}, *U*(*δt*)=e^{−iHδt}, where *H* is defined in equation (3.2). Here we quantized *V* , the potential well, with eight discrete points, corresponding to three qubits. Also, because it is tremendously difficult to implement the propagator operator directly, we used Trotter's formula [47] to decompose the propagator into a simpler form as follows:
3.3where the propagator *U*(*t*+*δt*,*t*) was related to the time interval between *t* and *t*+*δt*. Therefore, we decomposed it into unitary diagonal operators in either the position or momentum representation, which were associated into one coordinate system by executing QFT to one of the representations [48]—in this case, to the momentum one. Thus, *U*_{m} could be emulated in a simple manner, because the decomposed components were all unitary and diagonal in position representation, and one loop was accomplished by completing all these operators in order as shown in figure 3*c*, with an RF pulse sequence.

Indeed, we could obtain the same result by introducing an ancilla, but that would have required extra qubits, which were not available at the time of the experiment. Instead, we turned to the GRAPE approach, effective and successful on the NMR simulator with low dimensions of state functions. However, at a high dimension of the density matrix of the wave function, an ancilla provides a feasible solution to decompose the unitary operator, *U*, into simple additions and other arithmetical operations.

Besides, iterating this procedure would improve the accuracy of the final result. In practice, we repeated the process 25 times. The number of iterations was limited by the decoherence time of the NMR system and the operation time to implement each operator by applying RF pulses. In order to push to the limit of precision, we combined the individual RF pulse for each operator in one loop, into one single RF pulse with high fidelity using GRAPE, which is a facile task on a three-qubit system, and thus reduces the technical complexity of the experiment and the overall operation time, avoiding additional experimental errors and severe decoherence effects. Consequently, we reached highly accurate outcomes without losing fidelity. Between the iterations, we measured the overlap between the density matrix and that of the theoretical initial state and the final state, respectively, to perceive the transition from the reactant to the product. Although the measurements could be achieved with full state tomographies, the information we were interested in was extracted by population measurements, with the aid of a simple diagonalization technique, simplifying the experiment and reducing the required resources.

In addition, we would like to examine the difference between theory and experiment for the final result. By applying a full state tomography to the final-state density matrix at the end of the iterations, the experimental density matrix was found to agree with the theoretical one acquired in an eight-dimensional Hilbert space to a great degree, as shown in figure 5*b*, with a high fidelity of 0.957 out of 1, which resulted from the high fidelity of the GRAPE pulse. Furthermore, the experimental density matrix elements of the final state matched the theoretical results remarkably. Thus, simulated reaction dynamics can be studied by inspecting the probabilities of being in the reactant and product states, respectively, through the reaction, given the confidence on the experimental outcomes on the full density matrix level. Figure 4*b* illustrates the time dependence of the probabilities of being in either state attained from our quantum simulator. From the graph, the product-to-reactant ratio grows continuously with respect to time, with 77 per cent probability of product state when the simulation terminated.

As a result, our three-qubit system positively simulated a prototype laser-driven chemical reaction. Additionally, we point out that the spin decoherence time of our system is significantly longer than our simulation experiment, 30 ms, which is due to combining the gate operations and GRAPE pulses, and thus is due to the reduction of the number of operations. Further, a few sources could contribute to the minor discrepancy between the theory and the experiment, for instance, the unfaithfulness of the GRAPE pulses, and inhomogeneity in the RF pulses and in the static magnetic field.

## 4. Conclusion

Quantum simulation provides a huge prospect in quantum chemistry in terms of static molecule modelling and dynamic reaction simulation. Our proof-of-principle experiments verified those two points for simple molecules and reactions. However, with only a few tens of qubits, quantum simulation could reach beyond the limit of classical computing [21], i.e. simulating large molecules or complex reactions. This could be achieved when technical barriers are overcome, which occur in mapping eigenstates into qubits and decomposing the molecular evolution operator into quantum logic gates. In NMR systems, one method to encode the initial eigenstate into qubits is via adiabatic state preparation. However, this has problems when the system scales up, as calculating a conclusive mathematical analysis is hard at the moment. Nevertheless, progress has been made in this area, the scalability of the NMR quantum computer [36], and it is foreseeable that quantum simulation will have practical use for medium-sized molecules and reactions.

## Acknowledgements

This work was supported by National Natural Science Foundation of China (grants nos 10834005, 91021005 and 21073171), the CAS and the National Fundamental Research Program 2007CB925200.

## Footnotes

One contribution of 14 to a Theme Issue ‘Quantum information processing in NMR: theory and experiment’.

- This journal is © 2012 The Royal Society