Adiabatic Quantum Simulators

In his famous 1981 talk, Feynman proposed that unlike classical computers, which would presumably experience an exponential slowdown when simulating quantum phenomena, a universal quantum simulator would not. An ideal quantum simulator would be controllable, and built using existing technology. In some cases, moving away from gate-model-based implementations of quantum computing may offer a more feasible solution for particular experimental implementations. Here we consider an adiabatic quantum simulator which simulates the ground state properties of sparse Hamiltonians consisting of one- and two-local interaction terms, using sparse Hamiltonians with at most three-local interactions. Properties of such Hamiltonians can be well approximated with Hamiltonians containing only two-local terms. The register holding the simulated ground state is brought adiabatically into interaction with a probe qubit, followed by a single diabatic gate operation on the probe which then undergoes free evolution until measured. This allows one to recover e.g. the ground state energy of the Hamiltonian being simulated. Given a ground state, this scheme can be used to verify the QMA-complete problem LOCAL HAMILTONIAN, and is therefore likely more powerful than classical computing.


I. INTRODUCTION
Computer simulation of quantum mechanical systems is an indispensable tool in all physical sciences dealing with nanoscale phenomena. Except for specific and rare cases 1 , classical computers have not been able to efficiently simulate quantum systems, as in all known techniques at least one of the computational resources required to perform the simulation scales exponentially with the size of the system being simulated.
Numerous classical approximative methods, such as density functional theory (DFT 2 ) and quantum Monte Carlo (QMC 3 ) have been developed to address various aspects of the efficiency problem, but no known polynomially scaling methods are universally applicable. Each suffers from particular deficiencies such as the fermionic sign problem of QMC or the approximate exchange-correlation functionals of DFT. Quantum computers on the other hand, as conjectured by Feynman 4 , may be used to simulate other quantum mechanical systems efficiently. Feynman's conjecture was subsequently proven and expanded leading to the rapidly growing area of study known as quantum simulation 5-12 . Quantum simulation is expected to be able to produce classically unattainable results in feasible run times, using only a modest number of quantum bits. For example, calculating the ground state energy of the water molecule to the level of precision necessary for experimental predictions (≈ 1 kcal/mol) -a problem barely solvable on current supercomputers 13 -would require roughly 128 coherent quantum bits on a gate-model quantum computer 9 . However, this calculation would also require on the order of billions of quantum gates and thus quantum error correction would be essential 14 . The large number of gates and qubits needed for practical applications of quantum simulation has lead current experimental work to fine tune small systems 11,15,16 .
Theoretical quantum simulation falls into two main categories: dynamic evolution and static properties. Both categories rely heavily on the Trotter decomposition to handle non-commuting terms in the Hamiltonian when mimicking the unitary time propagator of the system to be simulated. To approximate evolution under a Hamilto- , one applies the sequence of unitary gates {e −iHit/n } k i=1 a total of n times. As number of repetitions n tends to infinity the error caused by the noncommutativity vanishes and the approximation converges to the exact result 5 . If each time slice requires a constant number of gates independent of the parameter t/n, then reducing the error by repeating the sequence n times can become prohibitively expensive for high accuracy applications.
Constructing a practical method of quantum simulation is a significant challenge (see in particular the recent work 14 which analyzes fault tolerance in a gate-model quantum simulation algorithm). Moving away from the gate model may offer a less resource-intensive, and consequently a more feasible solution. This letter addresses the problem by presenting a hybrid model of quantum simulation, consisting of an adiabatically controlled simulation register coupled to a single gate-model readout qubit. Our scheme can simulate the ground state properties of arbitrary spin graph Hamiltonians. Furthermore, it allows the measurement of the expectation value of any constant k-local observable using a k + 1-local measurement Hamiltonian. Even in the presence of considerable noise at least some information can be extracted.
In most quantum computing architectures, the natural interactions are 2-body. However, under certain conditions kbody interactions can be well approximated using techniques such as perturbative Hamiltonian gadgets [17][18][19] or average Hamiltonian methods 20 . In particular, reference 21 considered the mapping between a given n-qubit target Hamiltonian with k-body interactions onto a simulator Hamiltonian with two-body interactions.
Structure of this paper: We will continue by giving an overview of the simulator design, including algorithm initialization, adiabatic evolution and measurement. One of the key motivations for the proposed model is its resistance to certain types of noise which is outlined in Section III, and we conclude in Section IV. The supporting material in the Appendix further details the specifics of our method (including the full numeric simulation of the technique in the presence of noise).

II. SIMULATOR OVERVIEW
We consider the simulation of systems represented by finite collections of spins acted on by a time-independent Hamiltonian described by a graph G = (V, E) -e.g. the Heisenberg and Ising models. Each graph vertex v ∈ V corresponds to a spin acted on by a local Hamiltonian L v , and each edge e ∈ E to a pairwise interaction K e between the involved vertices. The Hamiltonian H T we wish to simulate is given by The simulator consists of an adiabatically controlled simulation register S with Hamiltonian H S , and a probe register P which will be acted on by gate operations and measured projectively. We will engineer the probe P such that it behaves as a controllable two-level system with orthonormal basis {|p 0 , |p 1 }. Without loss of generality, the probe Hamiltonian can be expressed as H P = δ|p 1 p 1 |, where δ is the spectral gap between the probe's ground and first excited states.
Initialization: We will first set the simulation register Hamiltonian H S to H S,I and prepare S and P in their respective ground states. The Hamiltonian H S,I has a simple ground state which can be (i) computed classically and (ii) prepared experimentally in polynomial time -such as a classical approximation to the ground state of the simulated system. The simulator Hamiltonian is thus initially given by Adiabatic evolution: According to the adiabatic theorem 22 , a quantum system prepared in an energy eigenstate will remain near the corresponding instantaneous eigenstate of the time-dependent Hamiltonian governing the evolution if there are no level crossings and the Hamiltonian varies slowly enough. By adjusting the simulator parameters, we adiabatically change H S from H S,I to H S,T , the fully interacting Hamiltonian of the system to be simulated.
Let us denote the ground state of H S,T as |s 0 . At the end of a successful adiabatic evolution P is still in its ground state |p 0 , and S is in (a good approximation to) the ground state of the simulated system, |s 0 . Hence the simulator is now in the ground state |g = |s 0 ⊗ |p 0 of its instantaneous Hamiltonian The computational complexity of preparing ground states of quantum systems has been studied [17][18][19] . It is possible to prepare a desired ground state efficiently provided that the gap between the ground and excited states is sufficiently large 22 . This depends on the adiabatic path taken. In general, finding the ground state energy of a Hamiltonian, even when restricted to certain simple models, is known to be complete for QMA, the quantum analogue of the class NP [17][18][19] . Hence, it is considered unlikely that even a quantum computer would be able to do this in polynomial time.
In fact there are physical systems such as spin glasses in nature which may never settle in their ground state. However, it is hoped (but not proven) that a host of realistic systems will satisfy the property of having a large energy gap on even simple linear adiabatic paths and thus be amenable to quantum simulation algorithms which rely on adiabatic state preparation.
Measurement: The measurement procedure begins by bringing S and P adiabatically into interaction. The simulator Hamiltonian becomes where the operator A corresponds to an observable of the simulated system that commutes with H S,T . Assuming that A can be decomposed into a sum of two-local operators, the interaction term H SP involves three-local interactions. These terms can be implemented using either Hamiltonian gadget techniques or average Hamiltonians (see the supporting material for details). Let us use a 0 := s 0 |A|s 0 to denote the expectation value we wish to measure. If a 0 + δ ≥ 0, then |g is also the ground state of H 2 (see the supporting material for proof). At time t = 0 we apply a Hadamard gate to the measurement probe which puts it into a superposition of its two lowest states. This is no longer an eigenstate of H 2 , and the system will evolve as where ω := (a 0 + δ)/ . We have thus encoded the quantity a 0 , a property of the ground state of H S,T , into the time dependence of the probe P . After a time t, we again apply a Hadamard pulse to the probe, resulting in the state and then measure the probe. The probability of finding it in the state |p 0 is If we have non-demolition measurements (see e.g. 23 ) at our disposal, then measuring the probe does not disturb the state of the simulator which can be reused for another measurement. One repeats the measurement with different values of t until sufficient statistics have been accumulated to reconstruct ω and hence a 0 -this is reminiscent of Ramsey spectroscopy 24 and hence should seem natural to experimentalists. In essence, we have performed Kitaev's phase estimation algorithm 25 , using the interaction Hamiltonian H SP instead of the controlled unitary.
If the ground state subspace of H S,T is degenerate and overlaps more than one eigenspace of A, or the simulation register S has excitations to higher energy states at the beginning of the measurement phase, the probability of finding the probe in the state |p 0 is given by a superposition of harmonic modes. For example, for a thermalized state with inverse temperature β, we obtain P 0 (t) = 1 2 1 + 1 where the first summation index runs over the energy eigenstates and the second over the eigenstates of A in which energy has value E k , and ω k,l = (a k,l + δ)/ .

III. EFFECTS OF NOISE
To assess the effects of noise on the simulation method, we performed a numerical simulation of the simplest nontrivial implementation of the hybrid simulator, consisting of two simulator qubits and one probe qubit, with a randomly chosen H T . Each qubit was coupled to its own bosonic heat bath with an Ohmic spectral density using the Born-Markov approximation 26 . The qubit-bath couplings were chosen such that the resulting single-qubit decoherence times T 1 and T 2 are compatible with recent superconducting flux qubit experiments with fully tunable couplings (see for example 27 ). The noise model is described further in the supporting material. The observable A was chosen to be H S,T , the simulated Hamiltonian itself, which means that the ground state energy was being measured.
We simulated measurements on 40 evenly distributed values of the time delay t. At each t i we performed 50 measurements, and averaged the results to get an estimate of P 0 (t i ). An exponentially decaying scaled cosine function was then fitted to this set of data points to obtain an estimate for ω and thus for s 0 . The fit was done using the MATLAB LSQNONLIN algorithm, guided only by fixed order-of-magnitude initial guesses for the parameter values.
The results of the simulation are presented in Fig. 1. The noise, together with the slight nonadiabaticity of the evolution, cause excitations out of the ground state which result in a signal consisting of multiple harmonic modes. However, the ground state mode still dominates and with a realistic level of noise and relatively modest statistics we are able to reconstruct ω to a relative precision of better than 0.01. This includes both the uncertainty due to finite statistics, and the errors introduced by the environment through decoherence and the Lamb-Stark shift. In an experimental implementation there may also be other noise sources not considered here related to the measurement process itself, as well as systematic (hardware) errors in the quantum processor (qubits, couplers etc.).

IV. CONCLUSION
The presented simulation scheme significantly differs from gate-model methods in that it does not require millions of coherent gate operations. Instead, an adiabatic control sequence is used, which requires much less complicated control hardware and is likely to be more robust against errors. A simple experimental implementation could be feasible with present-day technology. A direct comparison to a gate-model simulator using phase estimation to extract the result is difficult because such a design would have to use quantum error correction to be useful at the noise levels we have considered. At the measurement stage we do require single-qubit gate operations and measurements, but only on the probe qubit. Furthermore, these operations are relatively simple to implement compared to a full Trotter decomposition of the simulated Hamiltonian. It is well known that fault tolerance in the adiabatic model depends on the energy difference between the Hamiltonian's ground and first excited states. Several methods for the adiabatic model have been proposed which could in principle be used to increase the fault tolerance of our protocol [28][29][30] .
In order to simulate a system of n qubits with a 2-local Hamiltonian described by the graph G, ideally our scheme requires one probe qubit for the readout and n simulation qubits. Additionally, if Hamiltonian gadgets are used to implement 3-local interactions, one ancilla qubit is required for each two-local term in the simulated observable (represented by an edge in the corresponding graph). In practice the number of ancillas may be slightly higher if more involved types of gadgets are used to implement the 3-local interactions. The total number of qubits required thus scales as O(n) for sparse graphs and O(n 2 ) for maximally connected ones.
We will apply a Hamiltonian H p to the mediator qubit m as well as a Hamiltonian V to qubits 1-3 and m. The Hamiltonian H p + V has a ground state that is ǫ-close to the desired operator JA 1 ⊗ A 2 ⊗ A 3 ⊗ |0 0|.
where y is some Hamiltonian already acting on qubits 1 − 3 as well as a possible larger Hilbert space. Note the gadget above assumes A 2 i = 1, ∀i = 1, 2, 3. The so called self-energy expansion (A3) under appropriate conditions is known to provide a series approximation to the low-lying eigenspace of an operator. To verify that the Hamiltonian H p + V gives the desired approximation, one relies on expansion of the self energy to 4 th order (here the higher order terms give rise to effective interactions greater than second order): where the operator is written in the {|0 , |1 } basis of m. One considers the range |z| ≤ 2|J| + ǫ and notes that and applies the Gadget Theorem 17 . Before concluding this section, we note that care must be taken when adiabatically evolving gadgets. Reference 31 contains a proof that the linear path Hamiltonian is universal for adiabatic quantum computation. Universality (and hence a non-exponentially contracting gap) remains when the locality of the construction is reduced using perturbative gadgets 18 .

Average Hamiltonian Method
An alternate method of generating the special Hamiltonians we require is to make use of the time average of a series of simpler generating Hamiltonians 20 . It has long been known that by regularly switching between a set of fixed Hamiltonians {H i }, it is possible to approximate time evolution under any other Hamiltonian H, provided that H is contained within the algebra generated by {H i }. This fact lies at the heart of both average Hamiltonian theory and geometric control theory. Over the years many methods have been developed for making the approximations accurate to high order, however here we will focus only on approximations correct to first order.
In order to construct an average Hamiltonian we will make use of a first order approximation to the Baker-Campbell-Hausdorff formula: From this we obtain formulae for approximating the exponential of both the sum and the Lie bracket of A and B to first order.
These equations can be related to time evolution under some Hamiltonians, H i , by replacing A and B with operators of the form iH i t/ . Clearly by applying these rules recursively, it is possible to generate any Hamiltonian in the Lie algebra generated by the initial set of static Hamiltonians. We note that the combination of any pairwise entangling Hamiltonian, such as an Ising Hamiltonian together with tunable local Z and X fields is sufficient to generate the full Lie algebra su(2 N ) for an N qubit system, and so can be used to approximate an arbitrary Hamiltonian.
Although the time-varying Hamiltonian means that the system does not have a ground state in the normal sense, the average energy of the system is minimized when the system is within O(t) of the ground state of the target average Hamiltonian. As a result of this, if the time scale for switching between Hamiltonians is small compared to the time scale for the adiabatic evolution of the system, the system will behave as if it was experiencing the average Hamiltonian.

Noise model
The noise model used in our MATLAB simulation consists of a separate bosonic heat bath coupled to each of the qubits. The baths have Ohmic spectral densities, where the cutoff ω c was chosen to be above every transition frequency in the system, and the baths are assumed to be uncorrelated. Each bath is coupled to its qubit through an interaction operator of the form Using the Born-Markov approximation we obtain an evolution equation which is of the Lindblad form 26 . Denoting the level splitting of a qubit by ∆, we obtain the following uncoupled single-qubit decoherence times: where β = 1 kB T . Given T 1 , T 2 , T and ∆, we can solve the bath coupling parameters λ and α separately for each qubit, and then use the same noise model in the fully interacting case. The values used for the parameters are presented in Table I. The single-qubit decoherence times T 1 and T 2 were chosen to be compatible with recent coupled-qubit experiments with fully tunable couplings such as 27 . The bath temperature T and the energy scale ω 0 for of the effective Hamiltonian to be simulated were likewise chosen to match those found in typical superconducting qubit experiments. To simulate manufacturing uncertainties, the actual values of the T 1 and T φ parameters for the individual qubits were drawn from a Gaussian distribution with a small standard deviation.

Measurement
The pre-measurement system Hamiltonian is The operators H S and H P can be expanded in terms of their eigenvalues and eigenstates using the (possibly degenerate) spectral decomposition and correspondingly for H P . Let the states of systems S and P begin in their ground state subspaces, the full normalized state |g of the simulator belonging to the ground state subspace of the non-interacting Hamiltonian H 1 : |g ∈ span{|s 0,k ⊗ |p 0,l } k,l . (A13) S and P are brought adiabatically into interaction with each other. The Hamiltonian becomes where the operator A corresponds to an observable of the simulated system that commutes with H S , and is the projector to the ground state subspace of H P . Because A and H S commute, they have shared eigenstates: H S |s k,j = s k |s k,j , (A16) A|s k,j = a k,j |s k,j . (A17)

a. A ground state lemma
We will now show that Hamiltonians H 1 and H 2 have the same ground state subspace given that a min + p 1 − p 0 > 0 where a min is A's lowest eigenvalue. Hence, if A's lowest eigenvalue a min ≥ 0 then |g is the ground state of H 2 = H 1 + H SP as well. If A is not nonnegative, we can perform the transformation H ′ S = H S + a min I, H ′ P = H P − a min Π p0 , A ′ = A − a min I.
This leaves H 2 invariant, but makes A ′ nonnegative. As long as a min + p 1 − p 0 > 0 then {|p 0,m } m still span the ground state subspace of H ′ P and the above analysis remains valid.

b. Ground state degeneracy and excitations
If the ground state subspace of H S,T is degenerate and overlaps more than one eigenspace of A, or the simulation register S has excitations to higher energy states at the beginning of the measurement phase, we need a more involved analysis of the measurement procedure. Assume the pre-measurement-phase state of the simulator is given by = (e qHS ⊗ e qHP )(e qA ⊗ (1 1 P − Π p0 ) + 1 1 S ⊗ Π p0 ).