Improved Techniques for Preparing Eigenstates of Fermionic Hamiltonians

Modeling low energy eigenstates of fermionic systems can provide insight into chemical reactions and material properties and is one of the most anticipated applications of quantum computing. We present three techniques for reducing the cost of preparing fermionic Hamiltonian eigenstates using phase estimation. First, we report a polylogarithmic-depth quantum algorithm for antisymmetrizing the initial states required for simulation of fermions in first quantization. This is an exponential improvement over the previous state-of-the-art. Next, we show how to reduce the overhead due to repeated state preparation in phase estimation when the goal is to prepare the ground state to high precision and one has knowledge of an upper bound on the ground state energy that is less than the excited state energy (often the case in quantum chemistry). Finally, we explain how one can perform the time evolution necessary for the phase estimation based preparation of Hamiltonian eigenstates with exactly zero error by using the recently introduced qubitization procedure.


I. INTRODUCTION
One of the most important applications of quantum simulation (and of quantum computing in general) is the Hamiltonian simulation based solution of the electronic structure problem. The ability to accurately model ground states of fermionic systems would have significant implications for many areas of chemistry and materials science and could enable the in silico design of new solar cells, batteries, catalysts, pharmaceuticals, etc. [1,2]. The most rigorous approaches to solving this problem involve using the quantum phase estimation algorithm [3] to project to molecular ground states starting from a classically guessed state [4]. Beyond applications in chemistry, one might want to prepare fermionic eigenstates in order to simulate quantum materials [5] including models of high-temperature superconductivity [6].
In the procedure introduced by Abrams and Lloyd [7], one first initializes the system in some efficient-to-prepare initial state |ϕ which has appreciable support on the desired eigenstate |k of Hamiltonian H. One then uses quantum simulation to construct a unitary operator that approximates time evolution under H. With these ingredients, standard phase estimation techniques invoke controlled application of powers of U (τ ) = e −iHτ . With probability α k = | ϕ|k | 2 , the output is then an estimate of the corresponding eigenvalue E k with standard deviation σ E k = O((τ M ) −1 ), where M is the total number of applications of U (τ ). The synthesis of e −iHτ is typically performed using digital quantum simulation algorithms, such as by Lie-Trotter product formulas [8], truncated Taylor series [9], or quantum signal processing [10]. * Corresponding author: babbush@google.com Since the proposal by Abrams and Lloyd [7], algorithms for time-evolving fermionic systems have improved substantially [11][12][13][14][15][16][17][18]. Innovations that are particularly relevant to this paper include the use of first quantization to reduce spatial overhead [19][20][21][22][23][24] from O(N ) to O(η log N ) where η is number of particles and N η is number of single-particle basis functions (e.g. molecular orbitals or plane waves), and the use of post-Trotter methods to reduce the scaling with time-evolution error from O(poly(1/ )) to O(polylog(1/ )) [24][25][26]. The algorithm of [24] makes use of both of these techniques to enable the most efficient first quantized quantum simulation of electronic structure in the literature.
Unlike second quantized simulations which necessarily scale polynomially in N , first quantized simulation offers the possibility of achieving total gate complexity O(poly(η) polylog(N, 1/ )). This is important because the convergence of basis set discretization error is limited by resolution of the electron-electron cusp [27], which cannot be resolved faster than O(1/N ) using any singleparticle basis expansion. Thus, whereas the cost of refining second quantized simulations to within δ of the continuum basis limit is necessarily O(poly(1/δ)), first quantization offers the possibility of suppressing basis set errors as O(polylog(1/δ)), providing essentially arbitrarily precise representations.
In second quantized simulations of fermions the wavefunction encodes an antisymmetric fermionic system, but the qubit representation of that wavefunction is not necessarily antisymmetric. Thus, in second quantization it is necessary that operators act on the encoded wavefunction in a way that enforces the proper exchange statistics. This is the purpose of second quantized fermion mappings such as those explored in [28][29][30][31][32][33][34]. By contrast, the distinguishing feature of first quantized simulations is that the antisymmetry of the encoded system must be enforced directly in the qubit representation of the wavefunction. This often simplifies the task of Hamiltonian simulation but complicates the initial state preparation.
In first quantization there are typically η different registers of size log N (where η is the number of particles and N is number of spin-orbitals) encoding integers indicating the indices of occupied orbitals. As only η of the N orbitals are occupied, with η log N qubits one can specify an arbitrary configuration. To perform simulations in first quantization, one typically requires that the initial state |ϕ is antisymmetric under the exchange of any two of the η registers. Prior work presented a procedure for preparing such antisymmetric states with gate complexity scaling as O(η 2 ) [35,36].
In Section II we provide a general approach for antisymmetrizing states via sorting networks.
The circuit size is O(η log c η log N ) and the depth is O(log c η log log N ), where the value of c ≥ 1 depends on the choice of sorting network (it can be 1, albeit with a large multiplying factor). In terms of the circuit depth, these results improve exponentially over prior implementations [35,36]. They also improve polynomially on the total number of gates needed. We also discuss an alternative approach, a quantum variant of the Fisher-Yates shuffle, which avoids sorting, and achieves a sizecomplexity of O(η 2 log N ) with lower spatial overhead than the sort-based methods.
Once the initial state |ϕ has been prepared, it typically will not be exactly the ground state desired. In the usual approach, one would perform phase estimation repeatedly until the ground state is obtained, giving an overhead scaling inversely with the initial state overlap. In Section III we propose a strategy for reducing this cost, by initially performing the estimation with only enough precision to eliminate excited states.
In Section IV we explain how qubitization [37] provides a unitary sufficient for phase estimation purposes with exactly zero error (provided a gate set consisting of an entangling gate and arbitrary single-qubit rotations). This improves over proposals to perform the time evolution unitary with post-Trotter methods at cost scaling as O(polylog(1/ )). We expect that a combination of these strategies will enable quantum simulations of fermions similar to the proposal of [24] with substantially fewer T gates than any method suggested in prior literature.

II. EXPONENTIALLY FASTER ANTISYMMETRIZATION
Here we present our algorithm for imposing fermionic exchange symmetry on a sorted, repetition-free quantum array target. Specifically, the result of this procedure is to perform the transformation where π(σ) is the parity of the permutation σ, and we require for the initial state that r p < r p+1 (necessary for this procedure to be unitary). Although we describe the procedure for a single input |r 1 · · · r η , our algorithm may be applied to any superposition of such states. Our approach is a modification of that proposed in Refs. [35,36]; namely, to apply the reverse of a sort to a sorted quantum array. Whereas Refs. [35,36] provide a gate count of O(η 2 (log N ) 2 ), we can report a gate count of O(η log η log N ) and a runtime of O(log η log log N ).
This section proceeds as follows. We begin with a summary of our algorithm. We then explain the reasoning underlying the key step (Step 4) of our algorithm, which is to reverse a sorting operation on target. Next we discuss the choice of sorting algorithm, which we require to be a sorting network. Then, we assess the cost of our algorithm in terms of gate complexity and runtime and we compare this to previous work in Refs. [35,36]. Finally, we discuss the possibility of antisymmetrizing without sorting and propose an alternative, though more costly, algorithm based on the Fisher-Yates shuffle. Our algorithm consists of the following four steps: 1. Prepare seed. Let f be a function chosen so that f (η) ≥ η 2 for all η. We prepare an ancillary register called seed in an even superposition of all possible length-η strings of the numbers 0, 1, . . . , f (η) − 1.
If f (η) is a power of two, preparing seed is easy: simply apply a Hadamard gate to each qubit.
2. Sort seed. Apply a reversible sorting network to seed. Any sorting network can be made reversible by storing the outcome of each comparator in a second ancillary register called record. There are several known sorting networks with polylogarithmic runtime, as we discuss below.
3. Delete collisions from seed. As seed was prepared in a superposition of all length-η strings, it includes strings with repeated entries. As we are imposing fermionic exchange symmetry, these repetitions must be deleted. We therefore measure seed to determine whether a repetition is present, and we accept the result if it is repetition-free. We prove in Appendix A that choosing f (η) ≥ η 2 ensures that the probability of success is greater than 1/2. We further prove that the resulting state of seed is disentangled from record, meaning seed can be discarded after this step.

4.
Apply the reverse of the sort to target. Using the comparator values stored in record, we apply each step of the sorting network in reverse order to the sorted array target. The resulting state of target is an evenly weighted superposition of each possible permutation of the original values.
To ensure the correct phase, we apply a controlledphase gate after each swap.
Step 4 is the key step. Having prepared (in Step 1-Step 3) a record of the in-place swaps needed to sort a symmetrized, collision-free array, we undo each of these swaps in turn on the sorted target. We employ a sorting network, a restricted type of sorting algorithm, because sorting networks have comparisons and swaps at a fixed sequence of locations. By contrast, many common classical sorting algorithms (like heapsort) choose locations depending on the values in the list. This results in accessing registers in a superposition of locations in the corresponding quantum algorithm, incurring a linear overhead. As a result, a quantum heapsort requires O η 2 operations, not O(η). By contrast, no overhead is required for using a fixed sequence of locations. The implementation of sorting networks in quantum algorithms has previously been considered in Refs. [38,39].
Sorting networks are logical circuits that consist of wires carrying values and comparator modules applied to pairs of wires, that compare values and swap them if they are not in the correct order. Wires represent bit strings (integers are stored in binary) in classical sorting networks and qubit strings in their quantum analogues. A classical comparator is a sort on two numbers, which gives the transformation (A, B) → (min(A, B), max(A, B)). A quantum comparator is its reversible version where we record whether the items were already sorted (ancilla state |0 ) or the comparator needed to apply a swap (ancilla state |1 ); see Figure 1.
The positions of comparators are set as a predetermined fixed sequence in advance and therefore cannot depend on the inputs. This makes sorting networks viable candidates for quantum computing. Many of the sorting networks are also highly parallelizable, thus allowing low-depth, often polylogarithmic, performance.
Our algorithm allows for any choice of sorting network. Several common sort algorithms such as the insert and bubble sorts can be represented as sorting networks. However, these algorithms have poor time complexity even after parallelization. More efficient runtime can be achieved, for example, using the bitonic sort [40,41], which is illustrated for 8 inputs in Figure 2. The bitonic sort uses O(η log 2 η) comparators and O(log 2 η) depth, thus achieving an exponential improvement in FIG. 1. The standard notation for a comparator is indicated on the left. Its implementation as a quantum circuit is shown on the right. In the first step, we compare two inputs with values A and B and save the outcome (1 if A > B is true and 0 otherwise) in a single-qubit ancilla. In the second step, conditioned on the value of the ancilla qubit, the values A and B in the two wires are swapped.
Example of a bitonic sort on 8 inputs. The ancillae necessary to record the results as part of implementing each of the comparators are omitted for clarity. Comparators in each dashed box can be applied in parallel for depth reduction.
depth compared to common sorting techniques. Slightly better performance can be obtained using an odd-even mergesort [40]. The asymptotically best sorting networks have depth O(log η) and complexity O(η log η), though there is a large constant which means they are less efficient for realistic η [42,43]. There is also a sorting network with O(η log η) complexity with a better multiplicative constant [44], though its depth is O(η log η) (so it is not logarithmic).
Assuming we use an asymptotically optimal sorting network, the circuit depth for our algorithm is O(log η log log N ) and the gate complexity is The other two steps in our algorithm have smaller cost.
Step 1 has constant depth and O(η log η) complexity.
Step 3 requires O(η) comparisons because only nearestneighbor comparisons need be carried out on seed after sorting. These comparisons can be parallelized over two rounds, with complexity O(η log η) and circuit depth O(log log η). Then the result for any of the registers being equal is computed in a single qubit, which has complexity O(η) and depth O(log η). Thus the complexity of Step 3 is O(η log η) and the total circuit depth is O(log η). We give further details in Appendix B. Thus, our algorithm has an exponential runtime improvement over the proposal in Refs. [35,36]. We also have a quadratic improvement in gate complexity, which is O(η) for our algorithm but O(η 2 ) for Refs. [35,36].
Our runtime is likely optimal for symmetrization, at least in terms of the η scaling. Symmetrization takes a single computational basis state and generates a superposition of η! computational basis states. Each single-qubit operation can increase the number of states in the superposition by at most a factor of two, and two-qubit opera-tions can increase the number of states in the superposition by at most a factor of four. Thus, the number of oneand two-qubit operations is at least log 2 (η!) = O(η log η). In our algorithm we need this number of operations between the registers. If that is true in general, then η operations can be parallelized, resulting in minimum depth O(log η). It is more easily seen that the total number of registers used is optimal. There are O(η log η) ancilla qubits due to the number of steps in the sort, but the number of qubits for the system state we wish to symmetrize is O(η log N ), which is asymptotically larger.
Our quoted asymptotic runtime and gate complexity scalings assume the use of sorting networks that are asymptotically optimal. However, these algorithms have a large constant overhead making it more practical to use an odd-even mergesort, leading to depth O(log 2 η log log N ). Note that is possible to obtain complexity O(η log η log N ) and number of ancilla qubits O(η log η) with a better scaling constant using the sorting network of Ref. [44].
Given that the cost of our algorithm is dictated by the cost of sorting algorithms, it is natural to ask if it is possible to antisymmetrize without sorting. Though the complexity and runtime both turn out to be significantly worse than our sort-based approach, we suggest an alternative antisymmetrization algorithm based on the Fisher-Yates shuffle. The Fisher-Yates shuffle is a method for applying to a length-η target array a permutation chosen uniformly at random using a number of operations scaling as O(η). Our algorithm indexes the positions to be swapped, thereby increasing the complexity to O(η 2 ). Briefly put, our algorithm generates a superposition of states as in Step II of Ref. [36], then uses these as control registers to apply the Fisher-Yates shuffle to the orbital numbers. The complexity is O(η 2 log N ), with a factor of log N due to the size of the registers. We reset the control registers, thereby disentangling them, using O(η log η) ancillae. We provide more details of this approach in Appendix C.
To conclude this section, we have presented an algorithm for antisymmetrizing a sorted, repetition-free quantum register. The dominant cost of our algorithm derives from the choice of sorting network, whose asymptotically optimal gate count complexity and runtime are, respectively, O(η log η log N ) and O(log η log log N ). This constitutes a polynomial improvement in the first case and exponential in the second case over previous work in Refs. [35,36]. As in Ref. [36], our antisymmetrization algorithm constitutes a key step for preparing fermionic wavefunctions in first quantization.

III. FEWER PHASE ESTIMATION REPETITIONS BY PARTIAL EIGENSTATE PROJECTION REJECTION
Once the initial state |ϕ has been prepared, it typically will not be exactly the ground state (or other eigen-state) desired. In the usual approach, one would perform phase estimation repeatedly, in order to obtain the desired eigenstate |k . The number of repetitions needed scales inversely in α k = | ϕ|k | 2 , increasing the complexity. We propose a practical strategy for reducing this cost which is particularly relevant for quantum chemistry. Our approach applies if one seeks to prepare the ground state with knowledge of an upper bound on the ground state energyẼ 0 , together with the promise that E 0 ≤Ẽ 0 < E 1 . With such bounds available, one can reduce costs by restarting the phase estimation procedure as soon as the energy is estimated to be aboveẼ 0 with high probability. That is, one can perform a phase estimation procedure that gradually provides estimates of the phase to greater and greater accuracy, for example as in Ref. [45]. If at any stage the phase is estimated to be aboveẼ 0 with high probability, then the initial state can be discarded and re-prepared.
Performing phase estimation within error typically requires evolution time for the Hamiltonian of 1/ , leading to complexity scaling as 1/ . This means that, if the state is the first excited state, then an estimation error less than E 1 −Ẽ 0 will be sufficient to show that the state is not the ground state. The complexity needed would then scale as 1/(E 1 −Ẽ 0 ). In many cases, the final error required, f , will be considerably less than E 1 −Ẽ 0 , so the majority of the contribution to the complexity comes from measuring the phase with full precision, rather than just rejecting the state as not the ground state.
Given the initial state |ϕ which has initial overlap of α 0 with the ground state, if we restart every time the energy is found to be aboveẼ 0 , then the contribution to the complexity is 1/[α 0 (E 1 −Ẽ 0 )]. There will be an additional contribution to the complexity of 1/ f to obtain the estimate of the ground state energy with the desired accuracy, giving an overall scaling of the complexity of In contrast, if one were to perform the phase estimation with full accuracy every time, then the scaling of the complexity f , the method we propose would essentially eliminate the overhead from α 0 . In cases where α 0 is very small, it would be helpful to apply amplitude amplification. A complication with amplitude amplification is that we would need to choose a particular initial accuracy to perform the estimation. If a lower bound on the excitation energy,Ẽ 1 , is known, then we can choose the initial accuracy to beẼ 1 −Ẽ 0 . The success case would then correspond to not finding that the energy is aboveẼ 0 after performing phase estimation with that precision. Then amplitude amplification can be performed in the usual way, and the overhead for the complexity is 1/ √ α 0 instead of 1/α 0 .
All of this discussion is predicated on the assumption that there are cases where α 0 is small enough to warrant using phase estimation as part of the state prepa-ration process and where a bound meeting the promises ofẼ 0 is readily available. We now discuss why these conditions are anticipated for many problems in quantum chemistry. Most chemistry is understood in terms of mean-field models (e.g. molecular orbital theory, ligand field theory, the periodic table, etc.). Thus, the usual assumption (empirically confirmed for many smaller systems) is that the ground state has reasonable support on the Hartree-Fock state (the typical choice for |ϕ ) [46][47][48][49]. However, this overlap will decrease as a function of both basis size and system size. As a simple example, consider a large system composed of n copies of noninteracting subsystems. If the Hartree-Fock solution for the subsystem has overlap α 0 , then the Hartree-Fock solution for the larger system has overlap of exactly α n 0 , which is exponentially small in n.
It is literally plain-to-see that the electronic ground state of molecules is often protected by a large gap. The color of many molecules and materials is the signature of an electronic excitation from the ground state to first excited state upon absorption of a photon in the visible range (around 0.7 Hartree); many clear organics have even larger gaps in the UV spectrum. Visible spectrum E 1 −E 0 gaps are roughly a hundred times larger than the typical target accuracy of f = 0.0016 Hartree ("chemical accuracy") 1 . Furthermore, in many cases the first excited state is perfectly orthogonal to the Hartree-Fock state for symmetry reasons (e.g. due to the ground state being a spin singlet and the excited state being a spin triplet). Thus, the gap of interest is really For most problems in quantum chemistry a variety of scalable classical methods are accurate enough to compute upper bounds on the ground state energyẼ 0 such that E 0 ≤Ẽ 0 < E * , but not accurate enough to obtain chemical accuracy (which would require quantum computers). Classical methods usually produce upper bounds when based on the variational principle. Examples include mean-field and Configuration Interaction Singles and Doubles (CISD) methods [50].
As a concrete example, consider a calculation on the water molecule in its equilibrium geometry (bond angle of 104.5 • , bond length of 0.9584Å) in the minimal (STO-3G) basis set performed using OpenFermion [51] and Psi4 [52]. Hartree, which is about 370 times f . Thus, using our strategy, for α 0 > 0.003 there is very little overhead due to the initial state |ϕ not being the exact ground state. In the most extreme case for this example, that represents a speedup by a factor of more than two orders of magnitude. However, in some cases the ground state overlap might be high enough that this technique provides only a modest advantage. While the Hartree-Fock state overlap in this small basis example is α 0 = 0.972, as the system size and basis size grow we expect this overlap will decrease (as argued earlier).
Another way to cause the overlap to decrease is to deviate from equilibrium geometries [46,47]. For example, we consider this same system (water in the minimal basis) when we stretch the bond lengths to 2.25× their normal lengths. In this case, E 0 = −74.7505 Hartree, E * = −74.6394 Hartree, and α 0 = 0.107. The CISD solution provides an upper boundẼ 0 = −74.7248. In this case, E * −Ẽ 0 ≈ 0.085 Hartree, about 50 times f . Since α 0 > 0.02, here we speed up state preparation by roughly a factor of α −1 0 (more than an order of magnitude).

IV. PHASE ESTIMATION UNITARIES WITHOUT APPROXIMATION
Normally, the phase estimation would be performed by Hamiltonian simulation. That introduces two difficulties: first, there is error introduced by the Hamiltonian simulation that needs to be taken into account in bounding the overall error, and second, there can be ambiguities in the phase that require simulation of the Hamiltonian over very short times to eliminate.
These problems can be eliminated if one were to use Hamiltonian simulation via a quantum walk, as in Refs. [53,54]. There, steps of a quantum walk can be performed exactly, which have eigenvalues related to the eigenvalues of the Hamiltonian. Specifically, the eigenvalues are of the form ±e ±i arcsin(E k /λ) . Instead of using Hamiltonian simulation, it is possible to simply perform phase estimation on the steps of that quantum walk, and invert the function to find the eigenvalues of the Hamiltonian. That eliminates any error due to Hamiltonian simulation. Moreover, the possible range of eigenvalues of the Hamiltonian is automatically limited, which eliminates the problem with ambiguities.
The quantum walk of Ref. [54] does not appear to be appropriate for quantum chemistry, because it requires an efficient method of calculating matrix entries of the Hamiltonian. That is not available for the Hamiltonians of quantum chemistry, but they can be expressed as sums of unitaries, as for example discussed in Ref. [25]. It turns out that the method called qubitization [37] allows one to take a Hamiltonian given by a sum of unitaries, and construct a new operation with exactly the same functional dependence on the eigenvalues of the Hamiltonian as for the quantum walk in Refs. [53,54].
Next, we summarize how qubitization works [37]. One assumes black-box access to a signal oracle V that encodes H in the form: where |0 a is in general a multi-qubit ancilla state in the computational basis, 1 1 s is the identity gate on the system register and λ ≥ H is a normalization constant. For Hamiltonians given by a sum of unitaries, one constructs where A is an operator for state preparation acting as For U that is Hermitian, we can simply take V = U . This is the case for any local Hamiltonian that can be written as a weighted sum of tensor products of Pauli operators, since tensor products of Pauli operators are both unitary and Hermitian. More general strategies for representing Hamiltonians as linear combinations of unitaries, as in Eq. (4), are discussed in [55], related to the sparse decompositions first described in [56]. If U is not Hermitian, then we may construct a Hermitian V as where |± = 1 √ 2 (|0 ± |1 ). The multiqubit ancilla labelled "a" would then include this additional qubit, as well as the ancilla used for the control for SELECT-U. In either case we can then construct a unitary operator called the qubiterate as follows: The qubiterate transforms each eigenstate |k of H as where |0k ⊥ as has no support on |0 a . Thus, W performs rotation between two orthogonal states |0 a |k s and |0k ⊥ as . Restricted to this subspace, the qubiterate may be diagonalized as W |±k as = ∓e ∓i arcsin(E k /λ) |±k as (12) |±k as = 1 √ 2 |0 a |k s ± |0k ⊥ as .
This spectrum is exact, and identical to that for the quantum walk in Refs. [53,54]. This procedure is also simple, requiring only two queries to U and a number of gates to implement the controlled-Z operator (2 |0 0| a ⊗ 1 1 s − 1 1) scaling linearly in the number of controls. We may replace the time evolution operator with the qubiterate W in phase estimation, and phase estimation will provide an estimate of arcsin(E k /λ) or π − arcsin (E k /λ). In either case taking the sine gives an estimate of E k /λ, so it is not necessary to distinguish the cases. Any problems with phase ambiguity are eliminated, because performing the sine of the estimated phase of W yields an unambiguous estimate for E k . Note also that λ ≥ H implies that |E k /λ| ≤ 1.
More generally, any unitary operation e if (H) that has eigenvalues related to those of the Hamiltonian would work so long as the function f (·) : R → (−π, π) is known in advance and invertible. One may perform phase estimation to obtain a classical estimate of f (E k ), then invert the function to estimate E k . To first order, the error of the estimate would then propagate like In our example, with standard deviation σ phase in the phase estimate of W , the error in the estimate is Obtaining uncertainty for the phase of W requires applying W a number of times scaling as 1/ . Hence, obtaining uncertainty for E k requires applying W a number of times scaling as λ/ . For Hamiltonians given by sums of unitaries, as in chemistry, each application of W uses O(1) applications of state preparations and SELECT-U operations. In terms of these operations, the complexities of Section III have multiplying factors of λ.

V. CONCLUSION
We have described three techniques which we expect will be practical and useful for the quantum simulation of fermionic systems. Our first technique provides an exponentially faster method for antisymmetrizing configuration states, a necessary step for simulating fermions in first quantization. We expect that in virtually all circumstances the gate complexity of this algorithm will be nearly trivial compared to the cost of the subsequent phase estimation. Then, we showed that when one has knowledge of an upper bound on the ground state energy that is separated from the first excited state energy, one can prepare ground states using phase estimation with lower cost. We discussed why this situation is anticipated for many problems in chemistry and provided numerics for a situation in which this trick reduced the gate complexity of preparing the ground state of molecular water by more than an order of magnitude. Finally, we explained how qubitization [37] provides a unitary that can be used for phase estimation without introducing the additional error inherent in Hamiltonian simulation.
We expect that these techniques will be useful in a variety of contexts within quantum simulation. In particular, we anticipate that the combination of the three techniques will enable exceptionally efficient quantum simulations of chemistry based on methods similar to those proposed in [24]. While specific gate counts will be the subject of a future work, we conjecture that such techniques will enable simulations of systems with roughly a hundred electrons on a million point grid with fewer than a billion T gates. With such low T counts, simulations such as the mechanism of Nitrogen fixation by ferredoxin, explored for quantum simulation in [57], should be practical to implement within the surface code in a reasonable amount of time with fewer than a few million physical qubits and error rates just beyond threshold. This statement is based on the time and space complexity of magic state distillation (usually the bottleneck for the surface code) estimated for superconducting qubit architectures in [58], and in particular, the assumption that with a bil-lion T gates or fewer one can reasonably perform state distillation in series using only a single T factory. state gaps. We note that arguments relating the colors of compounds to their molecular excited state gaps in the context of quantum computing have been popular for many years, likely due to comments of Alán Aspuru-Guzik. DWB is funded by an Australian Research Council Discovery Project (Grant No. DP160102426).

AUTHOR CONTRIBUTIONS
DWB proposed the algorithms of Section II and the basic idea behind Section III as solutions to issues raised by RB. MK, AS and YRS worked out and wrote up the details of Section II and associated appendices. RB connected developments to chemistry simulation, conducted numerics, and wrote Section III with input from DWB. Based on discussions with NW, GHL suggested the basic idea of Section IV. CG helped to improve the gate complexity of our comparator circuits. Remaining aspects of the paper were written by RB and DWB with assistance from MK, AS and YRS.
In this appendix, we explain the most difficult-tounderstand step of our algorithm: the step in which we delete collisions from seed. There are two important points that require explanation. First, we have to show that the probability of failure is small. Second, we have to show that the resulting state of seed is disentangled from record, as we wish to uncompute record during the final step of our algorithm.
To explain these two points, we begin with an analysis of the state of seed after Step 1. The state of seed is We can decompose the state space of seed into two orthogonal subspaces: the 'repetition-free' subspace and its orthogonal complement. If we project the state of seed onto the repetition-free subspace, we obtain the unnormalized vector The norm of this vector is which is equal to 1 − C(f (η), η) in the terminology of Proposition A.1 in [59]. We sort the register in Step 2 before detecting repetitions in Step 3, because then it is only necessary to check adjacent registers. The probability of repetitions in unaffected by the sort, because it is unitary and does not affect whether there are repetitions. Therefore the probability of failure (detection of a repetition) in Step 3 is equal to C(f (η), η). Using Proposition A.1 in [59], the probability of failure is bounded as which is less than 1/2 for f (η) ≥ η 2 . The repetition-free outcome can therefore be achieved after fewer than two attempts on average. One can improve the success probability by using a larger function f or by using amplitude amplification. We now show that seed ⊗ record is in an unentangled state after Step 3. After Step 1, the state of seed ⊗ record projected to the repetition-free subspace can be represented (up to normalization) as 0≤ 0<...< η−1<f (η) σ∈Sη |σ ( 0 , . . . , η−1 ) seed |ι record .
(A9) This is a product state. Therefore, seed can be discarded after Step 3 without affecting record.
Appendix B: Quantum Sorting

Resource Analysis of Quantum Sorting Networks
In this section we expand on the resource analysis of implementing quantum sorting networks. We also illustrate that for small number of inputs to be sorted (up to η = 20), concrete bounds have been derived for optimized circuit depth as well as the number of comparators. This may be of interest and useful for implementing quantum simulations of small molecules, also in view of the observation that η ≈ 20 is nearly reaching a number of electrons for where classical simulations become intractable.
Optimizing sorting networks for small inputs is an active research area in parallel programming. Knuth [60] and later Codish et al. [61] gave networks for sorting up to 17 numbers that were later shown to be optimal in depth, and up to η ≤ 10 also optimal in the number of comparators. Optimizations for up to 20 inputs have recently been achieved, see Table 1 in [61]. In such optimizations one typically distinguishes between the optimal depth problem and the problem of minimizing the overall number of comparators. For illustration, the best known sorting networks for 20 numbers require depth 11 and 92 comparators, with lower bounds reported as 10 and 73 respectively. Efficient sorting networks can be produced by in-place merging of sorting networks with smaller sizes. However, this procedure necessarily produces some overhead.
For our resource analysis we assume that the quantum sorting network has η wires, where each wire represents a quantum register of length d (i.e., consists of d qubits). The resource requirement for implementing the quantum sort is obtained by taking the (classical) sorting network depth or the overall number of comparators involved and multiplying it by the corresponding resources needed to construct a comparator. As explained above, the latter requires one query to a comparison oracle, whose circuit implementation and complexity are provided in Appendix B 2, and a conditional swap applied to the compared registers of size d controlled by the single-qubit ancilla holding the result of the comparison.
The construction of the comparison oracle as well as the implementation of the conditional swaps both yield a network consisting predominantly of Toffoli, Not and CNot gates requiring O(d) elementary gate operations but only O(log d) circuit depth. Indeed, as shown in Appendix B 2, the comparison oracle can be implemented such that the operations can mostly be performed in parallel with only O(log d) circuit depth.
When implementing conditional swaps on two registers of size d as part of a comparator, all elementary swaps between the corresponding qubits of these registers must be controlled by the very same ancilla qubit, namely the one encoding the result of the comparison oracle. This suggests having to perform all the controlled swaps in sequence, as they all are to be controlled by the same qubit, which would imply depth scaling O(d) rather than O(log d). Yet the conditional swaps can also be parallelized. This can be achieved by first copying the bit of the ancilla holding the result of the comparison to Taking d = log N (the largest registers used in Step 4 of our sort-based antisymmetrization algorithm), conducting the quantum bitonic sort, for instance, thus requires O(η log 2 (η) log N ) elementary gates but only O(log 2 (η) log log N ) circuit depth, while the overall worst-case ancillary space overhead amounts to O(η log 2 (η) log N ).

Comparison Oracle
Here we describe how to implement reversibly the comparison of the value held in one register with the value carried by a second equally-sized register, and store the result (larger or not) in a single-qubit ancilla. We term the corresponding unitary process a 'comparison oracle'. We need to use it for implementing the comparator modules of quantum sorting networks as well as in our antisymmetrization approach based on the quantum Fisher-Yates shuffle. We first explain a naive method for comparison with depth linear in the length of the involved registers. In the second step we then convert this prototype into an algorithm with depth logarithmic in the register length using a divide and conquer approach.
Let A and B denote the two equally sized registers to be compared, and A and B the values held by these two registers. To determine whether A > B or A < B or A = B, we compare the registers in a bit-by-bit fashion, starting with their most significant bits and going down to their least significant bits. At the very first occurrence of an i such that To achieve a reversible comparison, we employ two ancillary registers, each consisting of d qubits, and each initialized to state |0 ⊗d , respectively. We denote them by While this algorithm works, it has the drawback that the bitwise comparison is conducted sequentially, which results in circuit-depth scaling O(d). It also uses more ancilla qubits than necessary. We can improve upon this. We can reduce the number of ancilla qubits by reusing some input bits as output bits, and we can achieve a depth scaling of O(log d) by parallelizing the bitwise comparison. To introduce a parallelization, observe the following. Let us split the register A into two parts: A 1 consisting of the first approximately d/2 bits and A 2 consisting of the remaining approximately d/2 bits. Split register B in the very same way into subregisters B 1 and B 2 . We can then determine which number is larger (or whether both are equal) for each pair (A 1 , B 1 ) and (A 2 , B 2 ) separately in parallel (using the method described above) and record the results of the two com- parisons in ancilla registers (A 1 , B 1 ), (A 2 , B 2 ). The least significant bits of these four ancilla registers can then be used to deduce whether A > B or A < B or A = B with just a single bitwise comparison. Thus, we effectively halved the depth by dividing the problem into smaller problems and merging them afterwards. We now explain a bottom-up implementation.
Instead of comparing the whole registers A and B, our parallelized algorithm slices A and B into pairs of bitsthe first slice contains A[0] and A [1], the second slice consists of A [2] and A [3], etc., and in the very same way for B.
The key step takes the corresponding slices of A and B and overwrites the second bit of each slice with the outcome of the comparison. The first bit of each slice is then ignored, so that the comparison results stored in the second bits become the next layer on which bitwise comparisons are performed. We denote the i'th bit forming the registers of the j th layer by A j [i] and B j [i]. The original registers A and B correspond to j = 0: A 0 ≡ A and B 0 ≡ B. The part of the circuit that implements a single bitwise comparison is depicted in Figure 3. We denote the corresponding transformation by 'Compare2', i.e. At each step, comparisons of the pairs of the original arrays can be performed in parallel, and produce two new arrays with approximately half the size of the original ones to record the results. Thus, at each step we approximately halve the size of the problem, while using a constant depth for computing the results. The basic idea is illustrated in Figure 4. This procedure is repeated for log d steps 3 until registers A fin := A log d and B fin := B log d both of size 1 have been prepared.
This parallelized algorithm is perfectly suited for comparing arrays whose length d is a power of 2. If d is not 3 All logarithms are taken to the base 2. 3. A circuit that implements Compare2, taking a pair of 2-bit integers and outputting a pair of single bits while preserving inequalities. The input pair is (x, y) = (x0 +2x1, y0 +2y1). The output pair is (x , y ) and will satisfy sign(x − y ) = sign(x − y). Output qubits marked "temp" store values that are not needed, and are kept until a later uncompute step where the inputs are restored. Each Fredkin gate within the circuit can be computed using 4 T gates and (by storing an ancilla not shown) later uncomputed using 0 T gates [62,63]. a power of 2, we can either pad A and B with 0s prior to their most significant bits without altering the result, or introduce comparison of single bits (using only the first two gates from the circuit in Figure 3 with targets on A j+1 and B j+1 registers respectively).
Formally, we can express our comparison algorithm as follows, here assuming d to be a power of 2: The key feature of this algorithm is that all the operations of the inner loop can be performed in parallel. Since one application of Compare2 requires only constant depth and constant number of operations, our comparison algorithm requires only depth O(log d).
x 5. A circuit that determines if two bits are equal, ascending, or descending. When the comparison is no longer needed, the results are uncomputed by applying the circuit in reverse order.
Our comparison algorithm constructed above can indeed be used to output a result that distinguishes between A > B, A < B and A = B. Observe that its reversible execution results in the ancillary single-qubit registers A fin and B fin generated in the very last step of the algorithm holding information about which number is larger or whether they are equal. Indeed, The three cases are separated into three control qubits by using the circuit shown in Figure 5. These individual control qubits can be used to control further conditional operations that depend on the result of the comparison.
For the purpose in our applications (comparator modules of quantum sorting networks or quantum Fisher-Yates shuffle), we only need to condition on whether A > B is true or false. Thus, we only need the first operation from the circuit in Figure 5 which takes a single qubit initialized to |0 and transforms it into the output of the comparison oracle. After the output bit has been produced, we must reverse the complete comparison algorithm (invert the corresponding unitary process), thereby uncomputing all the ancillary registers that have been generated along this reversible process and restoring the input registers A and B. A and B (holding values A and B) and a single-qubit ancilla q initialized to |0 . It reversibly computes whether A > B is true or false by executing the parallelized comparison process presented above. It copies the result (which is stored in A fin ) to ancilla q. Here the green boxes identify the array entry that has been swapped at each stage of the shuffle. Observe that the green boxes also label the largest value in the array truncated to position k.

Appendix C: Symmetrization Using The Quantum
Fisher-Yates Shuffle In this appendix we present an alternative approach for antisymmetrization that is not based on sorting, yielding a size-and depth-complexity O(η 2 log N ), but with a lower spatial overhead than the sort-based method. Our alternative symmetrization method uses a quantum variant of the well-known Fisher-Yates shuffle, which applies a permutation chosen uniformly at random to an input array input of length η. A standard form of the algorithm is given in [64].
Swap positions k and of input. end for The basic idea is illustrated in Figure 6 for η = 4.
There are two key steps that turn the Fisher-Yates shuffle into a quantum algorithm. First, our quantum implementation of the shuffle replaces the random selection of swaps with a superposition of all possible swaps. To achieve this superposition, the random variable is replaced by an equal-weight superposition 1 √ k+1 k =0 | in an ancillary register (called choice). At each step of the quantum Fisher-Yates shuffle, the choice register must begin and end in a fiducial initial state.
In order to reset the choice register, we introduce an additional index register, which initially contains the integers 0, . . . , η − 1. We shuffle both the length-η input register and the index register, and the simple form of index enables us to easily reset choice. The resulting state of the joint input ⊗ index register is still highly entangled; however, provided input was initially sorted in ascending order, we can disentangle index from input.
Our quantum Fisher-Yates shuffle consists of the following steps:

5.
Repeat. Increment k by one. If k < η, go to Step 2. Otherwise, proceed to the next step.
6. Disentangle index from input. For each k = = 0, 1, . . . , η − 1, subtract 1 from position of index if the element at position k in input is greater than the element at position in input. The resulting state of index is |0, 0, . . . , 0 , which is disentangled from input.
We present an overview of the algorithm in Figure 7. At the highest level, depicted in Figure 7a, we apply an initialization procedure to index, then η − 1 'Fisher-Yates' blocks (F Y k for k = 1, . . . , η −1), and finally a disentangling ('Detangle') procedure on index and input. Following the Detangle procedure, the ancillary registers choice and index are reset to their initial all-zero states and the input register has been symmetrized. In each Fisher-Yates block, depicted in Figure 7b, we apply the preparation operator Π k to choice, apply selected swaps on choice+index and choice+input, then apply a phase conditioned on choice to input, and finally reset the choice register. Preparing and resetting choice as well as executing swap are therefore part of each Fisher-Yates block and are thus each applied a total of η − 1 times (for each of k = 1, . . . , η − 1). Their gate counts and circuits depths must thus be multiplied by (η − 1). Disentangling index and input is the most expensive step, but it is executed only once, so it contributes only an additive cost to the overall resource requirement. In what follows, we explain each step of the algorithm and justify their corresponding resource contributions, which are briefly summarized here: Step 1 requires O(η log η) gates but has a negligible depth O(1).
Step 2 requires O(η) gates and has the same depth complexity.
Step 3 requires O(η log N ) gates and has also depth O(η log N ).
Step 4 requires O(η log η) gates but has only depth O(log η). As Step 2 to Step 4 are repeated η − 1 times, the total gate count before Step 6 is O(η 2 log N ). Finally, Step 6 requires O(η 2 log N ) gates and has depth O η 2 [log log N + log η] . Thus the total gate count of the quantum Fisher-Yates shuffle is O(η 2 log N ). Because most of the gates need to be performed sequentially, the overall circuit depth of the algorithm is also O(η 2 log N ).
Our complexity analysis is given in terms of elementary gate operations, a term we use loosely. Generally speaking, we treat all single-qubit gates as elementary and we allow up to two controls for free on each singlequbit gate. This definition of elementary gates includes several standard universal gate sets such as Clifford+T and Hadamard+Toffoli. A more restrictive choice of elementary gate set only introduces somewhat larger constant factors in most of the procedure. The exception is the application of Π k in the first step of F Y k , where we require the ability to perform controlled single-qubit rotations of angle arcsin +1 , where = 1, . . . , k. The Solovay-Kitaev theorem implies a gate-count overhead that grows polylogarithmically in the inverse of the error tolerance. We now proceed by analyzing each step to the quantum Fisher-Yates shuffle.

Initialization
The first step is to initialize choice in the state |0 ⊗η . This is assumed to have zero cost. The index register is set to the state |0, 1, . . . , η − 1 that represents the positions of the entries of input in ascending order. Because each of the η entries in index must be capable of storing any of the values 0, 1, . . . , η − 1, the size of index is η log η qubits. This step requires O(η log η) single-qubit gates that can be applied in parallel with circuit depth O(1).

Fisher-Yates Blocks
Each Fisher-Yates block has three stages: prepare choice, executed selected swaps, and reset choice. The exact steps depend on the encoding of the choice register; in particular, whether it is binary or unary.
We elect the conceptually simplest encoding of choice, which is a kind of unary encoding. We use η qubits (labelled 0, 1, . . . , η), define where X is the Pauli X applied to the qubit labelled . An advantage of our encoding for choice is that the selected swaps require only single-qubit controls. An obvious disadvantage is the unnecessary space overhead. Although one can save space with a binary encoding, the resulting operations become somewhat more complicated and hence come at an increased time cost. Our choice of encoding is made for simplicity.

a. Prepare choice
Our preparation procedure has two stages. First, we prepare an alternative unary encoding of the state which we name for its resemblance to the W-state 1 √ 3 (|001 + |010 + |100 ). Second, we translate the alternative unary encoding to our desired encoding. For a summary of the procedure, see Figure 8.
Next, we explain how to prepare |W k in the alternative unary encoding. The alternative encoding is We can prepare |W k in this encoding with a cascade of controlled rotations of the form Explicitly: Apply X to qubit 0. Apply R k to qubit 1. for = 1, . . . , k − 1 do Apply R k− controlled on qubit to qubit + 1. end for This is a total of k + 1 gates, k = 1 of which are applied sequentially.
Next we explain how to translate to the desired encoding. This is a simple procedure: for = k, . . . , 1 do Apply Not controlled on qubit to qubit − 1. end for The total number of CNot gates is k, and they must be applied in sequence. Thus the total gate count (and time-complexity) for preparing choice is O(k) = O(η).

b. Selected Swap
We need to implement selected swaps of the form where the Swap(c, k) operator acts on either target = index or target = input. Here the state of the choice register selects which entry in the target array is to be swapped with entry k. Our unary encoding of the choice register allows for a simple implementation of SelSwap; see Figure 9.
Observe that only the first k + 1 subregisters are involved of each choice, index and input, respectively. Also observe that, for each i = 0, 1, . . . , k, index[i] is of size log η whereas input[i] is of size log N . Hence, the circuit actually consists of k log η + k log N ordinary 3-qubit controlled-Swap gates that for the most part must be executed sequentially. As η ≤ N , we report O (η log N ) for both gate count and depth.  9. Implementation of the two selected swaps SelSwap k as part of F Y k , with the unary-encoded choice as the control register and index and input as target registers, respectively. As each wire of the target registers stand for several qubits, each controlled-Swap is to be interepreted as many bitwise controlled-Swaps.

c. Applying the controlled-phase
Applying the controlled-phase gate is straightforward. We select a target qubit in the input register -it does not matter which. Then, for each = 0, 1, . . . , k − 1, we apply a phase gate controlled on position of choice to the target qubit. The result is that input has picked up a phase of (−1) if choice specified a value strictly less than k. The total number of gates is k = O(η), while the depth can be made O(1).

d. Resetting choice register
The reason we execute swaps on both index and input is to enable reversible erasure of choice at the end of each Fisher-Yates block. This is done by scanning index to find out which value of k was encoded into choice. In general, we know that step k sends the value k to position of index, where is specified by the choice register. We thus erase choice by applying a Not operation to choice[ ] if index[ ] = k. This can be expressed as a multi-controlled-Not, as illustrated by an example in Figure 10. The control sequence of the multi-controlled-Not is a binary encoding of the value k. For compiling multiple controls, see Figure 4.10 in [65]. Each log η -fold-controlled-Not can be decomposed into a network of O(log η) gates (predominantly Toffolis) with depth O(log η). Because the k + 1 multifold-controlled-Nots (for ≤ k ≤ η − 1) can all be executed in parallel, resetting choice register thus requires a circuit with O(η log η) gates but only O(log η) depth.

Disentangling index from input
The last task is to clean up and disentangle index from input by resetting the former to the original state |0 ⊗η log η while leaving the latter in the desired antisymmetrized superposition. This can be achieved as follows.
We compare the value carried by each of the η subregisters input[ ] (labeled by position index = 0, 1, . . . , η−1) with the value of each other subregister input[ ] ( = ), thus requiring η(η − 1) comparisons in total. Note that these subregisters of input have all size log N . Each time the value held in input[ ] is larger than the value carried by any other of the remaining η − 1 subregisters input[ ], we decrement the value of the corresponding th subregister index[ ] of index by 1. In cases in which the value carried by input[ ] is smaller than input[ ], we do not decrement the value of index[ ]. After accomplishing all the η(η −1) comparisons within the input register and controlled decrements, we have reset the index register state to |0 ⊗η log η while leaving the input register in the antisymmetrized superposition state.
Each comparison between the values of two subregisters of input (each of size log N ) can be performed using the comparison oracle introduced in Appendix B 2.
The oracle's output is then used to control the 'decrement by 1' operation, after which the oracle is used again to uncompute the ancilla holding its result. The comparison oracle has been shown to require O(log N ) gates but to have only circuit depth O(log log N ).
Decrementing the value of the log η -sized index subregister index[ ] (for any = 0, 1, . . . , η − 1) by the value 1 can be achieved by a circuit depicted in Figure 11.