Preparation of many-body states for quantum simulation
Nicholas J. Ward, Ivan Kassal, and Alán Aspuru-Guzik∗
Department of Chemistry and Chemical Biology, Harvard University, Cambridge, MA 02138 While quantum computers are capable of simulating many quantum systems efficiently, the simulation algorithms must begin with the preparation of an appropriate initial state. We present a method for generating physically relevant quantum states on a lattice in real space. In particular, the present algorithm is able to prepare general pure and mixed many-particle states of any number of particles. It relies on a procedure for converting from a second-quantized state to its first-quantized counterpart. The algorithm is efficient in that it operates in time that is polynomial in all the essential descriptors of the system, such the number of particles, the resolution of the lattice, and the inverse of the maximum final error. This scaling holds under the assumption that the wavefunction to be prepared is bounded or its indefinite integral known and that the Fock operator of the system is efficiently simulatable.

arXiv:0812.2681v2 [quant-ph] 29 May 2009

I.

INTRODUCTION

Simulating quantum systems on a conventional computer requires resources that generally scale exponentially with the size of the system. Feynman proposed to solve this problem using a quantum machine that would be able to mimic the properties of the quantum system [1]. Subsequently, it has been demonstrated that quantum computers would be able to simulate the time-dependent Schrödinger equation for many systems of interest using resources that scale polynomially with the size of the system [2, 3, 4, 5, 6, 7, 8, 9]. However, all such simulations require the preparation of an appropriate initial state, which must be preparable to within a chosen error. In this work, we focus on the preparation of states on a gatemodel quantum computer. Our techniques can therefore complement those developed for the preparation of states in other models of quantum computation, such as adiabatic quantum computing [6, 10]. In general, we will call a state on n qubits “efficiently preparable” if it can be prepared, to within error ε, using poly(n, ε−1 ) elementary (one- and two-qubit) quantum gates. Unfortunately, the efficiently preparable states form only a small subset of all quantum states. This is because a general state on n qubits contains 2n amplitudes, and therefore one needs O(2n ) gates to prepare it [11]. Indeed, statepreparation algorithms are known that almost reach this lower bound [12, 13]. In this work, we show that if wavefunctions are represented on a grid in real space, then most quantum states of physical interest are efficiently preparable. This is of interest because efficient, grid-based simulation algorithms are known for physically realistic systems [3, 4, 5, 8]. Our work extends that of Zalka, who, in introducing realspace quantum simulation [4], also provided the first state preparation algorithm. However, his procedure is able to prepare only states of single particles or uncorrelated manyparticle systems. We show how to use Zalka’s single-particle wavefunctions as building blocks, permitting the preparation of general superpositions and mixed states of an arbitrary

number of particles. Our approach is motivated by electronicstructure theory, in that we we choose particularly convenient single-particle bases in which to expand more complicated states. We use the single-particle eigenstates to form Slater determinants (configurations), superpositions of which are used to express general many-particle states. Our scheme is essentially a method for translating states in second quantization to the corresponding states in first quantization. This has two advantages. First, many useful states that might be needed in first quantization are easily prepared in second quantization [14]. In particular, we can prepare eigenstates of operators if our scheme is combined with full configuration interaction (FCI) [15], an exact diagonalization method. FCI is classically an exponentially hard problem due to the exponential growth of the number of configurations with system size, but it can be computed on a quantum computer in polynomial time [6]. The quantum FCI operates in second quantization, and can compute, for example, the ground state wavefunction of a molecular system. The second benefit of our method is that it is often easier to simulate time-evolution in real space than in Fock space. For instance, every simulation in second quantization would require a separate set of basis-set–dependent operators and there might be some processes, such as ionization, which could not be adequately described using a small, localized basis set. In firstquantization, however, all problems of chemical interest can be efficiently simulated by direct simulation of the molecular Hamiltonian in real space [8]. This paper is organized as follows. We first consider the preparation of many-particle states in which all the particles are the same. There are three steps: the preparation of singleparticle eigenstates in a chosen basis, the preparation of manyparticle configurations, and finally the preparation of superpositions of configurations. We discuss the preparation of mixed states, after which we turn to systems with many different types of particles. Again, we consider the preparation of configurations, their superpositions, and mixed states. We close by showing that the algorithm is efficient in that its run-time is polynomial in the size of the system, the number of qubits used to encode the wavefunction, and the inverse of the maximum allowed error.

∗ Electronic

address: aspuru@chemistry.harvard.edu

2
II. ONE TYPE OF PARTICLE

ing a grid of 2l points, and a basis state φ(x) normalized over a length L, the algorithm first performs the transformation
2l −1

Our algorithm translates from second to first quantization, and therefore depends on the basis which is chosen for the representation of the second-quantized states. We require a finite orthonormal basis of functions {φ1 , . . . , φM }, which are ˆ the eigenstates of a known operator F on an M -dimensional, single-particle Hilbert space. In our analogy with electronicˆ structure theory, F would be the Fock operator for a single particle [15], and indeed we expect that the algorithm would ˆ be at its most useful if F is chosen as the Fock operator of ˆ an actual system. Although the form of F can be arbitrary, subject to a few restrictions below, we will take advantage of ˆ the analogy and refer to F as the Fock operator. We will also ˆ say, for example, that two eigenstates of F are degenerate if their energies (the corresponding eigenvalues) are the same. To ensure that the overall state-preparation algorithm scales ˆ efficiently, we require that F can be efficiently simulated on a quantum computer, i.e., that the simulation time scales polynomially with the size of the system. More precisely, if there are m particles occupying the M orbitals and the simulation is done on a grid of 2l sites (see below), then, for ˆ any t and any ε > 0, there exist a unitary U , composed of poly(m, M, l, t, ε−1 ) elementary (one- and two-qubit) quanˆ ˆ tum gates, such that U − e−iF t ≤ ε. Intuitively, this means

|0 → |φ =

φ x
x=0

L 2l

|x ,

where each integer-valued state |x is a position on the suitably scaled grid. This state is generated from the state |000 . . . by redistributing its amplitude l times across the eigenstates |x . To perform the redistribution correctly, we calculate the integrals
L (k+1) 2i

Ii,k =

k

L 2i

|φ(x)|2 dx |φ(x)|2 dx

L (k+2) 2i

,

(1)

kL 2i

for k = 0, . . . , 2i − 2 and i = 1, . . . , l. The fraction Ii,k is simply the probability that a particle in the (k + 1)th subdivision of size L/2i is also in its left half. If the denominator in Ii,k is zero, there is no amplitude to redistribute, so we can skip this step. The first split is realized by performing a rotation on the first qubit by arccos( I1,0 ), corresponding to the transformation |0, . . . → I1,0 |0, . . . + 1 − I1,0 |1, . . . .

that given an initial state, the final state generated by the acˆ tion of F for time t can be calculated with reasonable effort and reasonable error. Several classes of Hamiltonians are known to be efficiently simulatable, and together they ensure that most physically relevant Fock operators will also be efficiently simulatable. Very generally, an operator can be efficiently simulated if its matrix in a given basis is sparse [10, 16, 17]. In particular, this includes Hamiltonians that are sums of local operators, each of which acts on only a few degrees of freedom [2, 11]. In addition, many physically realistic real-space Hamiltonians (such as those for chemical systems) can be efficiently simulated [3, 4, 5, 8]. We finally note that the requirement that the basis be orthonormal may exclude certain commonly used basis sets. Many of the usually encountered bases are appropriate, such as plane waves or molecular orbitals, which diagonalize the molecular Hartree-Fock Hamiltonian. However, nonorthogonal bases, such as Gaussian wavepackets or atomic orbitals on more than one atomic center, are not suitable for state preparation using our procedure.

This splits the norm of the initial state so that the appropriate proportion is present on each half of the grid. The subsequent finer splits are carried out in superposition using controlled rotations on each qubit. For example, after the second iteration, the correct proportion of the norm is present in each quarter of the grid. After l iterations, one obtains the desired state. Note that adding a single qubit and the corresponding rotation doubles the precision of the grid. Consequently, the absolute value of the wavefunction can be efficiently approximated to any desired accuracy. Phases can be added where necessary through phase-kickback [20]. Given a procedure that can transform |x → ei arg φ(x) |x , we can complete the preparation of |φ as
2l −1 x=0

L φ x l 2

2l −1

|x → =

ei arg φ(xL/2 ) φ x
l

x=0 2l−1

L 2l

|x

φ x
x=0

L 2l

|x = |φ .

A.

Single-particle eigenstates

A single-particle basis function φ can be prepared on a grid by the state preparation method first proposed by Zalka [4] and rediscovered independently by both Grover and Rudolph [18] and Kaye and Mosca [19]. The algorithm first prepares the absolute value of the function, followed by the addition of the phases. Specifically, given a register of l qubits, represent-

The same algorithm can be straightforwardly generalized to a three-dimensional grid, where the position eigenstates are in Cartesian coordinates and the corresponding threedimensional integrals are used. In addition, particle spin can be represented using additional qubits. A particle with spin S requires ⌈log2 (2S + 1)⌉ qubits to store its z-projection mS . In particular, there is a natural mapping between the spin of 1 spin- 2 particles and the states of a single qubit. If the Fock operator is spin-free, the eigenstates will have separable spatial and spin degrees of freedom, making the complete singleparticle state |φ |mS . Preparing the spin part of this wavefunction is relatively easy, for it suffices to initialize the spin

3 register in an integer state. If, however, the eigenstate has correlation between the spatial and spin degrees of freedom, it can be prepared using the techniques in Secs. II D and III. That is, we treat the particle as if it were a composite system— composed of a spinless, spatial part and a spin—and prepare its eigenstate using the techniques below. In what follows, we will assume that our particles are fermions and we will note where the algorithm would need to be modified for bosons.
B. Computational complexity of integration

The preceding method for preparing single-particle states requires the evaluation of integrals (1). Since this must be performed in superposition, the integrals must be computed on the quantum computer: precomputing them classically would require an exponentially large look-up table. Consequently, the computational complexity of the state preparation procedure will depend on the the cost of computing the integrals [4, 18]. An integration procedure will, given a function φ : V ⊂ ˜ Rd → R (where V is a bounded region), supply an estimate I ˜ of the integral I = V φ(x)dx such that |I − I| ≤ εI with a certain fixed probability δ (we’ll call this condition the (εI , δ) absolute error). Integrals can be evaluated either analytically or numerically. If the indefinite integrals of the basis functions are known, the definite integrals over any box on the Cartesian grid can be computed. The values of the indefinite integrals themselves can usually be computed efficiently (i.e., with polynomial cost in the desired accuracy) because they usually contain simple mathematical functions. In particular, the time it takes to retrieve n digits of any elementary function is a polynomial in n [21], and likewise for compositions of elementary functions. If the indefinite integrals are either unknown or impractical to compute, numerical techniques can be used. In particular, any classical numerical technique can, in principle, be imple˜ mented on a quantum computer. For example, computing I by Monte Carlo requires, in the worst case [22], Φ−1 (1 − δ/2)
2

Furthermore, it is known that quantum computers are able to offer a quadratic speed-up over conventional probabilistic methods of integral evaluation. Quantum integration techniques [24, 25] rely on amplitude amplification [26] to achieve a computational complexity of O(ε−1 ). This has been proven I optimal by Nayak and Wu [27, 28]. These techniques have the same general applicability as classical Monte Carlo, and will likewise succeed for any bounded function. Furthermore, the state preparation scheme of Soklakov and Schack [29], which relies on amplitude amplification, also succeeds in O(ε−1 ) I time. The preceding assumes that the function that we seek to prepare does not depend substantially on the grid spacing. We would expect that of realistic wavefunctions, assuming that the grid spacing is smaller than the smallest wavelength of the system. A useful exception are Kronecker δ-functions, defined on a grid of 2l points as φ(x) = 2l/2 δx,x0 , where x0 is a constant vector. The variance of δ-functions grows exponentially in l, and therefore they cannot be integrated efficiently by Monte Carlo or prepared efficiently using the method of Soklakov and Schack [29]. However, they can still be prepared efficiently using our techniques because their indefinite integral, the Heaviside function, can be easily computed in time independent of l. It remains to be shown that an error in the evaluated integral translates to a comparable error in the prepared function. If the ˜ integrals (1) have a maximum error εI , that is, Ii,k − Ii,k ≤ ˜ εI , and the error in the final prepared state φ is εφ = 1 − ˜ φ φ , we find that εφ ≤ lεI /2. In the case l = 1, |φ = √ √ ˜ ˜ ˜ I |0 + 1 − I |1 and φ = I |0 + 1 − I |1 . Then, ˜ ˜ assuming 0 ≤ I ≤ 1, which is necessary for φ to be an acceptable state, εφ = 1 − ˜ ˜ I I − (1 − I)(1 − I) √ ≤ 1 − 1 − εI εI , ≈ 2

σ 2 /ε2 I

samples of φ for an (εI , δ) absolute error, where σ 2 is an estimate of the variance of φ over V and Φ(z) = z (2π)−1/2 −∞ exp(−u2 /2)du is the standard normal cumulative distribution function. In particular, if φ is bounded so that φL ≤ φ(x) ≤ φU for all x ∈ V , the number of required samples is limited [22] to Φ−1 (1 − δ/2)(φU − φL )/2εI
2

(1 − εI )l/2 ≥ 1 − lεI /2 (the last inequality holds for all εI if l ≥ 2), whence εφ ≤ lεI /2. That is to say, the error in the prepared state grows only polynomially with the error in the evaluated integral, a fact that we will use later on to establish the computational cost of the state preparation algorithm.

where we have assumed that εI ≪ 1. For larger l, a sim˜ ilar analysis applies qubit-wise: one finds that φ φ ≥

.
C. Many-particle eigenstates

That is, Monte Carlo integration of any finite-variance function requires time that scales as O(ε−2 ). Acceptable waveI functions need not be continuous or even finite [23] (and hence may have infinite variance), but such examples are rather contrived and rarely encountered in practice (but see below for δ-functions).

The next step is to use Zalka’s algorithm to prepare multiparticle configurations. That is, we wish to prepare the position-space representation of a second-quantization state |n1 n2 . . . nM 2nd (a Fock eigenstate), where ni is the occupation number of the basis orbital φi (1 ≤ i ≤ M ). The

4 position-space representation of |n1 n2 . . . nM 2nd will be a Slater determinant of the occupied orbitals, and it will be an eigenstate of the many-body Hartree-Fock Hamiltonian
m

Next we will transform this state into the superposition 1 √ |σ(1, . . . , m) m! σ∈Sm
B

ˆ H=
i=1

ˆ Fi ,

(2) as follows. First let B ′ [1] = B[1]. Then assign to B ′ [i] the B[i]th natural number not present in the set {B ′ [1], B ′ [2], . . . , B ′ [i − 1]}. This leaves the quantum computer in the state 1 √ |φj1 φj2 . . . φjm m!
A

ˆ where m is the total number of particles and Fi is the Fock ˆ operator F acting on the particle i [15]. We assume that the state |n1 n2 . . . nM 2nd has already been prepared by some previous algorithm. The M basis orbitals are occupied by m particles and we let j1 , . . . , jm be the indices of the occupied orbitals. We therefore wish to perform the transformation |n1 n2 . . . nM |n1 n2 . . . nM ⊗ |0 . . . 1st → 1 sgn(σ)|σ(φj1 φj2 . . . φjm ) 2nd ⊗ √ m! σ∈Sm
2nd

⊗

σ∈Sm

|σ(1, . . . , m)

B.

(4)

Register B now contains a symmetrized state and this symmetry can be transferred to register A by sorting B while performing the same swaps on the wavefunctions in A. This yields the symmetrized state 1 √ |σ(φj1 φj2 . . . φjm ) m! σ∈Sm
A

1st

= |n1 n2 . . . nM

2nd

⊗ |Φ

1st

,

(3)

⊗ |1, 2, . . . , m

B,

(5)

which takes the input state and prepares the appropriate firstquantized Slater determinant |Φ 1st , a superposition of all the permutations on the m occupied orbitals (Sm is the symmetric group on m elements and sgn denotes the signature). Here |0 . . . 1st contains m registers for the m first-quantized occupied orbitals |φj1 1st , . . . , |φjm 1st . Note that (3) is not in general a reversible operation, as multiple input states would be mapped to the same antisymmetrized result. To ensure the algorithm is reversible, we additionally require [30] that j1 < j2 < · · · < jm . The procedure can be slightly modified if bosons are in question: then sgn(σ) is to be omitted, and the ji must satisfy j1 ≤ j2 ≤ . . . ≤ jm . The transformation (3) is accomplished in two steps. First, the occupied single-particle basis orbitals are each prepared in a separate register, forming a Hartree product: |n1 n2 . . . nM
2nd

which is what we would keep if we were interested in preparing bosonic states. To instead obtain an antisymmetrized state, we need only count the number of exchanges made in the sort, and reverse the sign of the wavefunction if it is odd. If we now eliminate the register B, A contains the desired multi-particle state |Φ 1st . The original algorithm, introduced by Abrams and Lloyd, included an additional auxilliary register C, which would then be used as an intermediate for the sorting of A and B. We eliminate this step by sorting A and B together directly.

D. Superpositions

⊗ |0 . . . 1st → |n1 n2 . . . nM 2nd ⊗ |φj1 φj2 . . . φjm

1st .

The procedure can be modified in the case of bosons by counting the occupation of each orbital and preparing that many copies in separate registers. In the next step, the Hartree product is antisymmetrized, which produces the desired Slater determinant. To complete this step, we introduce an improved form of the antisymmetrization algorithm developed by Abrams and Lloyd [30]. The algorithm begins with the m wavefunctions |φji to be antisymmetrized in register A, and m ⌈log2 m⌉ qubits in register B (where each grouping of ⌈log2 m⌉ qubits constitutes a “quword”) initialized to |0 . Using a series of controlled rotations, B is converted to the state 1 √ m!
m i=1 m−1

We now generalize the algorithm to the preparation of superpositions of many-particle states. Given a superposition of second-quantization states |ni 2nd = |n1i n2i . . . nMi 2nd , with amplitudes αi , we wish to perform the transformation αi |ni ⊗|0 . . . → |0 . . .
2nd ⊗

2nd

1st

i

i

αi |Φi

1st

.

|i

B[1]

⊗

i=1

|i

B[2]

⊗ · · · ⊗ |1

B[m] ,

which is a superposition of m! unique states consisting of m quwords each, and B[i] denotes the ith quword in register B.

The superposition on the left might come from a variety of sources. For example, an easily-prepared equal superposition of Fock states would result in an equal superposition of realspace wavefunctions. Wang et al. provide an algorithm for preparing general superpositions of Fock states on a quantum computer [14]. Alternatively, a quantum electronic-structure algorithm could be used to efficiently produce a physically relevant superposition. For example, an FCI algorithm could specify the ground state of a chemical or other many-body system in terms of a superposition of Fock states [6]. As before, we begin by applying Zalka’s state preparation algorithm to the input state. Because this linear operation is carried out in superposition, it accomplishes the transforma-

5 tion αi |ni ⊗ |0 . . .
i

2nd

1st

i

→ ⊗ |φj1 i φj2 i . . . φjm i
1st .

αi |ni

2nd

values for each state in the superposition. These eigenvalues can then be used (for example in conjunction with a lookup table) to uniquely identify the wavefunction and subtract 1 from the corresponding occupation number vector of the second-quantization state. Because this is done in superposition for every single-particle wavefunction, the input state is converted to |0 2nd . This accomplishes the total transformation αi |ni ⊗ |0 . . . ⊗ →
1st

|k |ψ . An efficient quantum Fourier transform on the control qubits will now yield the first q digˆ its of the binary expansion of θ. If we choose U such that ˆ |φi 1st = e2πiEi |φi 1st , we can use phase estimation with U enough control qubits to obtain an approximation of the enerˆ ˆ gies Ei . In particular, the natural choice U = e−iHt , where ˆ H is the Hartree-Fock Hamiltonian (2), supplies the appropriˆ ate unitary for a suitable choice of the time t. Note that U ˆ is a sum of Fock opcan be simulated efficiently because H erators which are efficiently simulatable by assumption. The energy eigenvalues are stored in an additional register containing enough qubits to provide precision that distinguishes between nearby energies. ˆ In the case that the spectrum of H is degenerate, properties other than the energy of the states need to be used to distinguish them. If the degeneracy is caused by a symmetry of the Hamiltonian, the elements of the symmetry group can be used for this discrimination, as we outline in Sec. II E. If the degeneracies are accidental, other techniques are required, and we give some suggestions in Sec. II F. In addition, the techniques in Sec. II F can be used for distinguishing eigenvalues that are exponentially close together and therefore cannot be distinguished efficiently by phase estimation. ˆ Phase estimation using both U to find energy eigenvalues and appropriate symmetry operations to distinguish degenerate states will provide us with a unique combination of eigen-

Note that the single-particle wavefunctions are now entangled with the input state. For a multi-particle eigenstate, the situation was different because the resulting state was separable. Hence, to separate the first-quantized wavefunctions from the second-quantized ones, we must “uncompute” the secondquantized states. This must be accomplished using only manipulations on the register containing the first-quantized wavefunctions |φi 1st : if we can regenerate the input state from the wavefunctions, the input register can be reset to |0 2nd as desired. Given the one-to-one correspondence between a second quantization state and the corresponding first quantization wavefunctions, regenerating the input state amounts to the problem of identifying the wavefunctions |φi 1st given only the information contained in their first-quantized representation. For non-degenerate eigenstates, each state |φi 1st can be uniquely identified using its energy, which can be obtained through the phase estimation procedure [20, 31, 32]. In genˆ eral, given a unitary U and its eigenstate |ψ , the phase estimation algorithm finds the eigenvalue of |ψ . Specifically, ˆ ˆ since U |ψ = e2πiθ |ψ , we have U k |ψ = e2πikθ |ψ . By conˆ trolled applications of the powers of U to |ψ , controlled on 2q −1 2q −1 1 1 ˆ the state 2q/2 k=0 |k , one gets 2q/2 k=0 |k U k |ψ =
1 2q/2 2q −1 2πikθ k=0 e

2nd

1st

i

|0 . . .

2nd

i

αi |φj1 i φj2 i . . . φjm i

, (6)

which is a separable state. The antisymmetrization step can now proceed in superposition as usual, resulting in the final state |Ψ 1st = i αi |Φi 1st , as desired. This completes the state-preparation algorithm for a given superposition of multiparticle states.
E. Resolving degeneracies caused by symmetry

The procedure in Sec. II D assumes that it is possible to distinguish eigenstates based on their energy. If there are degenerate states, additional operations are required to distinguish them. Degeneracies in quantum states usually arise as a result of symmetry—degeneracies that do not are called “accidental” and we treat them separately in Sec. II F. For symmetrycaused degeneracy, distinguishing degenerate states requires an understanding of how they transform under the symmetry operations of the system. All of the wavefunctions |φi 1st are eigenstates of each symmetry operation within the point group, but degenerate wavefunctions will always have different eigenvalues for at least one of the operations. Phase estimation can still be used to obtain a unique set of eigenvalues, but in addition to finding the energies, we can distinguish the wavefunctions by symmetry. By applying phase estimation using an appropriate symmetry operation as the unitary operator, we obtain additional eigenvalues to distinguish degenerate states. Because there are only a limited number of symmetries that are possible in physical systems, it will rarely be necessary to use more than a few readout qubits to retrieve all the distinguishing eigenvalues. With the exception of systems with spherical, cubic, or icosahedral symmetry, which we treat below, all systems in three-dimensional space have a symmetry point group all of whose irreducible representations are one- or two-dimensional [33]. Wavefunctions transforming as the one-dimensional irreducible representations are non-degenerate, while the ones transforming as the twodimensional irreducible representations come in degenerate pairs. Distinguishing them, therefore, requires the determination of only one symmetry eigenvalue which is different for the two wavefunctions. This is most easily done in the case of point groups Cnv , C∞v , Dn , Dnh , D∞h , and Dnd , all of which contain a C2

6 axis or a reflection plane that has character zero in all of the two-dimensional irreducible representations. In this case, one of the two degenerate wavefunctions is invariant under the reflection or C2 rotation, while the other acquires a phase of −1. To distinguish them, one would use the reflection or the C2 rotation as the unitary of phase estimation with one readout qubit (note that these operations are easy to implement, being simple linear transformations). The readout qubit, ini√ tialized in the state (|0 + |1 ) / 2, would, under the action √ of the symmetry operation, be converted to (|0 ± |1 ) / 2, depending on the acquired phase. A Hadamard gate would then return |0 or |1 , perfectly discriminating between the two eigenfunctions. Symmetry groups Cn , Cnh , and S2n have, strictly speaking, only one-dimensional irreducible representations. However, there are pairs of representations that are complex conjugates of each other, meaning that the corresponding energy levels are degenerate due to time-reversal symmetry. These pairs of conjugate representations are called “separably degenerate [34]," and the corresponding wavefunctions can be distinguished using the principal symmetry axis Cn (or S2n in the S2n groups). In each case, under the action of Cn , one of the wavefunctions acquires a phase ω and the other ω ∗ , where ω = e2πi/n (there are also cases where the pairs ac∗ quire phases such as −ω and −ω ∗ , ω 2 and ω 2 , and so on, but these do not change the procedure outlined here). Phase estimation can, as usual, measure this phase up to a certain precision. However, since 1/n usually does not have a finite binary expansion, there will be an associated error in the phase estimation. This can be reduced below an arbitrary threshold by the addition of more readout qubits, as discussed in Sec. IV. This is especially true since real physical systems almost never have Cn axes with n > 8, meaning that only several qubits will be required for readout. The cubic and icosahedral groups, T , Th , Td , O, Oh , I, and Ih , all have three-dimensional irreducible representations (and I and Ih also have four- and five-dimensional ones). Fortunately, there are plenty of reflection planes and C2 axes which can be used for discrimination just as was done in the simpler groups above. Distinguishing three or four degenerate states requires two symmetry eigenvalue comparisons (and three in the case of five-fold degeneracy). Consequently, two readout qubits are required in these cases, one for each comparison (or three qubits in the five-fold degenerate case). Degenerate states of spherically symmetric systems, such as atoms, can be distinguished by energy and by their angular momentum quantum numbers ℓ and mℓ . The maximally symmetric case is the 1/r potential, where the conservation of the Laplace-Runge-Lenz vector implies that all states with equal principal quantum number n are degenerate. If our basis contains states with n ≤ nmax , we would require O(log2 nmax ) qubits for the discrimination of the angular momentum states (that is, O(log2 nmax ) qubits each for ℓ and mℓ ). While circumstances where one encounters states of extremely high angular momentum are rare, we can see that the discrimination can be performed efficiently. The phase estimation in this case would use discrete rotations as its unitary operator. A similar approach was suggested by Zalka for the related problem of implementing unitary representations of SU(2) [35].
F. Resolving accidental degeneracies and exponentially close eigenstates

In Sec. II E, we outlined a procedure for distinguishing states that are degenerate because of symmetry. However, the eigenstates might also be accidentally degenerate or exponentially close in energy so that they cannot be efficiently distinguished by phase estimation. In those cases, it is not possible to distinguish between the (near-)degenerate states using the symmetry-based procedure. One way around these problems is to transform to another basis where the (near-)degeneracy does not arise. A way of accomplishing this is to use a perturbed Fock operator ˆ ˆ ˆ ˆ F ′ = F + V , where V is a small, efficiently simulatable perturbation that breaks the (near-)degeneracies. In a finite ˆ basis, V must also be small to ensure that the new basis can adequately describe the target state. The new eigenfunctions are obtained from the old using perturbation theory, as are the new coefficients of the state that we wish to prepare. This change of basis can be done efficiently on a classical computer, before proceeding as normal with the state preparation algorithm. For the purposes of phase estimation, the new Fock operator can be efficiently simulated by operator splitting beˆ ˆ cause both F and V are efficiently simulatable. A drawback of this procedure is that the perturbation may destroy certain desirable symmetries of the system. In some ˆ cases, this can be avoided if we choose V ∝ |φi φi |, where |φi is one of the (near-)degenerate eigenstates. In that case, ˆ ˆ F ′ and F would have the same eigenstates and no change of ˆ basis would be needed. Of course, it is possible that V in this form is not efficiently simulatable, in which case this scheme would not be efficient.
G. Mixed states

The previous sections outline the procedure for preparing general pure states, which in the chosen basis read |Ψ
1st

=
i

αi |Φi

1st .

(7)

From now on, we drop the subscript 1st for clarity. We now wish to prepare a mixed state with density operator ρ= ˆ
i

pi |Ψi Ψi | ,

where |Ψi are arbitrary pure states of the form (7) and the probabilities pi add up to 1. This scheme could be used for the preparation of thermal states, in which case one would choose |Ψi to be the Hamiltonian eigenstates and pi = e−βEi /Z, where β = 1/kB T and Z is the partition function. Our approach to the thermalization problem is therefore different from that of Terhal and DiVincenzo, who prepare thermal states by simulating an external bath [36].

7 We assume that each |Ψi can be efficiently specified using some specification |ξi (for example, |Ψi is the ξi th eigenstate of the Hamiltonian). We begin by preparing the state √ i pi |ξi . This can be done using the procedure in Sec. II A if we order the ξi ’s so that they may be thought of as a function on a one-dimensional grid. We then run the entire state-preparation algorithm in superposition, preparing the appropriate |Ψi conditional on the value of the |ξi . This yields the state √ pi |ξi |Ψi , |Ξ =
i

by preparing the appropriate state in separate registers as was done in Sec. II C. Creating |Θ itself can be done in analogy to the preparation of superpositions in Sec.II D. We start by efficiently specifying |Θ using occupation number vectors of the |ΦA,i and the |ΦB,i , namely
i,j

αi,j |nA,i

2nd |nB,j 2nd

⊗ |0A

1st |0B 1st .

the density operator of which is ρΞ = ˆ
i,i′

√

pi pi′ |ξi ξi′ | ⊗ |Ψi Ψi′ | .

Tracing out the specification register, we get the desired density operator ρ = Trξ ρΞ ˆ ˆ √ = pi pi′ |Ψi Ψi′ | Tr |ξi ξi′ |
i,i′

=
i

pi |Ψi Ψi | .

We then complete the state preparation, in superposition, as we did in Sec. II D, treating each register separately. Doing so produces |Θ . There are many circumstances in which the ability to prepare states such as these would be valuable. For instance, in chemical dynamics it is necessary to treat the nuclei and the electrons separately. If we restricted our state preparation to simple product states such as |ΦA,i |ΦB,j , we would get a state in the Born-Oppenheimer approximation, which is often a good approximation to the initial states of reactants participating in chemical reactions. However, as the procedure for preparing |Θ shows, quantum computers could just as easily prepare non–Born-Oppenheimer states in which there is correlation between electronic and nuclear degrees of freedom. Many-particle mixed states can likewise be prepared by following the procedure in Sec. II G separately for each type of particle.
IV. ERRORS AND THE COMPUTATIONAL COST

In practical terms, tracing out the specification register amounts to doing nothing at all. That is, each |Ψi is entangled to a different |ξi , meaning that the |Ψi ’s evolve separately under time evolution, as they would if they were independent members of an ensemble. One can notice that density operators diagonal in the |Φi basis can be prepared more directly. In the previous Sec. II D, we had to “disentangle” the first- and second-quantized states. If we had instead simply traced out the input register, we would have obtained a mixed state diagonal in the |Φi basis.
III. MANY TYPES OF PARTICLES

For the state preparation algorithm to be considered efficient, the time it takes to execute it must scale as a polynomial in the sizes of the input. More precisely, it should scale as a polynomial in l, the number of qubits used to store the wavefunction and m, the number of occupied single-particle orbitals, which is the best descriptor of the total size of the system. In this section, we first show that pre-existing errors are amplified at most linearly by subsequent steps of the algorithm. We then use this fact to obtain the total computational cost of preparing an arbitrary quantum state.
A. Errors

In Sec. II, we outlined an algorithm for the preparation of arbitrary many-particle states (pure or mixed) of a system of identical particles. However, one often wants to consider systems of more than one type of particle, or treat particles of the same kind, but separated in space, as different (the latter approach might be useful, for example, in computing electron transfer matrix elements for large molecules) [37]. We consider the case of two types of particles, with the generalization to more types being clear. One wants to prepare an arbitrary two-particle state |Θ = αi,j |ΦA,i |ΦB,j ,

i,j

where |ΦA,i is a many-particle eigenstate of particles of type A, and |ΦB,j is an eigenstate of particles of type B. Each element |ΦA,i |ΦB,j of this superposition is easily created

Assuming that the quantum gates are executed perfectly— or that the gate errors are corrected using efficient error correction algorithms—there are five sources of error in the state preparation algorithm: 1. Preparation of single-particle eigenstates. Zalka’s method that we adopt in Sec. II A requires evaluation of the integrals (1). We have addressed the computational cost of integral evaluation in Sec. II B, where we show that the procedure can be accomplished in time polynomial in ε−1 if the I wavefunction’s indefinite integral is known or, more generally, if the wavefunction is bounded. The resulting error in the prepared single-particle eigenstate is εφ ≤ lεI /2. 2. Assembly of many-particle eigenstates. Many-particle eigenstates (3) inherit the errors present in the single-particle

8 eigenstates |φi that are used to assemble them. Supposing ˜ that the prepared states φi approximate the true states |φi ˜ ˜ product φj1 . . . φjm suffers an error εΦ = 1 − = 1−
m

with error εφ,i = 1 −

˜ φi φi

, then the prepared Hartree

˜ ˜ φj1 . . . φjm φj1 . . . φjm
m i=1

(1 − εφ,ji )

plete. The second-quantized register, which should be uncomputed during the procedure, should be measured at the end. If |0 . . . is observed, phase estimation will have succeeded. Otherwise, a misidentification will have occurred, and the procedure ought to be repeated. This simple, classical error correction introduces only a constant overhead and ensures that phase estimation does not contribute to the error in the final prepared state. 5. Assembly of mixed states. In the notation of Sec. II G, ˜ if the prepared states Ψi approximate the true states |Ψi with an error εΨ,i = 1 −

≤

i=1

εφ,ji ≤ mεφ ≤ mlεI /2,

where εφ = max εφ,i . Since the total error grows as a polynomial in both the single-state error and the number of occupied states, the assembly of Hartree products amplifies the pre-existing errors only linearly in m. The remaining step, the antisymmetrization of the Hartree product, does not introduce additional errors. 3. Preparation of superpositions. The parallel statepreparation that is used to perform the transformation (6) does not introduce any additional errors with the exception of the possible failures of phase estimation, discussed below. Nevertheless, we should see how pre-existing errors propagate ˜ ˜ through this step. If the prepared state is Ψ = i αi Φi , we see that it suffers an error with respect to the target state εΨ = 1 − = 1− |αi |
2

˜ Ψi Ψi , and assuming per√ fect preparation of the state i pi |ξi , the final prepared ˜ ˜ ˆ pi Ψi Ψi . If we assume that mixed state will be ρ = ˜
i

˜ ˆ ˜ Ψi Ψj = 0 for i = j, then ρ suffers an error [11]
1/2

ερ = 1 − Tr = 1 − Tr = 1−
i

ˆ ρρ ρ ˜ˆ ˆ ˜
1/2

p2 i
i

˜ Ψi

˜ Ψi Ψi

˜ Ψi Ψi

˜ Ψi

pi (1 − εΨ,i )

≤ εΨ ≤ mlεI /2, where εΨ = max εΨ,i . That is, the assembly of mixed states does not magnify the pre-existing errors. Overall, we see that errors introduced in any stage of the state preparation algorithm are not amplified more than polynomially by subsequent stages. The final error in the prepared state is ε = ερ ≤ mlεI /2, meaning that the error scales linearly with the size of the system m and the error of the integration procedure, as well as logarithmically with the grid size 2l .
B. Computational cost

˜ Φi Φi

i

i

|αi |2 (1 − εΦ,i )

≤ εΦ ≤ mlεI /2, where εΦ = maxεΦ,i and where we have assumed that ˜ ˜ Φi Φj = 0 for i = j. In other words, the error in Ψ is limited by the error of its components. 4. Discrimination of states in a superposition. The preparation of superpositions described in Secs. II D-II F and III relies on phase estimation as a means of distinguishing states. Since the eigenenergies will rarely have finite binary expansions, there will be errors introduced at this step. If two phases differ at the nth bit and we perform phase estimation with q = n + p qubits, the probability of an incorrect identification is 1/2(2p − 2), meaning that the success probability will be 1 − εPE provided we implement phase estimation with p = ⌈log (2 + 1/2εPE)⌉ additional qubits [11]. The additional overhead, logarithmic in ε−1 , does not compromise the PE efficiency. The same arguments apply to the phase estimation of eigenvalues of Cn belonging to states in separably degenerate irreducible representations of groups Cn , Cnh , and S2n (see Sec. II E). The symmetry eigenvalues that are useful for states in the other point groups are always ±1, and can be perfectly resolved using phase estimation with a single readout qubit. In addition, failures of state discrimination based on phase estimation can be detected after the state preparation is com-

There are three time-consuming steps in the state preparation algorithm. The first is the evaluation of the integrals (1) and the resulting single-qubit rotations, the second is the phase-estimation that is used to distinguish states in the superposition (see Sec. II D), and the final is the antisymmetrization procedure described in Sec. II C. We characterize the cost of each step in turn. In the previous section, we have seen that the total error of the prepared state will be ε ≤ mlεI /2. Therefore, if we want to ensure a maximum error ε, we must choose εI = 2ε/ml, implying that O(mlε−1 ) time is required for each integration (see Sec. II B). The integration procedure itself is called ml times: for each of the m occupied orbitals, l qubits have to be rotated correctly. Therefore, the total time required for all the qubit rotations is O(m2 l2 ε−1 ). The cost of the phase-estimation procedure that is used to distinguish the eigenstates cannot be given precisely because we have not made any assumptions about the nature of the

9 ˆ Fock operator F other than that it is efficiently simulatable, that is, running in time poly(m, M, l, ∆−1 ) (here, ∆ is the precision at which the simulation needs to be run, i.e., it is half the gap between the closest two eigenstates, which we assumed is not exponentially small). Simulating the entire Hartree-Fock Hamiltonian requires the simulation of the Fock operator acting separately on each particle, meaning that the total simulation requires mpoly(m, M, l, ∆−1 ) time. In addition to this, two quantum Fourier transforms (QFTs) are required on the readout register of the phase estimation. If q qubits are used for the readout (see Sec. IV A.4), the QFTs require O(q 2 ) time. It should be noted that the required q is determined only by needed precision in the phase estimation, and that it does not depend strongly on m, l, or M . Therefore, the cost of the QFTs can be treated as essentially a constant overhead. Furthermore, there is the cost of looking up the state’s energy in the look-up table; a simple binary search requires O(log2 M ) time per register, for a total cost of O(m log2 M ). But this, too, is a negligible cost in comparison to mpoly(m, M, l, ∆−1 ), which we conclude is the asymptotic cost of the eigenstate discrimination portion of the state preparation algorithm. The bottleneck of the antisymmetrization procedure used to produce fermionic states (or the symmetrization for bosonic ones) is the sort that takes state (4) to (5). Sorting register B by a comparison sort requires Ω(m log m) swaps. These swaps must also be performed on each of the corresponding l qubits of register A, for a total cost of Ω(lm log m). For large systems, this expression will be dominated by the scalings of the integral evaluation and the phase estimation. Based on the foregoing, the total computational cost of the state preparation algorithm is O(m2 l2 ε−1 + mpoly(m, M, l, ∆−1 )) = poly(m, M, l, ε−1 , ∆−1 ), an expression polynomial in all the basic descriptors of the system. This allows us to conclude that the algorithm, as described above, is efficient.

V.

CONCLUSION

We have outlined a quantum algorithm for the preparation of physically realistic quantum states on a lattice. In particular, we have gone beyond previous proposals by describing a method for preparing any pure or mixed state of any number of particles. This is achieved by using Zalka’s method for preparing single-particle states and then combining those into manyparticle states. The assembly of many-particle states requires that we be able to distinguish them on a quantum computer, a task that we address using phase estimation. We also provided symmetry-based solutions for degenerate cases, where phase estimation using a single operator is insufficient to distinguish the states. Accidentally degenerate states can be distinguished by adding a perturbation to the system Hamiltonian. Our algorithm is efficient, with a run-time of poly(m, M, l, ε−1 , ∆−1 ), subject only to the requirements that the wavefunction be bounded or that its indefinite integral be known and that the Fock operator be efficiently simulatable.
Acknowledgments

We acknowledge support from the Army Research Office under contract W911NF-07-0304. NJW thanks the Harvard College Research Program and IK the Joyce and Zlatko Balokovi´ Scholarship. c

[1] [2] [3] [4] [5] [6] [7] [8] [9]

[10] [11] [12] [13] [14] [15]

R. Feynman, Inter. J. Theor. Phys. 21, 467 (1982). S. Lloyd, Science 273, 1073 (1996). S. Wiesner, quant-ph/9603028 (1996). C. Zalka, Proc. Roy. Soc. A 454, 313 (1998). D. A. Lidar and H. Wang, Phys. Rev. E 59, 2429 (1999). A. Aspuru-Guzik, A. D. Dutoi, P. J. Love, and M. HeadGordon, Science 309, 1704 (2005). H. Wang, S. Kais, A. Aspuru-Guzik, and M. R. Hoffmann, Phys. Chem. Chem. Phys. 10, 5388 (2008). I. Kassal, S. P. Jordan, P. J. Love, M. Mohseni, and A. AspuruGuzik, Proc. Natl. Acad. Sci. 105, 18681 (2008). B. P. Lanyon, J. D. Whitfield, G. G. Gillet, M. E. Goggin, M. P. Almeida, I. Kassal, J. D. Biamonte, M. Mohseni, B. J. Powell, M. Barbieri, A. Aspuru-Guzik, A. G. White, Submitted (2008). D. Aharonov and A. Ta-Shma, SIAM J. Comput. 37, 47 (2007). M. A. Nielsen and I. L. Chuang, Quantum Computation and Quantum Information (Cambridge University Press, 2000). M. Möttönen, J. J. Vartiainen, V. Bergholm, and M. M. Salomaa, Quant. Inf. Comp. 5, 467 (2005). V. Bergholm, J. J. Vartiainen, M. Möttönen, and M. M. Salomaa, Physical Review A 71, 052330 (2005). H. Wang, S. Ashhab, and F. Nori, 0902.1419 (2009). A. Szabo and N. Ostlund, Modern Quantum Chemistry: Introduction to Advanced Electronic Structure Theory (Dover, 1996).

[16] D. W. Berry, G. Ahokas, R. Cleve, and B. C. Sanders, Comm. Math. Phys. 270, 359 (2007). [17] A. Childs, Ph.D. thesis, Massachusetts Institute of Technology (2004). [18] L. Grover and T. Rudolph, quant-ph/0208112 (2002). [19] P. Kaye and M. Mosca, quant-ph/0407102 (2004). [20] R. Cleve, A. Ekert, C. Macchiavello, and M. Mosca, Proc. Roy. Soc. A 454, 339 (1998). [21] J. M. Borwein and P. B. Borwein, Pi and the AGM: A Study in Analytic Number Theory and Computational Complexity (Wiley-Interscience, 1998). [22] G. S. Fishman, Monte Carlo: Concepts, Algorithms, and Applications (Springer, 1995). [23] A. Peres, Quantum Theory: Concepts and Methods (Springer, 1995). [24] L. K. Grover, in Proceedings of the 30th Annual ACM Symposium on the Theory of Computing (ACM, 1998), pp. 53–62. [25] D. S. Abrams and C. P. Williams, quant-ph/9908083 (1999). [26] G. Brassard, P. Høyer, M. Mosca, and A. Tapp, quantph/0005055 (2000). [27] A. Nayak and F. Wu, in Proceedings of the 31st Annual ACM Symposium on the Theory of Computing (ACM, 1999), pp. 384– 393. [28] E. Novak, J. Complexity 17, 2 (2001). [29] A. N. Soklakov and R. Schack, Phys. Rev. A 73, 012307 (2006).

10
[30] [31] [32] [33] D. S. Abrams and S. Lloyd, Phys. Rev. Lett. 79, 2586 (1997). A. Y. Kitaev, quant-ph/9511026 (1995). D. S. Abrams and S. Lloyd, Phys. Rev. Lett. 83, 5162 (1999). F. A. Cotton, Chemical Applications of Group Theory (Wiley, 1990), 3rd ed. [34] P. R. Bunker and P. Jensen, Molecular Symmetry and Spectroscopy (NRC Research Press, 2006), 2nd ed. [35] C. Zalka, quant-ph/0407140 (2004). [36] B. M. Terhal and D. P. DiVincenzo, Phys. Rev. A 61, 022301 (2000). [37] L. Y. Zhang, R. A. Friesner, and R. B. Murphy, J. Chem. Phys. 107, 450 (1997).