Projected Cooling Algorithm for Quantum Computation

In the current era of noisy quantum devices, there is a need for quantum algorithms that are efficient and robust against noise. Towards this end, we introduce the projected cooling algorithm for quantum computation. The projected cooling algorithm is able to construct the localized ground state of any Hamiltonian with a translationally-invariant kinetic energy. The method can be viewed as the quantum analog of evaporative cooling. We start with an initial state with support over a compact region of a large volume. We then drive the excited quantum states to disperse and measure the remaining portion of the wave function left behind. The method can be used in concert with other techniques such as variational methods and adiabatic evolution to achieve better performance than existing approaches for the same number of quantum gates per qubit. For the nontrivial examples we consider here, the improvement is substantial. The only additional resource required is performing the operations in a volume significantly larger than the size of the localized state.

The quantum many-body problem is a profound challenge that pervades nearly all branches of quantum physics. While our own interests are in ab initio methods for nuclear structure and reactions [1,2], the difficulties that arise at strong coupling are similar for other fields, from strongly-correlated electrons to degenerate atomic gases, from quantum spin models to interacting relativistic quantum fields. For computational methods that construct quantum wave functions, the problem is the exponential growth of dimensionality with system size. For computational methods based on Monte Carlo simulations, the obstacle is the Monte Carlo sign problem, where positive and negative contributions cancel to produce exponentially small signals. Quantum computing has emerged as a new computational paradigm that offers the possibility of evading both of these fundamental problems. By allowing for arbitrary quantum superpositions of products of qubits, one can store exponentially more information than classical bits without the need for stochastic sampling. Furthermore, qubits naturally evolve with real-time dynamics, a feature that is extremely challenging to reproduce with classical computing. With these opportunities, however, there are enormous challenges to realizing their promise. All of the currently available quantum computing devices suffer from short decoherence times and significant readout errors. Therefore it is necessary to develop quantum algorithms that are efficient and robust against noise. To address this need, in this work we introduce a method called the projected cooling algorithm.
We start by reviewing several other algorithms for constructing the ground state of a quantum Hamiltonian. One approach is quantum phase estimation [3,4], which starts with some initial state |ψ I and performs successive operations of the time evolution operator U for some time duration ∆t. By constructing superpositions of the states U n |ψ I , one can find approximate eigenvectors of U which are then also eigenvectors of the Hamiltonian H. The main challenge with this approach is that one needs coherent evolution over long time intervals that scale inversely with the desired energy resolution.
A second approach is based on quantum variational methods [5,6]. Here one starts with some initial state |ψ I and computes expectation values of H using wave functions of the form U (θ i ) |ψ I , where U is a unitary operator with control parameters θ i . The ground state wave function is achieved when the expectation value of H is at its global minimum. As with variational approaches used in classical computing, this method is very efficient for constructing approximations to ground state wave functions. However there is no known systematic route for improving the approximation to arbitrary accuracy for general many-body systems.
A third approach is the quantum adiabatic evolution algorithm [7,8]. The adiabatic theorem states that if we start in an eigenstate of some time-dependent Hamiltonian H(t), then we remain in an eigenstate of H(t) so long as the rate of change is sufficiently slow. The strategy is to start with some simple Hamiltonian H I with ground state |ψ I that is easy to prepare, and then slowly evolve to the desired Hamiltonian H with ground state |ψ 0 . The difficulty is that for nontrivial problems, the wave functions |ψ I and |ψ 0 are usually very different, and the time evolution in H(t) must be done very slowly to avoid large errors caused by traversals across avoided level crossings. In order to mitigate some of these issues, variational methods have been developed based on approximations to the adiabatic evolution method [9].
A fourth approach is the recently developed spectral comb method [10]. Here the system of interest is entangled with an auxiliary system of qubits with its own set of tunable energy levels. By varying the properties of the auxiliary system, avoided level crossings are produced that transfer energy out of the system of interest. As with the other approaches, it provides an efficient approximation to the ground state.
However, improving the scaling properties for accurate calculations of larger systems will require further development.
The projected cooling algorithm is a new approach that can be regarded as the quantum analog of evaporative cooling. Rather than evaporating hot gas molecules, we start with some initial state |ψ I with support over a compact region, ρ, and drive excited quantum states to disperse away from ρ. We then measure the remaining portion of the wave function left behind. The algorithm is able to construct the localized ground state of any Hamiltonian with a translationally-invariant kinetic energy. To illustrate, in the following we consider three different examples which we call Models 1A, 1B, and 2. We start with a Hamiltonian defined on a one-dimensional chain of 2L + 1 qubits at sites n = −L, · · · L. We take the vacuum to be the tensor product state where all qubits are |0 , and from this vacuum state we can define particle excitations. In the one-particle subspace, we let |[n] be the tensor product state where qubit n is |1 and the rest are |0 . In the one-particle space, our Hamiltonian is defined as where the kinetic energy term K n ′ ,n is δ n ′ ,n − 1 2 δ n ′ ,n+1 − 1 2 δ n ′ ,n−1 , and V n is the single-particle potential energy on site n. For the first model we consider, Model 1A, we take the interaction term to be V n = −δ 0,n .
We define the compact region ρ to correspond with the qubits n = −R, · · · R where R ≪ L. We define P to be the projection operator that projects onto the subspace where all particle excitations are contained entirely in ρ. Therefore P |[n] = 0 for |n| > R, and P |[n] = |[n] for |n| ≤ R.
Let |ψ 0 be the ground state of H. For our Model 1A Hamiltonian, |ψ 0 is a localized bound state and is the only bound state of H. All of the other states are continuum states that extend to infinity in the limit L → ∞. Let U (t) = e −iHt be the time evolution operator. We are using dimensionless units for all quantities and taking = 1. The key result underpinning the projected cooling method is the fact that in the large volume limit L → ∞, the projected time evolution operator P U (t)P has a stable fixed point proportional to P |ψ 0 . As the time evolution operator U (t) acts on P |ψ I , the excited continuum states leave the compact region ρ. In the limit of large t, the only part of the wave function that remains upon projection by P is from the bound state |ψ 0 . This statement is true for any Hamiltonian with a translationally-invariant kinetic energy and exactly one localized bound state. If P |ψ I has well-defined quantum numbers associated with an exact symmetry of H, then the stable fixed point property applies to cases where there is exactly one localized bound state in the symmetry sector of interest.
For any two states |x and |y , we define the normalized overlap to be | x|y |/ x|x y|y . Let O(t) be the normalized overlap between P |ψ 0 and P U (t)P |ψ I . In Fig. 1 we show O(t) for five randomly chosen initial states |ψ I versus time t for Model 1A. In this example we take R = 5 and take L large enough to prevent reflection waves returning from the boundary. We see that in all cases the normalized overlap function approaches 1, demonstrating that P |ψ 0 is a stable fixed point.
In order to apply projected cooling to a Hamiltonian with more than one bound state, we will need to use a time dependent Hamiltonian H(t). In the first stage of the time evolution, we let H(t) be a Hamiltonian with only one bound state, but one that has good overlap with the ground state of H. One simple and effective way to achieve this is to multiply the kinetic energy operator by a scale factor greater than 1 until all the bound excited states are pushed into the continuum. This augmentation of the kinetic energy operator is also very helpful to accelerate the time evolution at early times, as measured by the number of quantum gate operations per qubit. Once the time-evolved state |ψ(t) is closely tracking the ground state of H(t), we can then use adiabatic evolution to gradually reach the desired ground state of H.
As an example we consider Model 1B, which is defined in exactly the same manner as Model 1A except that we take V n = −1.6δ 0,n − 1.5(δ 2,n + δ 3,n ) − 1.4δ −2,n . This change is enough to produce four bound states. We first compute results using the full time evolution operator, U (t, t − ∆t) = e −iH(t)∆t . We then also use the Trotter approximation [11] and break apart the Hamiltonian into pieces where A n ′ ,n is the off-diagonal part of K n ′ ,n when min(n ′ , n) is even, B n ′ ,n is the off-diagonal part of K n ′ ,n when min(n ′ , n) is odd, and D n ′ ,n is the diagonal part of K n ′ ,n . The Trotterized time evolution operator is then In addition to the Trotter approximation, we also consider the effect of stochastic noise upon the time evolution. After each ∆t time step, we multiply each component of the evolved wave function by a factor 1 + z, where z is a complex Gaussian random variable centered at zero with root-mean-square values for the real and imaginary parts both equal to ε/ √ 2.
For the adiabatic evolution calculations, we start with initial Hamiltonian H I = V . For the initial state we use the ground state of H I , which is the point-like wave function |ψ 1 . For the projected cooling calculations, we can use any initial state contained within the region ρ. In addition to the point-like wave function |ψ 1 I , we also use the smeared wave function In Panel A of Fig. 2 we show the normalized overlap O(t) between the evolved wave function and the exact wave function over the interior region ρ versus the number of time steps N t for Model 1B. We take R = 5, L = 25, and the time step ∆t is 0.3. This corresponds to an interior region ρ with 11 dimensions in the one-particle sector. AE corresponds to adiabatic evolution, while PC corresponds to projected cooling. Full evolution denotes evolution using the full time-dependent Hamiltonian for each time step, while Trotter evolution denotes the Trotter approximation. Point initial indicates the initial state |ψ 1 I , while spread initial indicates the initial state |ψ 2 I . The quoted numerical error corresponds to the value of the parameter ε. More details of these calculations can be found in the Supplemental Materials. We see that standard adiabatic evolution has difficulties finding the ground state, achieving an overlap of no more than 0.35. In contrast, the projected cooling algorithm is able to achieve an overlap of at least 0.94 in 40 time steps or less for all cases, even with the errors due to Trotter approximation and stochastic noise of size ε = 0.05. For our next model, Model 2, we consider a Hamiltonian defined on two linked one-dimensional chains of 2L + 1 qubits each, n 1 = −L, · · · L and n 2 = −L, · · · L. We again define the vacuum to be the tensor product state where all qubits are |0 . This time we consider the two-particle sector and define |[n 1 , n 2 ] as the tensor product state where qubit n 1 on the first chain is |1 , qubit n 2 on the second chain is |1 , and all other qubits are |0 . For the interior region ρ, we take the sites where max(|n 1 |, |n 2 |) ≤ R. We again will use the values R = 5, L = 25, and ∆t = 0.3. This corresponds to an interior region with 121 dimensions where the kinetic energy term K n ′ ,n is the same as in Model 1A and Model 1B. For the interactions we take V n = −1.0δ 0,n + 0.2δ 1,n − 0.9(δ 2,n + δ 3,n ) − 0.3δ −1,n for the single-particle potential energy and W n 1 ,n 2 = −0.2δ n 1 ,n 2 for the two-particle interaction. For this model there are again four localized bound states that remain localized in region ρ. We use the point-like initial state |ψ 1 I = |[0, 0] for adiabatic evolution. For projected cooling we consider |ψ 1 I as well as the smeared initial state As we did for Model 1B, we simulate noise after each time step by multiplying each component of the evolved wave function by a factor 1 + z, where z is a Gaussian random complex variable whose size is parameterized by ε. We refer the reader to the Supplemental Materials for details of the adiabatic evolution times steps is shown on the right in Panel C. For these plots we used full time evolution with error ε = 0 and the spread initial state |ψ 2 I for the projected cooling calculation. We see that the improvement obtained using projected cooling is quite significant.
The results presented here are typical of the performance one can obtain using projected cooling. Because of the fixed point properties of projected cooling when H(t) has only one bound state, the method is flexible, efficient, and resilient against small errors. The projected cooling algorithm is able to construct the localized ground state of any Hamiltonian with a translationally-invariant kinetic energy, and the only additional resource required is using a volume that is significantly larger than the size of the ground state. While more work needs to be done to optimize the algorithm, it can now be applied to quantum computations of the structure of atomic nuclei as well as other bound quantum systems [12].

Implementation of the model Hamiltonians
We can implement Models 1A and 1B using a single one-dimensional chain of qubits. The underlying qubit Hamiltonian has the form where Similarly, we can implement Model 2 using two linked one-dimensional chains of qubits. The terms in the Hamiltonian have the form where

Time evolution of Model 1B
For Model 1B, our one-particle Hamiltonian has the form H = K + V with where the kinetic energy term is K n ′ ,n = δ n ′ ,n − 1 2 δ n ′ ,n+1 − 1 2 δ n ′ ,n−1 , and the interaction is V n = −1.6δ 0,n − 1.5(δ 2,n + δ 3,n ) − 1.4δ −2,n . For the Trotter approximation we break apart the Hamiltonian into pieces n is the off-diagonal part of K n ′ ,n when min(n ′ , n) is even, B n ′ ,n is the off-diagonal part of K n ′ ,n when min(n ′ , n) is odd, and D n ′ ,n is the diagonal part of K n ′ ,n . For notational clarity when discussing time dependent operators, we writeH,K,V , etc., to denote these static operators. The Hamiltonian has the form H = K + V + W with [n ′ 1 , n ′ 2 ]| H |[n 1 , n 2 ] = K n ′ 1 ,n 1 δ n ′ 2 ,n 2 + K n ′ 2 ,n 2 δ n ′ 1 ,n 1 + (V n 1 + V n 2 + W n 1 ,n 2 )δ n ′ 1 ,n 1 δ n ′ 2 ,n 2 , where the kinetic energy term is K n ′ ,n = δ n ′ ,n − 1 2 δ n ′ ,n+1 − 1 2 δ n ′ ,n−1 . For the interactions we take V n = −1.0δ 0,n + 0.2δ 1,n − 0.9(δ 2,n + δ 3,n ) − 0.3δ −1,n and W n 1 ,n 2 = −0.2δ n 1 ,n 2 . For the Trotter approximation we break apart the Hamiltonian into pieces where A [i] n ′ i ,n i is the off-diagonal part of the kinetic energy for particle i when min(n ′ i , n i ) is even, B [i] n ′ i ,n i is the off-diagonal part of the kinetic energy when min(n ′ i , n i ) is odd, and D is the diagonal part of the kinetic energy. When discussing time dependent operators, we writeH,K,V , etc., to denote these static operators.
For the adiabatic evolution from initial time t = 0 to final time t = t F , we use the time-dependent Hamiltonian H(t) = t t FK +V +W .