Quantum Simulation of Resonant Transitions for Solving the Eigen-problem of an Effective Water Hamiltonian

TL;DR

This paper demonstrates a quantum algorithm using a four-qubit simulator to efficiently find energy spectra and eigenstates of a water molecule's effective Hamiltonian, overcoming classical computational limitations.

Contribution

The work introduces an experimental quantum algorithm leveraging resonant transitions to solve eigen-problems of molecular Hamiltonians with a small quantum simulator.

Findings

Successfully obtained the energy spectrum of a water molecule's effective Hamiltonian.

Demonstrated state transformation accessing large Hilbert space regions.

Efficiently prepared specific eigenstates based on measured eigenenergies.

Abstract

It is difficult to calculate the energy levels and eigenstates of a large physical system on a classical computer because of the exponentially growing size of the Hilbert space. In this work, we experimentally demonstrate a quantum algorithm which could solve this problem via simulated resonant transitions. Using a four-qubit quantum simulator in which two qubits are used as ancillas for control and measurement, we obtain the energy spectrum of a 2-qubit low-energy effective Hamiltonian of the water molecule. The simulated transitions allow the state of the quantum simulator to transform and access large regions of the Hilbert space, including states that have no overlap with the initial state. Furthermore, we make use of this algorithm to efficiently prepare specific eigenstates on the simulator according to the measured eigenenergies.

Figures

Equations

$$\mathcal{H}=-\frac{\omega }{2}{\sigma }_z\otimes {\mathbf{I}}_2^{\otimes \left(n+1\right)}+{\mathbf{I}}_2\otimes {\mathcal{H}}_R+c{\sigma }_x\otimes A,$$

$${\mathcal{H}}_R=|0⟩⟨0|\otimes {\mathcal{H}}_{E_0}+|1⟩⟨1|\otimes {\mathcal{H}}_S.$$

$$P_j(\omega )=\frac{(2c|d_j|{)}^2}{{\Omega }_{0j}^2}{\sin }^2\left(\frac{{\Omega }_{0j}\tau }{2}\right),j=1,2,...,N,$$

$$|\Psi ⟩=$$

$$+e^{i{\varphi }_1}\sqrt{P_j({\omega }_j)}|0⟩|1⟩|E_j⟩,$$

$$\mathcal{H}=-$$

$$+{\mathbf{I}}_2\otimes \left(|0⟩⟨0|\otimes E_0|00⟩⟨00|+|1⟩⟨1|\otimes {\mathcal{H}}_S\right),$$

$$H_{\mathrm{NMR}}=\sum _{j=1}^4\pi {\nu }_j{\sigma }_z^j+\sum _{1⩽j<k⩽4}\frac{\pi }{2}{J}_{jk}{\sigma }_z^j{\sigma }_z^k,$$

$${\mathcal{H}_{S}}=\left({\begin{array}[]{*{20}{c}}{-83.9566}&{-0.0820}&{0.0458}&{0.0594}\\ {-0.0820}&{-83.4080}&{0.0110}&{0.0767}\\ {0.0458}&{0.0110}&{-82.5661}&{0.1323}\\ {0.0594}&{0.0767}&{0.1323}&{-82.4800}\end{array}}\right)$$

Full text

Calculating molecular energy eigenvalues and eigenstates is an essential task in quantum physics and chemistry. Normally, one has to solve the time-independent Schrödinger equation of the system for this purpose. However, numerical simulation of the Schrödinger equation on a classical computer is inefficient, because the computational cost grows exponentially with the size of the simulated system Levine (2000); Thøgersen and Olsen (2004). In quantum chemistry, computational algorithms are usually based on approximate methods that provide rough estimates for chemical properties, and their accuracy varies with the nature of the chemical elements involved. Most of these mathematical algorithms, such as the Variation method Szabo and Ostlund (1996) and Davidson’s algorithm Ernest (1975), rely on the diagonalization of the Hamiltonian matrix in an iterative manner with a proper trial state. For large systems, especially with complex interactions, trial states that sufficiently overlap with the desired eigenstates are quite difficult to predict, which generally leads to slow convergence in the iterative computation.

A quantum computer could efficiently simulate a quantum system by using resources that scale only polynomially with the size of the system Feynman (1982). In recent years, quantum simulation of a variety of physical properties has been studied and demonstrated in several experiments Georgescu et al. (2014); O’Malley et al. (2016); Paesani et al. (2017); Kandala et al. (2017); Du et al. (2010); Lanyon et al. (2010, 2011); Kandala et al. (2018). For instance, the quantum phase estimation algorithm (PEA), which can be used to obtain some eigenvalues of a Hamiltonian with exponential speedup over classical algorithms Aspuru-Guzik et al. (2005); Dobšíček et al. (2007), has been implemented to determine the ground-state energy of the Hydrogen molecule Du et al. (2010); Lanyon et al. (2010). The success probability of the PEA depends on the overlap between the trial state and the desired energy eigenstate. The initial guess of the wave function is usually based on either polynomial-scaling classical $ab~{}initio$ methods Veis and Pittner (2010); Wang et al. (2008) or the output of the adiabatic state preparation (ASP) method Aspuru-Guzik et al. (2005). For complicated or excited many-body systems, it is very difficult to guess and efficiently prepare a good trial state for almost any desired eigenstate on a quantum computer.

Here, we report a nuclear magnetic resonance (NMR) demonstration of a quantum algorithm for determining the energy levels and eigenstates of a physical Hamiltonian in a given energy range via the simulation of controlled resonant transitions Wang et al. (2012); Berkley et al. (2013); Babbush et al. (2014); Wang (2016). The algorithm does not require preparing an approximation for any particular eigenstate of the simulated Hamiltonian. Instead we start the experiments from a simple reference state and simulated the resonance transitions from this initial state to the eigenstates of the physical Hamiltonian of interest. In the experiment we determine the eigenenergies of the water molecule and then use this information to prepare the corresponding eigenstates on the quantum simulator.

The basic idea for applying the algorithm to a physical system is schematically shown in Fig. 1. The quantum simulator comprises a probe qubit coupled to a $(n+1)$-qubit quantum register $R$. The register $R$ consists of one ancillary qubit and $n$ system qubits that encode the simulated system. The probe qubit exhibits a dynamical response only when it is resonant with a transition between a pair of energy eigenstates of the register $R$ with the same energy separation. The non-invasive measurement via the probe qubit is crucial for the use of the algorithm for eigenstate preparation and for performing a sequence of iterations through which the state of the simulator can move far away in the Hilbert space from the initial state.

The Hamiltonian of the simulator has the form

[TABLE]

where ${\mathbf{I}_{2}}$ is the two-dimensional identity operator, and $\sigma_{x,y,z}$ are Pauli matrices. The first term in Eq. (1) denotes the Hamiltonian of the probe qubit which has an energy gap $\omega$ ($\hbar=1$). The second term represents the Hamiltonian of register $R$ and is given by

[TABLE]

Here $\mathcal{H}_{E_{0}}$ is defined as $\mathcal{H}_{E_{0}}=E_{0}\left|E_{0}\right\rangle\left\langle E_{0}\right|$ with $E_{0}$ serving as an energy reference point and the eigenstate $\left|E_{0}\right\rangle$ being the input state. This state does not need to have substantial overlap with any particular target eigenstate. The transitions that take place during the implementation of the algorithm can transform the state of the simulator to states that have no overlap with the initial state. It should also be noted here that if one has a good guess state, e.g. one obtained using coupled-cluster or tensor-network methods, it can be used as an initial state to enhance the efficiency of the algorithm. ${\mathcal{H}_{S}}$ is the Hamiltonian of the simulated system whose Hilbert space size is $N=2^{n}$ and $\left|{{E_{j}}}\right\rangle$ denotes the $j$th eigenstate with eigenenergy $E_{j}$, $j=1,...,N$.

The third term in Eq. (1) describes the interaction between the probe qubit and the register $R$. The coefficient $c$ is the coupling strength, which should be set to a value that is much smaller than $\omega$ to suppress off-resonant transitions. The operator $A$ is defined as $A=\sigma_{x}\otimes B$ where $B$ is the operator that will drive the transformation in the system qubits.

In the initial state, we prepare the probe qubit in its excited state $\left|1\right\rangle$ and the $(n+1)$-qubit register $R$ in the reference state $\left|0E_{0}\right\rangle$. Then, the quantum simulator evolves under the Hamiltonian $\mathcal{H}$ with a given $\omega$ for time $\tau$. When the energy gap between $E_{0}$ and $E_{j}$ equals the gap of the probe qubit, i.e., $\omega=E_{j}-E_{0}$, the simulator is on resonance, which leads to oscillations where the probe qubit is de-excited from $\left|{1}\right\rangle$ to $\left|{0}\right\rangle$ while register $R$ is excited from $\left|0E_{0}\right\rangle$ to $\left|{{1E_{j}}}\right\rangle$. The transition probability $P$ is given by

[TABLE]

where $d_{j}={\left\langle{{1E_{j}}}\right|A\left|{{0E_{0}}}\right\rangle}={\left\langle{{E_{j}}}\right|B\left|{{E_{0}}}\right\rangle}$ and $\Omega_{0j}={\sqrt{(2c|d_{j}|)^{2}+{{\left({{E_{j}}-{E_{0}}-{\omega}}\right)}^{2}}}}$. After the time evolution we measure the probability $P(\omega)$ of the probe qubit to be in state $\left|0\right\rangle$. Absorption peaks appear in $P(\omega)$ at $\omega$ values that are close to $E_{j}-E_{0}$. Hence, the observation of a peak at frequency $\omega$ can be seen as the measurement of one energy eigenvalue of $\mathcal{H}_{S}$ with the value ${E_{j}}={E_{0}}+{\omega}$. It should be noted that the performance of the algorithm can be optimized by choosing appropriate values for the coupling strength $c$, transition operator $B$ and evolution time $\tau$ Wang (2016). For example, one could first run the algorithm with a relatively large value of $c$, which would lead to broad and easily detectable resonances. One could then rerun the algorithm with smaller values of $c$ and focus on the regions where the broad resonances were observed to determine the locations of the resonance peak centers with higher accuracy. It should also be noted that, if one chooses a relatively large value of $c$, one must be careful to choose other parameters, e.g. the Trotter time step, accordingly, e.g. to keep errors from growing beyond the desired measurement accuracy.

By setting $\omega$ at the resonant frequency $\omega_{j}=E_{j}-E_{0}$ and $c\ll\left|{E_{j}}-{E_{i}}\right|$, which suppresses off-resonant transitions, the final state of the quantum simulator is approximately of the form:

[TABLE]

where the phases $\phi_{0}$ and $\phi_{1}$ do not lead to any observable effects in the algorithm. In this situation the transition probability ${P_{j}(\omega_{j})}={\sin^{2}}\left({{c|d_{j}|\tau}}\right)$, and hence resonant Rabi oscillations in the time domain could be observed if one varies the time $\tau$ with fixed $\omega$. By choosing the appropriate $\tau$ that leads to the maximal $P$, which is 1 on resonance, the desired eigenstate $|E_{j}\rangle$ can be prepared deterministically.

In the experiment, we implement the algorithm for the water molecule as an example. Considering the electronic structure of the water molecule and multi-reference-configuration interaction (MRCI) calculations Szabo and Ostlund (1996); Wang et al. (2008, 2012), we encode the low-energy Hamiltonian $\mathcal{H}_{S}$ of the water molecule in a four-dimensional Hilbert space. Details of the Hamiltonian matrix calculation are described in the Supplementary Material Sup . In this article, energies and time are recorded in units of Hartree and Hartree*-1*, respectively. We set $\left|E_{0}\right\rangle=\left|00\right\rangle$. Hence, the Hamiltonian of the quantum simulator in Eq. (1) has the form

[TABLE]

in which the three terms represent the Hamiltonian of the probe, the interaction Hamiltonian and the Hamiltonian of register $R$, respectively. Here, $A={\sigma_{x}}\otimes B={\sigma_{x}}\otimes H_{d}\otimes H_{d}$ acts as the transition operator with the Hadamard operator $H_{d}$, which can transform register $R$ from the reference state $\left|0\right\rangle\otimes\left|00\right\rangle$ to the state ${\left|1\right\rangle}\otimes\frac{1}{2^{2}}\sum\nolimits_{k=0}^{3}{\left|k\right\rangle}$ covering all the computational basis states in the subspace of $\left|1\right\rangle\!\!\left\langle 1\right|\!\otimes\!\mathcal{H}_{S}$.

The experiments were carried out on a Bruker AV-400 NMR spectrometer at room temperature. We used ${}^{13}\rm C$-iodotrifluoroethylene ($\rm{C}_{2}\rm{F}_{3}\rm{I}$) dissolved in $d$-chloroform as the quantum simulator, with a ${}^{13}\rm{C}$ nuclear spin acting as the probe qubit, a ${}^{19}\rm{F}_{1}$ nuclear spin as the ancillary qubit and ${}^{19}\rm{F}_{2}$ and ${}^{19}\rm{F}_{3}$ nuclear spins as qubits for the simulated system. The natural Hamiltonian of this four-nuclear-spin liquid-state NMR system in the doubly rotating frame is

[TABLE]

where $\nu_{j}$ is the chemical shift of spin $j$, and $J_{jk}$ represents the coupling constant between spin $j$ and spin $k$. The molecular structure and Hamiltonian parameters of this sample are shown in Fig. 2(a). The procedures of the experimental implementation are described as follows.

(i) Preparation of the initial state. Starting from the thermal equilibrium state of the four-qubit quantum simulator, we employ the line-selective method to prepare the pseudopure state (PPS) Gershenfeld and Chuang (1997); Peng et al. (2001) $\rho=\left({1-\varepsilon/16}\right){{\mathbf{I}}_{2}^{\otimes 4}}+\varepsilon\left|{0000}\right\rangle\left\langle{0000}\right|$, in which $\varepsilon\approx{10^{-5}}$ represents the polarization. Then, by applying a $\pi$ radio-frequency pulse along the $Y$ axis to the probe qubit, we prepare the simulator in the initial state $\left|\Psi\right\rangle_{0}=\left|{1000}\right\rangle$ where the probe qubit is in the excited state $\left|1\right\rangle$ and register $R$ remains in $\left|{000}\right\rangle$ serving as the reference state.

(ii) Implementation of the evolution operator. We set the reference energy at $E_{0}=-84.20$, the scanning frequency range of the probe qubit as $\left[{{\omega_{\min}},{\omega_{\max}}}\right]=\left[{0,2}\right]$, and the parameters in the Hamiltonian $\mathcal{H}$ at $c=0.006$ and $\tau=1000$. By dividing the frequency range into $100$ intervals of width $\Delta\omega=0.02$, we form a scanning data set for the probe qubit with ${\omega_{k}}={\omega_{\min}}+{k}\Delta\omega$, $k=1,2,...,100$. The quantum circuit to implement the evolution operator $U=e^{-i\mathcal{H}\tau}$ using the Trotter formula Reiher et al. (2017); Casanova et al. (2012); Poulin et al. (2015) is shown in Fig. 2(b) Oritz et al. (2001); Tseng et al. (1999), keeping in mind that some studies have suggested that alternative methods such as sparse-matrix Berry et al. (2007) and Taylor-series methods Berry et al. (2015) might be more efficient for larger systems. In the experiment, the evolution operator for each $k$ is implemented using a shaped pulse which has a duration of 20 ms and is optimized by using the gradient ascent pulse engineering technique Khaneja et al. (2005). All the shaped pulses have theoretical fidelities above 99.5% and are robust against the inhomogeneity of radio-frequency pulses.

(iii) Measurement of the energy spectrum. The final state of the simulator after evolution is approximately given by Eq. (4). When we trace out the subspace of the register $R$, the final density matrix of the probe qubit is reduced to {\rho^{\prime}}=\left({\begin{array}[]{*{20}{c}}P_{j}(\omega)&0\\ 0&{1-P_{j}(\omega)}\end{array}}\right). By applying the readout pulse $Y_{\frac{\pi}{2}}$ to the probe qubit, the ensemble average probability $P_{j}(\omega)$ of the probe qubit in state $\left|0\right\rangle$ is measured from the NMR signal spectrum (Fig. 3(b)).

By scanning the frequency $\omega$ of the probe qubit through the discrete set $\omega_{k}$, we obtain the energy spectrum of the water molecule shown in Fig. 3(a). Figure 3(b) shows the NMR signal spectrum in two different cases. For the resonant frequency setting $\omega_{a}=0.22$ (the first peak in Fig. 3(a)), the probe qubit is de-excited from state $\left|1\right\rangle$ to state $\left|0\right\rangle$ with probability $0.4531$, which indicates that one of the eigenstates of the system has the eigenvalue $E_{1}=E_{0}+\omega_{a}=-83.98$. In contrast, the spectrum at the non-resonant frequency $\omega_{b}=1.80$ shows that the quantum simulator remains in the state $\left|1000\right\rangle$ with the transition probability near zero, indicating that there is no resonance at this point. In total, there are four peaks observed in the energy spectrum, and thus four energy eigenvalues of the water molecule are measured as $-83.98$, $-83.40$, $-82.66$, $-82.38$, which are in good agreement with theoretical expectations: $-83.9731$, $-83.4010$, $-82.6604$, $-82.3763$. The difference between the theoretical and experimental values, which is below $0.01$ for all four energies, is limited by $\Delta\omega$. More accurate results can be obtained by focusing on the resonance peaks obtained in the first run of the algorithm and repeating the algorithm with smaller values of $c$ and $\Delta\omega$.

(iv) Preparation of eigenstates. As discussed above, once an energy $E_{j}$ has been identified, one can run the algorithm with appropriately set parameters to deterministically prepare the state $|E_{j}\rangle$. Alternatively, one can also resort to a heralded state preparation approach, keeping only those instances in which the probe and ancillary qubits are in the state $|01\rangle_{12}$. The heralded approach can be useful for discarding some of the unsuccessful instances in the presence of decoherence and imperfections in the implementation. Here, we have implemented this protocol, as we demonstrate using a tomography procedure  Sup ; Lee (2002). Figure 4 shows that the experimentally reconstructed density matrix of the final state for different values of $\omega$. At the resonant frequency $\omega_{a}=E_{1}-E_{0}=0.22$, the system makes a transition from the reference state $\left|E_{0}\right\rangle$ in the subspace of ${\left|{10}\right\rangle}_{12}$ to the ground state $\left|E_{1}\right\rangle$ in the subspace of ${\left|{01}\right\rangle}_{12}$ with the probability $P_{1}(\omega_{a})=0.4531$. After renormalization, the density matrix in the subspace of ${\left|{01}\right\rangle}_{12}$ has a closeness cor of $99.32\%$ with the ground state of the water molecule, indicating that the system is faithfully prepared in the ground state when the transition occurs. The excited state $\left|E_{4}\right\rangle$ is prepared by setting the frequency at $\omega_{c}=E_{4}-E_{0}=1.82$, with a closeness of $99.34\%$ and a success probability of $0.8264$. For comparison, when we set the frequency outside the resonance regions, such as $\omega_{b}=1.80$ or $\omega_{d}=1.84$, the simulated system remains in its initial state $\left|10E_{0}\right\rangle$ (Fig.4(b,d)).

The deviations of the experimental final states from the expected eigenstates result mainly from two factors: (1) the fluctuations of the strength of the NMR signal that are around $1\%$, (2) the imperfection of GRAPE pulses with theoretical infidelities up to $\sim$ 0.5% which leads to a fluctuation of the spectra of around $2\%$. The high fidelities in the tomography data confirm the feasibility of the quantum algorithm for preparing eigenstates.

As we have mentioned above, the input state $\left|E_{0}\right\rangle$ of our algorithm does not need to be a good approximation to the eigenstate of interest. Moreover, since the eigenstates can be prepared deterministically and non-destructively in our algorithm, one can use the eigenstates that are prepared in the first run of the algorithm as input states for further transitions. By allowing a sequence of multiple transitions, the quantum simulator can explore a large portion of the Hilbert space, including states that have little or no overlap with the initial state. In contrast, when using the PEA for example, the algorithm only projects the initial state onto different energy eigenstates that have substantial overlap with the initial state. The transition operator $B$ is one of the inputs that can be varied in the algorithm. While we have chosen a simple operator in this experiment, the operator could be designed based on any information that we might have about the target state, or it could be varied in different runs to reduce the probability of missing a transition with a very small matrix element.

Considering the scaling of the algorithm, two factors should be considered: (1) the number of experiments that need to be performed to obtain an energy spectrum in a given scanned energy range. (2) the time needed to run the circuit in each experiment. The first factor, i.e. the number of experiments, scales as $O(1/\epsilon^{2})$ where $\epsilon$ denotes the desired uncertainty of the eigenvalue $E_{j}$. The second factor is determined by the computational cost of implementing the unitary operation $e^{-iH\tau}$, which obeys the rules of implementing the time evolution operator of a realistic system, with a standard additional overhead for the ancillary and probe qubits. Although the exact scaling laws of simulating a realistic system are not fully known and can vary depending on the details of the implementation Kivlichan et al. (2018); Heyl et al. (2018), one can safely say that the computational resources required scale polynomially with system size, as opposed to the exponential scaling expected for classical computers. The resources required also scale polynomially with the desired accuracy $\epsilon$ of the eigenvalues Sup .

To conclude, we have experimentally demonstrated a quantum algorithm for calculating the eigenenergies and eigenstates of a known Hamiltonian without having a good initial approximation for any of its eigenstates. As an example, we successfully obtained the lowest four energy eigenvalues of the water molecule by simulating resonant transitions between the system and a probe qubit, and prepared the corresponding eigenstates with high fidelity on a four-qubit NMR quantum simulator. Our quantum algorithm provides an additional option in the set of quantum algorithms for simulating eigenvalue problems and preparing eigenstates of large, complicated Hamiltonians.

I Supplementary Material

II The Hamiltonian of the Water molecule

The electronic Hamiltonian of the water molecule in the form of second quantization Levine (2000) is $H_{S}=\sum\limits_{p,q}\left\langle p\right|\left({T+{V_{N}}}\right)\left|q\right\rangle a_{p}^{+}{a_{q}}-\frac{1}{2}\sum\limits_{p,q,r,s}\left\langle p\right|\left\langle q\right|{V_{e}}\left|r\right\rangle\left|s\right\rangle a_{p}^{+}a_{q}^{+}{a_{r}}{a_{s}}$ where $|p\rangle$, $|q\rangle$, $|r\rangle$, $|s\rangle$ denote the one-particle states, ${a_{p}^{+}}$ is its fermionic creation operator, and operators $T$, $V_{N}$, and $V_{e}$ are the one-particle kinetic energy operator, nuclear attraction operator and two-particle electron repulsion operator, respectively. The Hartree-Fock wave function for the electronic structure of the water molecule is ${(1{a_{1}})^{2}}{(2{a_{1}})^{2}}{(1{b_{2}})^{2}}{(3{a_{1}})^{2}}{(1{b_{1}})^{2}}$. The number of qubits needed to represent the water molecule on a quantum simulator can be optimized in terms of the state-mapping technique presented in Ref.Wang et al. (2008). Considering the ${}^{1}A_{1}$ symmetries and using the STO-$3$G basis set and freezing the first two $a_{1}$ orbitals, the $b_{1}$ orbital and the first $b_{2}$ orbital, we construct a model space that includes the $3a_{1}$, $4a_{1}$ and $2b_{2}$ orbitals for the low-energy Hamiltonian $\mathcal{H}_{S}$ of water molecule, and perform configuration interaction (CI) calculations Szabo and Ostlund (1996). Accordingly, the multi-reference-configuration interaction (MRCI) space is composed of four configuration state functions which require two qubits to represent the state of the water molecule, and the Hamiltonian matrix is given by Ernest (1975):

[TABLE]

whose eigenvalues are $-83.9731$, $-83.4010$, $-82.6604$ and $-82.3763$. In this paper, energies and time are recorded in units of Hartree and Hartree*-1*, respectively.

In the general case, one can always write down the Hamiltonian of a chemical system in the second quantization formalism and map it to operators that act on qubits via the Jordan-Wigner transformation. The Hamiltonian contains only a polynomial number of local interaction terms, which can be transformed to operators that only act on a few qubits, and can be implemented efficiently. Thus, one can safely say that the computational resources required for simulating a chemical system on a quantum computer scale polynomially with system size, as opposed to the exponential scaling expected for classical computers.

III The tomography procedure

We implement a tomography procedure to read the final state in the subspaces of ${\left|{01}\right\rangle}_{12}$ and ${\left|{10}\right\rangle}_{12}$. Since the natural abundance of this NMR sample with just one ${}^{13}\text{C}$ carbon is about $1\%$, we read out all qubits via the ${}^{13}\rm{C}$ channel to distinguish these molecules against the large background. This partial tomography is achieved by applying eight combined readout pulses: $Y_{\frac{\pi}{2}}\text{III}$, $Y_{\frac{\pi}{2}}\text{III}\ast\text{SWAP}_{12}$, $Y_{\frac{\pi}{2}}\text{III}\ast\text{SWAP}_{13}$, $Y_{\frac{\pi}{2}}\text{III}\ast\text{SWAP}_{14}$ for measuring all diagonal elements, $\text{SWAP}_{13}$, $\text{SWAP}_{14}$ for detecting single-quantum coherences, and $\text{III}Y_{\frac{\pi}{2}}\ast\text{SWAP}_{13}$, $\text{III}X_{\frac{\pi}{2}}\ast\text{SWAP}_{13}$ for measuring zero- and double-quantum coherences, in which I is the unity operator; $X_{\frac{\pi}{2}}$ and $Y_{\frac{\pi}{2}}$ represent the $\pi/2$ rotations along the $x$ and $y$ axis; $\text{SWAP}_{ij}$ denotes the swap operation between the $i$th and $j$th qubits.

IV The Scalability of the algorithm

In our algorithm, we note that the state of the system register changes as a result of the resonant transitions. As a result, by using multiple transitions the simulator can explore a large portion of the Hilbert space, including states that have little or no overlap with any guess state. In contrast, when using the PEA for example, the algorithm is repeated $N$ times, giving the energies of states that have $N|\langle E_{j}|E_{guess}\rangle|^{2}>1$. The algorithm only projects the initial state onto different energy eigenstates that have substantial overlap with the initial state. Thus, the PEA algorithm is sensitive to the overlap between the input state and the measured eigenstates as are other similar algorithms. The algorithm implemented in our work has the advantage that even if the direct transition from the initial input state to a given eigenstate is very weak, one can use eigenstates that are prepared using the algorithm as input states for further transitions. The transformation of the initial state as a result of the resonant transitions is therefore an important advantage of our algorithm.

Considering the scaling of the algorithm, as it essentially performs quantum simulations of the dynamics, it follows the same scaling laws (as a function of system size) as quantum simulation. It therefore gives an exponential speedup over classical algorithms when dealing with large, complicated Hamiltonians. Two factors should be considered in relation to the efficiency of the algorithm: (1) the number of experiments that need to be performed to obtain an energy spectrum in a given scanned energy range. (2) the time needed to run the circuit in each experiment. The total number of experiments is equal to $KL$, while $K$ represents the number of experiments required to have a statistical estimation of the transition probability $P_{j}(\omega)$ for each frequency $\omega$, and $L$ represents the number of frequency points in the spectrum. If one wishes to determine the eigenvalue with uncertainty $\epsilon$, the total number required will scales as $O(1/\epsilon^{2})$. We note that one does not need to cover the entire frequency domain with the same intervals of frequency points in the spectrum, i.e. one can zoom in on the resonances and keep a lower density of frequency points elsewhere. Hence, depending on the problem, one might be able to obtain better scaling than $O(1/\epsilon^{2})$.

The second factor, i.e. the time needed to run the algorithm, is determined by the computational cost of implementing the unitary operation $e^{-iH\tau}$. In the experiment, the system-probe coupling strength is chosen such that $c\left|d_{j}\right|\sim\epsilon$. The evolution would then need to be performed for a time $\pi/(2c\left|d_{j}\right|)$, and hence the time will scale as $\tau\sim 1/(2\epsilon)$. The implementation of $e^{-iH\tau}$ obeys the rules of implementing the time evolution operator of the quantum system of interest, i.e. $e^{-iH_{s}\tau}$, with a standard additional overhead for the control and probe qubits. Since $H_{S}$ is a physical Hamiltonian, it typically includes only one and two-body interactions and could be represented by $O(N^{2})$ second-quantized terms, where $N$ is the number of spin orbitals Babbush et al. (2018). Although the exact scaling laws of the quantum simulation are not fully known and can vary depending on the details of the implementation Kivlichan et al. (2018); Heyl et al. (2018), one can safely say that the computational resources required for simulating a physical system scale polynomially with system size, as opposed to the exponential scaling expected for classical computers.

Furthermore, while we have used an NMR setup, the algorithm is based on the circuit model of quantum computing. It can therefore be implemented on other platforms that support the circuit model. For example, there are now superconducting devices with 10-100 qubits that can be used to implement the algorithm for larger molecules. In those devices it is the coherence time rather than the number of qubits that limits the size of the problem that can be addressed. Judging from recent experiments, implementing our algorithm with 5-10 system qubits seems to be a reasonable estimate for the system size that can be tackled on present-day quantum computers.”