Quantum algorithm for spectral projection by measuring an ancilla iteratively

TL;DR

This paper introduces a quantum algorithm that projects a system onto eigenstates of any Hermitian operator using iterative ancilla measurements, useful for quantum annealing and simulating many-body systems.

Contribution

It presents a novel iterative spectral projection algorithm utilizing ancilla measurements, applicable to quantum annealing and many-body quantum simulations.

Findings

Algorithm successfully projects onto eigenstates following the Born rule.

Simulation demonstrates effectiveness in driving spin chains to ground states.

Potential integration as a subroutine in quantum optimization procedures.

Abstract

We propose a quantum algorithm for projecting a quantum system to eigenstates of any Hermitian operator, provided one can access the associated control-unitary evolution for the ancilla and the system, as well as the measurement of the controlling ancillary qubit. Such a Hadamard-test like primitive is iterated so as to achieve the spectral projection, and the distribution of the projected eigenstates obeys the Born rule. This algorithm can be used as a subroutine in the quantum annealing procedure by measurement to drive the system to the ground state of a final Hamiltonian, and we simulate this for quantum many-body spin chains.

Tables

Table 1. Table 1 : Comparison of the standard QPE  NielsenChuang , the iterative QPE (iQPE)  Dobsicek2007 and our spectral projection algorithm (SPA). QFT stands for quantum Fourier transform. Both the QPE and iQPE have fixed accuracy set by the choice of highest power in U 2 k superscript 𝑈 superscript 2 𝑘 U^{2^{k}} and during the procedure the highest power cannot be changed; the QPE is fixed by the total number of ancillas and the iQPE needs to fix the highest power in the beginning of the procedure. Both the QPE and iQPE require precise execution of c − U 2 k 𝑐 superscript 𝑈 superscript 2 𝑘 c-U^{2^{k}} for all k < t g 𝑘 subscript 𝑡 𝑔 k<t_{g} . On the other hand, our SPA uses c − e − i ​ h ^ ​ Δ ​ t 𝑐 superscript 𝑒 𝑖 ^ ℎ Δ 𝑡 c-e^{-i\hat{h}\Delta t} and the range of Δ ​ t Δ 𝑡 \Delta t can be fixed, but the accuracy can still be improved by running more iterations. Our SPA does not require Δ ​ t Δ 𝑡 \Delta t to be exact 2 k superscript 2 𝑘 2^{k} , and in fact it can be somewhat arbitrary. The drawback of our SPA is that the number of required iterations for achieving a fixed accuracy can vary from run to run.

Methods	Projection capability	Phase estimation capability	Accuracy limitation
QPE	yes	yes	no. of ancillas & power in $c-U^{2^k}$; requires QFT
iQPE	yes	yes	power in $c-U^{2^k}$; requires no QFT
SPA	yes	yes	no. of iterations; requires no QFT

Equations

$$|{\psi }_m^{\prime }⟩$$

$$=\frac{1}{\sqrt{2}}[\alpha +(-1{)}^m{\mathrm{e}}^{-i\xi }\beta e^{-i\hat{t}\Delta t}]|\psi ⟩,$$

$$\displaystyle p_{m}=||\psi^{\prime}_{m}||^{2}=\frac{1}{2}\big{[}1+2(-1)^{m}{\rm Re}(\alpha^{*}\beta e^{-i\xi}\langle\psi|e^{-i\hat{h}\Delta t}|\psi\rangle)\big{]}.$$

$$p_{j,m}$$

$$=$$

$$\Delta E_{(m)}=\frac{2(-1)^{m}\big{(}{\cal R}_{h}-\langle h\rangle{\cal R}_{1}\big{)}}{1+2(-1)^{m}{\cal R}_{1}}$$

$${\mathcal{R}}_1$$

$${\mathcal{R}}_h$$

$$\Delta E_{(m)}=\frac{2(-1{)}^m\mathrm{Im}({\alpha }^{\ast }\beta )}{1+2(-1{)}^m\mathrm{Re}({\alpha }^{\ast }\beta )}⟨\Delta h^2⟩\Delta t,$$

$$\Delta E_{(m)}=-\frac{{\rm Re}(\alpha^{*}\beta)\big{(}\langle\hat{h}^{3}\rangle-\langle\hat{h}^{2}\rangle\langle\hat{h}\rangle\big{)}(\Delta t)^{2}}{1+2(-1)^{m}{\rm Re}(\alpha^{*}\beta)}.$$

$$\Delta E=\left\{\begin{array}[]{ll}\frac{1}{4}(\langle\hat{h}^{3}\rangle-\langle\hat{h}^{2}\rangle\langle\hat{h}\rangle)(\Delta t)^{2},&m=0,\\ \phantom{}&\phantom{}\\ \frac{\langle\hat{h}^{3}\rangle}{\langle\hat{h}^{2}\rangle}-\langle\hat{h}\rangle,&m=1,\end{array}\right.$$

$$\overline{\delta V_E}$$

$$\overline{\delta V_{E}}=\langle\psi|\hat{h}|\psi\rangle^{2}-\big{[}p_{0}\langle\tilde{\psi}^{\prime}_{0}|\hat{h}|\tilde{\psi}^{\prime}_{0}\rangle^{2}+p_{1}\langle\tilde{\psi}^{\prime}_{1}|\hat{h}|\tilde{\psi}^{\prime}_{1}\rangle^{2}\big{]}.$$

$$\overline{\delta V_{E}}=\frac{-4}{1-4{\cal R}_{1}^{2}}\big{(}{\cal R}_{1}\langle h\rangle-{\cal R}_{h}\big{)}^{2},$$

$$\displaystyle\overline{{\delta V_{E}}}=-\frac{4\,{\rm Im}(\alpha^{*}\beta)^{2}}{1-4{\rm Re}(\alpha^{*}\beta)^{2}}\big{(}\langle\psi|(\Delta h)^{2}|\psi\rangle\big{)}^{2}(\Delta t)^{2}.$$

$$\overline{{\delta V_{E}}}=-\frac{{\rm Re}(\alpha^{*}\beta)^{2}\big{(}\langle\hat{h}^{3}\rangle-\langle\hat{h}^{2}\rangle\langle\hat{h}\rangle\big{)}^{2}(\Delta t)^{4}}{1-4{\rm Re}(\alpha^{*}\beta)^{2}},$$

$$\displaystyle\overline{{\delta V_{E}}}=-\frac{\big{(}\langle\hat{h}^{3}\rangle-\langle\hat{h}^{2}\rangle\langle\hat{h}\rangle\big{)}^{2}}{4\langle\hat{h}^{2}\rangle}(\Delta t)^{2}.$$

$$|{\psi }_m^{\prime }⟩\approx \frac{\alpha +(-1{)}^m\beta }{\sqrt{2}}\left[1+\frac{-i\hat{h}\Delta t-\frac{1}{2}(\hat{h}\Delta t{)}^2}{1+(-1{)}^m\alpha /\beta }\right]|\psi ⟩.$$

$${\hat{P}}_m=\frac{-i\hat{h}\Delta t-\frac{1}{2}(\hat{h}\Delta t{)}^2}{1+(-1{)}^m\alpha /\beta }+\frac{1}{2}\frac{(\hat{h}\Delta t{)}^2}{[1+(-1{)}^m\alpha /\beta {]}^2}.$$

$$e^{{\hat{Q}}_1(\alpha ,\beta ,\Delta t)}{e}^{{\hat{Q}}_0(\alpha ,\beta ,\Delta t)}{e}^{{\hat{Q}}_1(\alpha ,\beta ,\Delta t)}\cdots e^{{\hat{Q}}_0(\alpha ,\beta ,\Delta t)},$$

$$\sum _{m=0,1}{p}_m{\hat{Q}}_m=-\frac{\mathrm{Im}({\alpha }^{\ast }\beta {)}^2\Delta t^2}{1-4\mathrm{Re}({\alpha }^{\ast }\beta {)}^2}[(\hat{h}-⟨\hat{h}⟩{)}^2-⟨\hat{h}{⟩}^2].$$

$$|\psi ⟩\to e^{-c(\alpha ,\beta )\Delta t^2(\hat{h}-⟨\hat{h}⟩{)}^2}|\psi ⟩.$$

$$|{\psi }_m^{(k)}⟩=\frac{1}{\sqrt{2}}\left[{\alpha }_k|{\psi }^{(k-1)}⟩+(-1{)}^m{\beta }_k{U}_k|{\psi }^{(k-1)}⟩\right],$$

$$H_5=\left(\begin{array}{ccccc}-0.0763231& \ast & \ast & \ast & \ast \\ -0.51328+0.0732759i& 0.691614& \ast & \ast & \ast \\ 0.516039+0.20004i& -0.884252-0.248885i& -0.495554& \ast & \ast \\ -0.379429+0.303255i& 0.0981619-0.603679i& -0.484382-0.134895i& 0.921927& \ast \\ 0.0142526+0.421276i& 0.635987+0.0817911i& -0.450215-0.808964i& 0.6387+0.188711i& 0.736562\end{array}\right),$$

$$|{\psi }_0⟩=(0.506424,-0.370456+0.164849i,-0.444258+0.194814i,-0.0372888-0.33439i,-0.475495-0.0671035i{)}^T.$$

$$\{0.554875,0.0729256,0.262368,0.00841186,0.10142\}.$$

$$H_{\mathrm{TFI}}(g)=\sum _{i=1}^{N_q}\left[g{\sigma }_i^x{\sigma }_{i+1}^x-(1-g){\sigma }_i^z\right].$$

$$|\psi ⟩=\sum _{i;E_i=-E}{a}_i|E_i⟩+\sum _{j;E_j=E}{b}_j|E_j⟩,$$

$$⟨\hat{h}⟩$$

$$⟨{\hat{h}}^2⟩$$

$$⟨{\hat{h}}^3⟩$$

$$|\Psi (T)⟩=\hat{T}{e}^{-i{\int }_0^{T}H(g(t))dt}|G(0)⟩,$$

$$H_{\mathrm{XzY}}(g,r)$$

Full text

I Introduction

The measurement postulate of quantum mechanics states that when measuring an observable $\hat{o}$, only its eigenvalues $o_{n}$ will be observed and the state of the system will be projected to the corresponding eigenstate $|o_{n}\rangle$, for which $\hat{o}|o_{n}\rangle=o_{n}|o_{n}\rangle$, immediately after the measurement. Furthermore, the Born rule prescribes the probability of such an outcome for an initial quantum state $|\psi_{0}\rangle$ as $p_{n}=|\langle o_{n}|\psi_{0}\rangle|^{2}$. Whether one can derive the rule and hence remove it from the postulates of quantum mechanics is still of fundamental interest Masanes2018 . From the perspective of quantum information processing, general construction of such spectral projection is also of practical importance. For example, Ref. Poulin2018 constructs a quantum walk approach to achieve this and emphasizes its utility in carrying out a key step of the quantum simulated annealing (QSA) algorithm for optimization problems Somma2008 . The latter can be used as an alternative to the adiabatic quantum computation (AQC) Adia1 ; Adia2 . In fact, the standard quantum phase estimation (QPE) NielsenChuang and its variants Kitaev2002 ; Aspuru-Guzik2005 ; Dobsicek2007 can also achieve approximate spectral projection when the system is not in an eigenstate.

The QPE is crucial in many quantum information processing applications NielsenChuang , including factoring and, more relevant to the present paper, the quantum-walk spectral measurement in Ref. Poulin2018 , as well as related methods for preparing a thermal Gibbs state Poulin2009 ; Temme2012 ; Yung2012 ; Moussa2019 . The standard QPE uses ${\cal O}(t_{g})$ controlled unitary gates of the form $c-U^{2^{k}}$ (for $k=0$ to $t_{g}-1$) to encode the $t_{g}$ binary digits of the phase value (in unit of $2\pi$) and it requires ${\cal O}(t_{g}^{2})$ gates in the inverse quantum Fourier transform to retrieve the phase NielsenChuang . Regarding the accuracy of QPE, in order to have the phase accurate in $m$ binary digits with the success probability of at least $1-\epsilon$, the total number of ancillary qubits needed is $t_{g}=m+\log(2\epsilon+1/2\epsilon)$ NielsenChuang . In other words, using $t_{g}$ ancillary qubits allows the phase value to be accurate in $t_{g}-\log(2\epsilon+1/2\epsilon)$ binary digits. The accuracy in the phase is thus limited by the number of available ancillas employed in representing the value of the phase, and when used as spectral projection subroutine, the eigenstate the system is projected by the QPE to is only approximate. The unitary $U$ may be implemented by $e^{-i\hat{o}\Delta t}$, and in the QPE, the power in the unitary $U$ needs to go as large as $2^{t_{g}-1}$; equivalently, the timing $\Delta t$ needs to be made accurate to $2^{k}$ (for $k=0$ to $t_{g}-1$). Maintaining the stability of $U$ and coherence of the quantum register when carrying out the QPE is important for noisy intermediate-scale quantum processors.

Here, we apply a simple iterative approach to achieve the spectral projection of an associated observable $\hat{o}$, and in each step of the iteration only one ancilla is used as the control to enact a unitary evolution (${\rm c-}e^{-i\Delta t\hat{o}}$) on the system, conditioned on the ancillary state being $|1\rangle$. Then only the ancilla is measured in the Pauli X basis. After sufficient number of steps have been carried out (see below), the system is projected to an eigenstate of the operator $\hat{o}$. We demonstrate by numerical simulations that our procedure can lead to spectral projection by varying the parameter $\Delta t$ and the ancilla’s state parameter.

To understand that repeated application of the primitive eventually leads to spectral projection, we provide two perspectives. First, we show that on average the energy variance of the system will decrease; see Eq. (13). If the energy variance decreases to zero, then an eigenstate is reached. Second, an intuitive picture of our procedure emerges: at each step, the measurement of the ancillary qubit gives rise to a random walk in the operator action, i.e. with either $e^{\hat{Q}_{0}}$ or $e^{\hat{Q}_{1}}$ acting on the system. The choice of which operators depends on the measurement outcome; see Fig. 1 below. The key notable difference from the conventional random walk is that the outcome probability is state dependent. However, we calculate the average random-walk action $p_{0}\hat{Q}_{0}+p_{1}\hat{Q}_{1}$ per step that is valid in the small $\Delta t$ limit, and find that it leads to a map, see Eq. (19), that when repeated will drive the system to an eigenstate. Both viewpoints validate that our procedure can lead to spectral projection, as eigenstates have no energy (or observable-value) variance and are fixed points of the iterative procedure.

We emphasize that the time $\Delta t$ here, unlike in the QPE, does not need to be exactly of the form $2^{k}$. Thus, in some sense the protocol for spectral projection does not require exact timing and can tolerate fluctuations and imprecision in timing. In addition, the range of $\Delta t$ used needs not span over many orders of magnitudes related to the accuracy of the QPE, i.e., $\max\{\Delta t\}/\min\{\Delta t\}$ can be much smaller than $2^{t_{g}-1}$. Moreover, the ancilla state does not need to be in the $|\pm\rangle$ state right before the controlled unitary and it can be in almost any pure state. As seen below, we can also used a fixed $\Delta t$ in our procedure to achieve the spectral projection.

Given that spectral projection can be achieved, one immediate question is what governs the distribution of the projected eigenstates. For this we show that the distribution of this eigenstate projection obeys the Born rule. Fundamentally, our algorithm can be regarded as a procedure to achieve the effect described in the measurement postulate. As an application, we simulate the use of our spectral projection algorithm in two spin-chain models, and demonstrate that ground states at different transverse field strengths can be successfully obtained, when there is a gap in the Hamiltonian throughout the parameter range of interest.

Our initial motivation for this study comes from the incentive to devise a simple quantum version of Lanczos algorithm. An approach was recently proposed in Ref. Motta2019 by implementing an effective unitary evolution $e^{-ih_{\rm eff}\Delta\tau}$ to simulate the effect of imaginary time evolution $e^{-h\Delta\tau}$ on a quantum state. We wish to develop an alternative approach that does not require the searching of the effective Hamiltonian $h_{\rm eff}$. However, we could not make the procedure to work due to high-order effect, and we describe such a failed attempt in the Appendix. However, it was by analyzing this that leads us to the spectral projection algorithm and the understanding why the attempt failed.

The remainder of the paper is organized as follows. In Sec. II we discuss a primitive that slightly generalizes the Hadamard test by using a general ancillary state. By repeating this primitive with sufficient number of times, we argue that it will project the system to an eigenstate. In Sec. III we describe the approach to classically simulate the above procedure and verify by simulations that it indeed leads to an eigenstate or spectral projection algorithm. There, we use random Hermitian matrices for illustration and also demonstrate that such spectral projection obeys the Born rule for the final distribution of projected eigenstates. In Sec. IV we give illustrations of our spectral projection algorithm using the quantum transverse-field Ising spin chain. In Sec. V we discuss the effect of decoherence. In Sec. VI we illustrate the use of our spectral projection algorithm in the quantum annealing for two different spin chains. Finally, in Sec. VII, we make some concluding remarks.

II The primitive and the alogrithm for spectral projection

The primitive that our algorithm is based on is similar to the Hadamard test and will be described below. The algorithm itself is a repeated application of such a primitive. We will provide analysis to support that our algorithm can achieve spectral projection.

II.1 The Hadamard test and the primitive

The basic idea of our approach is to entangle a system with an ancilla qubit prepared in a certain state, and then measure the ancilla in a chosen basis, similar to the so-called Hadamard test. This is commonly used in many quantum information processing protocols NielsenChuang . We will describe a slightly varied primitive, in which the ancilla needs not be in the $|+\rangle$.

Let the system be in an initial state $|\psi\rangle$ and an ancilla in $|A\rangle=\alpha|0\rangle+\beta|1\rangle$ with $|\alpha|^{2}+|\beta|^{2}=1$. We entangle the ancilla (as the control) and the system (as the target) by the controlled operation ${c-U}\equiv|0\rangle\langle 0|\otimes\mathds{1}+|1\rangle\langle 1|\otimes U$, where $U=\exp{(-i\hat{h}\Delta t)}$ is the unitary evolution under a Hamiltonian $\hat{h}$ within a duration $\Delta t$. We then measure the controlling ancillary qubit in the basis $(|0\rangle\pm\mathrm{e}^{i\xi}|1\rangle)/\sqrt{2}$, with the $\pm$ associated with the measurement outcome $m=0$ or $1$, respectively. This is equivalent to measuring the observable $\cos(\xi)\sigma_{x}+\sin(\xi)\sigma_{y}$ on the ancilla. We shall see below that we can take $\xi=0$ without loss of generality, and thus the measurement will correspond to the Pauli X basis, and the primitive is illustrated in Fig. 1.

The measurement of the ancilla then collapses the system to the unnormalized state:

[TABLE]

and the corresponding probability of obtaining the outcome $m$ is

[TABLE]

Here we see that the phase factor $e^{-i\xi}$ from the measurement basis can be absorbed into the ancilla’s initial state parameter $\beta$, and thus we can set $\xi=0$ from now on without loss of generality, resulting in the ancilla measurement to the fixed the Pauli X, whose eigenstates are simply $|\pm\rangle\equiv(|0\rangle\pm|1\rangle)/\sqrt{2}$.

Eigenstates are fixed points of the primitive. It is easy to see that for any eigenstate $|E_{j}\rangle$ with eigenenergy $E_{j}$, the post-measurement state is still $|E_{j}\rangle$, but the probability of getting the $m$-th outcome is

[TABLE]

The probabilities for ‘0’ and ‘1’ outcomes add up to unity: $p_{j,0}+p_{j,1}=1$. Moreover, their difference $p_{j,0}-p_{j,1}=2{\rm Re}(\alpha^{*}\beta e^{-iE_{j}\Delta t})=2|\alpha\beta|\cos(\phi-E_{j}\Delta t)$ can be used to determine $E_{j}\Delta t$ up to an overall sign and multiples of $2\pi$, where $\alpha^{*}\beta=|\alpha\beta|e^{i\phi}$. To uniquely determine $E_{j}$, one can use a different set of $(\alpha,\beta)$ and $\Delta t$ to obtain different distributions for estimation. Note that in order to achieve optimal determination we can maximize $|\alpha\beta|$, which is achieved when $|\alpha|=|\beta|=1/\sqrt{2}$ and corresponds to using an ancillary state $|A\rangle=(|0\rangle+e^{i\phi}|1\rangle)/\sqrt{2}$. The choice of $|+\rangle$ is the typical ancillary state in the Hadamard test.

Suppose we have two different energy eigenstates with distinct energies $E_{j}\neq E_{k}$, generically the two distributions are different, $p_{j,m}\neq p_{k,m}$, unless the choice of $(\alpha,\beta)$ and $\Delta t$ coincidentally make ${\rm Re}(\alpha^{*}\beta e^{-iE_{j}\Delta t})={\rm Re}(\alpha^{*}\beta e^{-iE_{k}\Delta t})$. Hence, by accumulating enough statistics, one can determine whether the two eigenstates have the same energy or not. One can use different ‘distance’ measures, such as the relative entropy to quantify the distinguishability. Quantities such as the Chernoff bound can also be used to quantify the likelihood of deviating from the average values and thus the degree of distinguishability for a finite number of measurements performed.

For the system’s initial state being $|\psi\rangle=\sum_{j}c_{j}^{(0)}|E_{j}\rangle$, then the probability of getting outcome $m$ in the ancilla’s measurement can be shown to be $p_{m}^{(0)}=\sum_{j}|c_{j}^{(0)}|^{2}p_{j,m}$, which is a convex mixture of the extremal distributions $p_{j,m}$ from the eigenstates. This means that the knowledge of the distribution ${p}_{m}^{(0)}$ is not sufficient to infer uniquely the compositions of the eigenstates. If we are given only a copy of $|\psi\rangle$ and if it is not in an eigenstate, then it is not possible to estimate $p_{m}^{(0)}$ by any measurement.

Energy change.

First, we can ask how much the energy has changed after one such a primitive step: $\Delta E_{(m)}\equiv\langle\tilde{\psi^{\prime}}_{m}|\hat{h}|\tilde{\psi^{\prime}}_{m}\rangle-\langle\psi|\hat{h}|\psi\rangle$, where $|\tilde{\psi^{\prime}}_{m}\rangle\equiv|\psi^{\prime}_{m}\rangle/\sqrt{p_{m}}$ is the normalized post-measurement state. By using the expressions for the post-measurement state $|\psi^{\prime}_{m}\rangle$ (1) with $\xi=0$ and the probability $p_{m}$ (2), we can calculate $\Delta E_{(m)}$ explicitly and arrive at (see Appendix B for derivations)

[TABLE]

where parameters ${\cal R}_{1}$ and ${\cal R}_{h}$ are defined as

[TABLE]

In order to obtain some intuition of the above expression, we can expand it to the first nonvanishing order. We find that for ${\rm Im}(\alpha^{*}\beta)\neq 0$, the lowest nonvanishing contribution occurs at the first order in $\Delta t$,

[TABLE]

where the expectation $\langle\cdots\rangle$ is evaluated w.r.t. $|\psi\rangle$, e.g. $\langle\hat{h}\rangle\equiv\langle\psi|\hat{h}|\psi\rangle$, and the $\langle(\Delta h)^{2}\rangle\equiv\langle\psi|\hat{h}^{2}|\psi\rangle-\langle\psi|\hat{h}|\psi\rangle^{2}$ is the energy variance of the state $|\psi\rangle$. Thus, generically the change in the energy after one step is proportional to the energy variance before the application of the primitive.

We note that, however, when ${\rm Im}(\alpha^{*}\beta)=0$, the change in energy is in the second order,

[TABLE]

Of course, the exception is when ${\rm Re}(\alpha^{*}\beta)=\pm 1/2$, which corresponds to the case of the ancillary state being $|\pm\rangle$. In the case of using $|+\rangle$ of the ancilla, the change in energy to the first nonvanishing contribution is

[TABLE]

The above two outcomes are switched, if the ancilla’s state is $|-\rangle$.

It is interesting to observe that in general the energy will always change if $|\psi\rangle$ is not an eigenstate, except when the ancilla’s state satisifies ${\rm Im}(\alpha^{*}\beta)=0$ (e.g. in the $|\pm\rangle$ state) and the system satisfies $\langle\hat{h}^{3}\rangle-\langle\hat{h}^{2}\rangle\langle\hat{h}\rangle=0$, then the energy will not change.

II.2 The algorithm

The algorithm is simply a procedure that repeats the above primitive many times. In each step, the ancilla parameters $(\alpha,\beta)$ and the duration $\Delta t$ can be different. In the following, we provide argument to support that our algorithm indeed will achieve spectral projection, by analyzing two quantities that characterize the average effect, in terms of the energy variance and the average action of a random walk.

Before we proceed to the analysis, we need to first ask the question: how do we know the algorithm has produced a converged eigenstate? Assume that the system converges to an eigenstate $|E_{j}\rangle$. Applying the $c-U$ gate to the ancilla (initially in $\alpha|0\rangle+\beta|1\rangle$) and the system leaves the system intact but changes the relative phase in the ancilla: $\alpha|0\rangle+\beta e^{-iE_{j}\Delta t}|1\rangle$. Following the idea in the so-called eigenstate witness method Santagati2018 , one can perform quantum state tomography on the identically prepared ancillary qubit after applying the control unitary. A single-qubit tomography involves measurement in Pauli X, Y, and Z bases, as for a general one-qubit mixed state $\rho_{a}=(I+\sum_{i=x,y,z}r_{i}\sigma_{i})/2$, its parameters can be obtained from measurement: $r_{i}={\rm Tr}(\rho_{a}\sigma_{i})$. In our algorithm, we perform X measurement, but measurement in Y can also achieve spectral projection, as we have argued that the measurement phase $\xi$ in $|0\rangle\pm e^{i\xi}|1\rangle$ was conveniently aborbed in the ancilla parameter $\beta$. We can also perform Pauli $Z$ measurement, which will project the ancilla and the system to $|0\rangle|\psi\rangle$ or $|1\rangle e^{-i\hat{h}\Delta t}|\psi\rangle$, and it does not affect the spectral projection. If the state tomography shows that it remains in a pure state, then the system must be in an eigenstate and hence it is converged. One can use the purity of the resultant ancillary state as a measure for convergence. Moreover, from the tomography, the quantity $E_{j}\Delta t$ can also be determined (up to multiple of $2\pi$). By using different sets of $\Delta t$, one can then uniquely determine $E_{j}$.

Perspectives from probability distribution. As we shall see below, successive coupling of individual ancillas with the system and measurement on ancillas help to drive the system to an eigenstate, therefore arriving at the final distribution $p_{m}=p_{j^{*},m}$ for some $j^{*}$ labeling the eigenstate $|E_{j^{*}}\rangle$. The procedure gives rise to a sequence of 0/1 outcomes, namely a bit string $0,1,1,0,\dots$ from the ancillas’ measurement. From the persective of probability distribution, it starts with $p_{m}^{(0)}$, and after successive application of the primitive, the distribution flows: $p_{m}^{(0)}\rightarrow p_{m}^{(1)}\rightarrow\cdots p_{m}^{(n)}$, and after $n$ steps, $p_{m}^{(n)}$ will be close to some $p_{j^{*},m}$. (In terms of coins, there are $n+1$ different coins.) We stress, however, that in the process we cannot obtain the distributions $p_{m}^{(k)}$ from measurement as we are given only one copy of the system, but we only obtain a sequence of 0 and 1. From this sequence we can only estimate the average probability of getting 0 and 1, i.e. $\overline{p_{m}}$. In the case of the eigenstates, such average distribution can be used to distinguish whether two given eigenstates have different energy or not (as the eigenstate does not change, and hence the ‘coins’ are identical, as we have discussed previously). However, for an arbitrary initial state of the system, does the knowledge of $\overline{p_{m}}$ guarantee the projection?

Given that our algorithmic procedure enables the projection to eigenstates (as argued and numerically demonstrated below), then from the bit string and the knowledge of the initial state of the system, one can indeed infer the distributions $p^{(k)}_{m}$, as well as whether and what energy eigenstate is arrived and what the corresponding eigenenergy is by direct classical simulations. However, for large system sizes, classical simulations will not be possible. How do we argue that our procedure indeed leads to spectral projection? How do we explain that in the limit of long seqence $\overline{p_{m}}$ will eventually flow to an extremal or fixed-point distribution $p_{j^{*},m}$? In the following, we provide two physically motivated approaches to understand the spectral projection.

Energy variance. If the energy variance of a quantum state is zero, the state is an energy eigenstate, i.e., $V_{E}(\psi)=0\leftrightarrow\hat{h}|\psi\rangle=E|\psi\rangle$, where $V_{E}(\psi)\equiv\langle\psi|(\Delta h)^{2}|\psi\rangle=\langle\psi|\hat{h}^{2}|\psi\rangle-\langle\psi|\hat{h}|\psi\rangle^{2}$. Thus, energy variance is an important indicator to how close the state has converged to an eigenstate. Carrying out the primitive yields the outcome ‘0’ with probability $p_{0}$ and the normalized post-measurement system state $|\tilde{\psi}^{\prime}_{0}\rangle\equiv|\psi^{\prime}_{0}\rangle/\sqrt{p_{0}}$, and the outcome ‘1’ with probability $p_{1}$ and the normalized post-measurement system state $|\tilde{\psi}^{\prime}_{1}\rangle/\equiv|\psi^{\prime}_{1}\rangle/\sqrt{p_{1}}$. Given the probabilistic nature due to measurement, it is thus natural consider the average change of the energy variance after one step:

[TABLE]

Since the expectation value of $\hat{h}$ and its function such as $\hat{h}^{2}$ are conserved, the above change can be simplified to be

[TABLE]

By using the expressions for the post-measurement state (1) and the probability (2), we can calculate $\overline{\delta V_{E}}$ explicitly and arrive at (see Appendix B for derivations)

[TABLE]

where the two parameters ${\cal R}_{1}$ and ${\cal R}_{h}$ are defined previously as in Eqs. (6) and (7), respectively. Given that $|\alpha^{*}\beta|\leq 1/2$, the parameter $R_{1}$ satisfies $R_{1}^{2}\leq 1/4$. Thus, $\overline{\delta V_{E}}\leq 0$ and generically $\overline{\delta V_{E}}<0$.

In order to obtain some intuition about $\overline{\delta V_{E}}$, we can expand it in series of $\Delta t$, and we find that when ${\rm Im}(\alpha\beta)\neq 0$, it is nonvanishing at the second order:

[TABLE]

The factor, c(\alpha,\beta)\equiv{{\rm Im}(\alpha^{*}\beta)^{2}}/\big{[}{1-4{\rm Re}(\alpha^{*}\beta)^{2}}\big{]} is maximized with a value $1/4$ when $\beta/\alpha=e^{i\phi}$ and $|\alpha|=|\beta|=1/\sqrt{2}$. Such a choice represents a maximum ‘rate’ of change in the energy variance. Moreover, the change is also proportional to the square of the energy variance of the system’s state before carrying out the primitive.

We note that, however, when ${\rm Im}(\alpha^{*}\beta)=0$,

[TABLE]

except when $\alpha=1/\sqrt{2}$ and $\beta=\pm 1/\sqrt{2}$, the average change in the energy variance is

[TABLE]

We observe that the quantity $\langle\hat{h}^{3}\rangle-\langle\hat{h}^{2}\rangle\langle\hat{h}\rangle$ has previously appeared in the change of the energy (10).

The above analysis suggests that we should usually choose ancillary parameters such that ${\rm Im}(\alpha^{*}\beta)\neq 0$, except when $\alpha=1/\sqrt{2}$ and $\beta=\pm 1/\sqrt{2}$, so as to make the average energy variance decrease in ${\cal O}(\Delta t^{2})$. If the energy variance continues to decrease closely to zero, then an energy eigenstate is approached. We remark that as demonstrated below it is not necessary to use the same ancillary state and time duration $\Delta t$ in every step of the procedure. Varying ancilla’s state away from $|\pm\rangle$ can be useful to avoid the system state to get stuck in states that have $\langle\hat{h}^{3}\rangle-\langle\hat{h}^{2}\rangle\langle\hat{h}\rangle=0$. See also below in Sec. IV for further discussions on this.

Random-walk approach. As the procedure outputs a pure state if the input is also pure, a question arises as to how we can analytically understand how the system is eventually driven to an eigenstate? Let us analyze the post-measurement states $|\psi^{\prime}_{m}\rangle$ by expanding it to the second order in $\Delta t$,

[TABLE]

We can rewrite the above equation to find the exponentiated action on $|\psi\rangle$, i.e., $|\psi^{\prime}_{m}\rangle\sim e^{\hat{P}_{m}}|\psi\rangle$ and ignore the overall constant. As shown in Appendix B, we find that to the second order in $\Delta t$

[TABLE]

As $\hat{P}_{m}$ is a polynomial of $\hat{h}$, one can separate it into two commuting parts: one Hermitian and the other anti-Hermitan, $\hat{P}_{m}=(\hat{P}_{m}+\hat{P}_{m}^{\dagger})/2+(\hat{P}_{m}-\hat{P}_{m}^{\dagger})/2=:\hat{Q}_{m}+i\hat{R}_{m}$. As the part $i\hat{R}_{m}$ is anti-Hermitian, its corresponding action $e^{i\hat{R}_{m}}$ is a unitary, and it does not modify the relative weight in the decomposition of energy eigenstates, so we can ignore it when we consider eigenstate projection. Thus, we focus on $|\psi^{\prime}_{m}\rangle\sim e^{\hat{Q}_{m}}|\psi\rangle$, where $\hat{Q}_{m}=(\hat{P}_{m}+\hat{P}_{m}^{\dagger})/2$.

After a long sequence of iterations, we will have a long product of operators $e^{\hat{Q}}$’s (which commute with one another) acting on the initial state $|\psi\rangle$, such as

[TABLE]

which looks like a sequence of ‘random walk’ using the two operators in the exponent. However, the key difference from a typical random walk is that there is a quantum state that changes after every step and the probability of moving to the left or right $p_{0/1}$ is state dependent, as in Eq. (2).

Here, as an approximation for the average action $e^{\overline{\hat{Q}(\alpha,\beta,\Delta t)}}$, which is valid in the limit $\Delta t\rightarrow 0$, we ignore the subsequent state dependence and use the initial $p_{0/1}(\psi)$ to calculate the average in the exponent: $p_{0}(\psi)\cdot\hat{Q}_{0}+p_{1}(\psi)\cdot\hat{Q}_{1}$, and we arrive at

[TABLE]

Thus, the average one-step action gives rise to a map on the system:

[TABLE]

The factor $c(\alpha,\beta)\geq 0$ is defined earlier and is maximized with a value $1/4$ when $\beta/\alpha=e^{i\phi}$ and $|\alpha|=|\beta|=1/\sqrt{2}$. This represents the optimal choice of ancillary parameters to maximize the converge rate, consistent with results presented earlier.

The meaning of the above equation is that the procedure tends to suppress components of eigenstates that have eigenvalues further away from $h_{\psi}\equiv\langle\psi|\hat{h}|\psi\rangle$. As one repeatedly applies the primitive, the state $|\psi\rangle$ itself will change and hence so will the expectation value $\langle\psi|\hat{h}|\psi\rangle$, with the latter eventually approaching the energy eigenvalue and the system state approaching the corresponding eigenstate. The random-walk analysis gives similar conclusion as that by the change in the average energy variance. We note that when ${\rm Im}(\alpha^{*}\beta)=0$ and the ancilla’s state not being $|\pm\rangle$, we need to carry out the expansion to the fourth order, but we do not perform the calculation here.

III Procedure for classical simulations

The primitive looks similar to the Hadamard test and consists a controlled-unitary action on the ancilla and the system, as well as a subsequent measurement on the ancillary qubit. Since the effect is to update the state vector of the system, for classical simulations of this process, we only need to compute two (un-normalized) wave functions $|\psi_{m}^{(k)}\rangle$ and their norm squares $p_{m}^{(k)}\equiv\langle\psi_{m}^{(k)}|\psi_{m}^{(k)}\rangle$ at each step, say, $k$-th,

[TABLE]

given the state, $|\psi^{(k-1)}\rangle$, of the system from the end of the previous step, the parameters $\alpha_{k}$ and $\beta_{k}$, and the unitary $U_{k}(\Delta t_{k})=e^{-i\Delta t_{k}\hat{h}}$.

One then decides to update the state $|\psi^{(k)}\rangle=|\psi^{(k)}_{m}\rangle/\sqrt{p_{m}}$ by choosing $m=0$ or $m=1$ with probability $p_{m}^{(k)}$. With a suitable choice of $\{(\alpha_{k},\beta_{k})\}$ and $\{\Delta t_{k}\}$, the long-iterated state $|\psi^{(k\gg 1)}\rangle$ will converge to some eigenstate $|E_{n}\rangle$, as illustrated below.

Simulating this procedure for spectral projection also provides us a quantum-inspired classical algorithm to obtain (randomly) excited states, whose accuracy does not depend on other lower lying levels. The costly part is applying $e^{-i\hat{h}\Delta t}$ to a state vector. However, for the purpose of a short-range Hamiltonian, one can use the Trotter decomposition and the individual $e^{-i\hat{h}_{j}\Delta t}$ from $\hat{h}=\sum_{j}\hat{h}_{j}$. Tensor-network representations can also be useful.

To obtain the entire set of eigenstates, we need to simulate the spectral projection as many times as the Hilbert space dimension. One can start with the system in an arbitrary initial $|\psi_{0}\rangle$ state. Run the procedure to obtain some eigenstate $|\phi_{1}\rangle$, then subtract the portion of $|\phi_{1}\rangle$ from $|\psi_{0}\rangle$: $|\psi^{[1]}\rangle=|\psi_{0}\rangle-\langle\phi_{1}|\psi_{0}\rangle|\phi_{1}\rangle$ and use the normalized version of $|\psi^{[1]}\rangle$ as the input of the procedure. Repeat this until one exhausts all the eigenstates that have nonzero overlap in $|\psi_{0}\rangle$. For the remaining eigenstates having zero overlap with $|\psi_{0}\rangle$, we can generate another random state and remove the components of all previously found eigenstates and use the resultant state as the input. In this way, we can eventually exhaust all eigenstates. The benefit of this method is that the accuracy of each eigenstate is independent of one another.

III.1 Simulation results: illustrative examples

Let us illustrate the algorithm by considering the system to be five-level, i.e. qudit with $d=5$. We generate a $5\times 5$ Hermitian matrix $H_{5}$, with $(H_{5})_{ij}=(H_{5})_{ij}^{*}=x+yi$ and $x$ and $y$ uniformly sampled from the range $[-1,1]$ (except $y=0$ for the diagonal elements). Here we only display its elements in the diagonals and below:

[TABLE]

whose eigenvalues $E_{i}$’s, sorted from the smallest to largest, are $\{-1.51593,-0.700576,0.388005,1.0888,2.51793\}$. We also randomly generate a 5-component normalized vector to be the initial state,

[TABLE]

The state $|\psi_{0}\rangle$ has an expected energy being $-0.525913$, with the probabilities $|\langle E_{i}|\psi_{0}\rangle|^{2}$ in the five eigenstates being, respectively,

[TABLE]

Next, we explore various combinations of $(\alpha,\beta)$ and $\Delta t$ in our classical simulations. Given that the optimal choice of $(\alpha,\beta)$ is such that $|\alpha\beta|=1/2$, i.e., within the one-parameter family $(1,e^{i\phi})/\sqrt{2}$, we first discuss the choice of the phase $\phi$ in this family. We have carried our a few simulations and displayed the results in Fig. 2.

(I) Iterations with fixed $\Delta t$ and $\phi$. With sufficient number of iterations, even fixing $\Delta t=1.0$ and $\phi=0$ (i.e. the standard Hadamard test), eigenstates can be reached with increasing accuracy as the number of iterations increases. In the simulations, we terminate the iteration once the energy variance $\langle(\Delta h)^{2}\rangle$ has reached below $10^{-10}$. The specific example run takes as long as 365 steps and converges to the eigenenergy $E_{3}=0.388005$.

(II) Iterations with fixed $\Delta t$ but $\phi$ from a given set. Bying fixing $\Delta t=1.0$ but choosing $\phi$ from $k\pi/4$ (with $k=0,1,\dots,7$), the specific example run takes 67 steps to converge to the eigenenergy $E_{2}=-0.700576$.

(III) Iterations with fixed $\Delta t$ but random choice of $\phi\in[0,2\pi)$. In the previous choice, $\phi$ is chosen from a set of values, here we consider choosing $\phi$ randomly from $[0,2\pi)$. In the example run, it takes 68 steps to converge to the eigenenergy $E_{1}=-1.51593$.

(IV) Iterations with varying $\Delta t$. In the previous three cases, we do not need to change $\Delta t$. But by allowing $\Delta t$ to vary, the efficiency can be improved. For example, by recycling $\Delta t$ from the set $\{100,100/3,100/3^{2},100/3^{3},100/3^{4},100/3^{5}\}$ and using random $\phi$, it takes 43 steps to converge to the eigenenergy $E_{5}=2.51793$.

We have also repeated the simulations using the above four types of choices but for $\alpha=\sqrt{3}/2$ and $\beta=e^{i\phi}/2$, and the results are shown in Fig. 3. We see that even without the optimal ancillary parameters, spectral projection can still be achieved. The steps it take to converge are not significantly larger than those using the optimal choice of the ancilla.

Let us compare our procedure to the QPE, in which the control-unitary needs to go as large power as $c-U^{2^{t_{g}-1}}$, in order to gain accuracy in $m$ binary digits, i.e. accurate up to $2^{-m}$, where $m=t_{g}-\log(2\epsilon+1/2\epsilon)$ and $1-\epsilon$ is the lower bound on the success probability of the QPE. To achieve the same accuracy as $2^{-33}\approx 10^{-10}$ in spectral projection by the QPE, one needs the number of ancillary qubits $t_{g}$ to be more than 33, and the power in $U$ differs in magnitude by $2^{33}$. In contrast, the ratio of the largest $\Delta t$ to the smallest used in our simulation (IV above) is only $3^{5}\approx 2^{8}$. In the above (I)-(III), $\Delta t$ is fixed, but it takes more steps to converge to the desired accuracy.

In the QPE, the power of the unitary $U^{2^{k}}$ needs to be precise in order for the algorithm to work. The procedure that we propose here does not require precise $\Delta t$. We have tested that the ability for the spectral projection does not depend on the precise values of $\Delta t$ as above, and other sequences can be used. For example, a different sequence is used in Fig. 4 as an example by perturbing the previous set of $\Delta t$, and spectral projection is still achieved.

In all of the above simulations, in addition to the energy value, the energy variance $\langle(\Delta h)^{2}\rangle$ is also recorded as the procedure is carried out. We have seen that on average, the energy variance indeed decreases.

III.2 Distribution of eigenstates: the Born rule

Given that the iterations based on the primitive in Sec. II lead to a procedure for projecting a system to eigenstates of an Hermitian operator $\hat{h}$, here, we investigate the distribution of eigenstates when this procedure is repeated many times. We again take the $H_{5}$ Hamiltonian (21) and the same initial state (22) of the system and carry out simulations for our spectral projection algorithm.

As seen from the results in Fig. 5 using 10,000 repetitions of the procedure, the distribution of the eigenstates agrees well with the Born rule, which predicts that the probability of the $n$-th eigenstate $|E_{n}\rangle$ is $p_{n}=|\langle E_{n}|\psi_{0}\rangle|^{2}$. That the Born rule applies can be explained as followed. Since the controlled unitary ${c-U}=|0\rangle\langle 0|\otimes\openone+|1\rangle\langle 1|\otimes e^{-i\Delta t\hat{h}}$ commutes with the Hamiltonian $\hat{h}$ of the system, and hence with any eigenstate projector $|E_{n}\rangle\langle E_{n}|$. Therefore the expectation value of the observable $|E_{n}\rangle\langle E_{n}|$ must be conserved and equals $|\langle E_{n}|\psi_{0}\rangle|^{2}$. Under the assumption and as observed above that the procedure leads to eigenstate projection, then the distribution $\{q_{n}\}$ of the projected eigenstates should remain the same as the initial distribution, i.e. $q_{n}=|\langle E_{n}|\psi_{0}\rangle|^{2}$.

We note that there is nothing special about the Hamiltonian $\hat{h}$, and our proposed algorithm works for any Hermitian operator. The Born rule will also apply. One may also regard our procedure as a method to realize the statement in the measurement postulate.

Number of iterations. In addition to the Born rule, we also investigate how many iterations are needed to reach a desired accuracy, e.g. $\langle(\Delta h)^{2}\rangle<10^{-10}$. In the same simulation for the study of the Born rule above, we also keep track of the number of iterations in each run it takes to reach that accuracy. The results are shown in Fig. 6 using histograms. As observed, the number of required iterations is not narrowly peaked and this reflects the randomness in the ancilla measurement outcome and the state dependence in the outcome probability.

IV Spectral projection algorithm applied to the transverse-field Ising model

Here we consider physical models, such as the Ising model in a transverse field (with the periodic boundary condition),

[TABLE]

Our parameterization is slightly different from that in the literature. The spin-spin coupling strength is $J=g$ (antiferromagnetic if $J>0$) and the external field is $B=(1-g)$. In Fig. 7, we take $g=0.5$ (the critical point in the large $N_{q}$ limit), $N_{q}=5$ and the initial state $|\psi_{0}\rangle=|+-+-+\rangle$ and simulate the spectral projection procedure. The values of $\phi$ are randomly chosen and those of $\Delta t$ are listed in the caption. We see in Fig. 7 that spectral projection can be achieved with accuracy of $10^{-10}$ by using $\Delta t$ that ranges less than three orders of magnitude. To use the QPE for spectral projection, it will require the unitary controlled by the ancilla to raise to at least $2^{32}$, which is far from practical at present.

We also compare the distribution of projected eigenstates in the simulation with the ideal Born rule. In the case of degeneracy, we assign the portion according to the overlap square with these degenerate eigenstates. This again confirms the Born rule of our spectral projection procedure in a spin model.

The use of the ancillary state $|A\rangle=|\pm\rangle$. In our simulations for the Ising model we have encountered cases where the use of $|\pm\rangle$ in the ancillary state has caused the system to flow to certain class of states which under further iterations do not change the energy, despite that they were not eigenstates. But we have not observed such phenomena in the random Hamiltonian case explored earlier. This can be explained by the expressions in Eq. (10), which shows that when $\langle\hat{h}^{3}\rangle-\langle\hat{h}^{2}\rangle\langle\hat{h}\rangle=0$, the energy does not change to lowest order in $\Delta t$. This occurs when

[TABLE]

as one can verify that

[TABLE]

and, hence, the above condition is satisfied. The state does change under the iteration, but not the magnitudes $|a_{i}|^{2}$ and $|b_{j}|^{2}$ (after proper normalization). In the case of the transverse-field Ising model, there are eigenstates of opposite energies, and thus it can happen the system is driven to states of the form (25). In our example from random Hermitian matrices, there are no eigenstates of opposite energies.

V Effect of decoherence

Our method in general does not protect against decoherence. Let us consider a simple depolarizing channel $D(\rho_{i})=(1-\epsilon)\rho_{i}+\epsilon\openone_{i}/2$ apply to every system qubit. Here we assume the ancillary control qubit is relatively error free, and this reminds us of the assumption in the so-called DQC1 quantum computing model Knill1998 , where only one qubit is clean. Due to the depolarizing channel, the state remains in the original un-decohered state with a probability approximately $(1-\epsilon)^{N_{q}}$, where $N_{q}$ is the number of qubits in the system, but the remaining portion $1-(1-\epsilon)^{N_{q}}$ can contribute substantially to the energy change and its variance. If the decoherence is applied at each step, then our procedure will have an error of order at least $1-(1-\epsilon)^{N_{q}}\approx N_{q}\epsilon$ for small $\epsilon$ at each step of the iteration. This is confirmed in our numerical simulations, as shown by the record of energy variance in Fig. 8.

However, we imagine a contrived scenario that the depolarizing channel acts only, e.g., every 30 steps. Then the procedure can achieve better accuracy in between two strikes of the decoherence. This is illustrated in Fig. 9. The ‘disruptions’ due to decoherence are visible, especially in terms of the upward jump in the energy variance in Fig. 9b. If there can be sufficient number of iteration steps carried out before decoherence takes place, then the system can converge close to an eigenstate. Of course, the depolarizing channel takes the system out of the eigenstate and the subsequent iterative spectral projection procedure may take the system towards another eigenstate. Since the decoherence process does not commute with the system’s Hamiltonian, our procedure in the presence of decoherence may serve at best a robust way of finding arbitrary eigenstates, rather than a robust way of spectral projection.

VI Spectral projection as a subroutine in the quantum annealing algorithm

We begin by describing the idea of quantum annealing and related algorithms. One of the first proposed quantum annealing methods is to use imaginary-time Schrödinger’s equation proposed by Finnila et al. FirstQA . The one that is close to the modern AQC Adia1 ; Adia2 is proposed by Kadowaki and Nishimori KadowakiNishimori , where the Hamiltonian is the combination of the time-independent Ising model and a time-dependent transverse field. The evolution of the quantum state was discussed in terms of real-time Schrödinger’s equation that takes the system in the ground state of the large-field limit towards that of the zero-field limit. The AQC similarly has a Hamiltonian $H(g)$ that interpolates between a simple Hamiltonian $H(g=0)$ with an easily prepared ground state $|G(0)\rangle$ and the final Hamiltonian $H(g=1)$ that encodes the solution of certain problem in the ground state of $H(g=1)$. Provided the minimum gap of $H(g)$ is not too small, then evolving under the Hamiltonian via a suitable path $g(t)$ will take the initial ground state very close to the final ground state at the end of the evolution,

[TABLE]

where $\hat{T}$ indicates that the integration is time-ordered, and $T$ is the total time duration.

The key idea of the QSA by Somma et al. Somma2008 is to exploit the quantum Zeno effect and replace the unitary evolution by measurement in the eigenbasis of $H(g_{i})$, in a successive sequence of discrete $g_{i}$ ($0<g_{1}<g_{2}<\dots<g_{T}=1$). If the overlap of successive ground state $|\langle G(g_{k})|G(g_{k+1})\rangle|\geq 1-\mu^{2}$ is sufficiently close to unity, then by the quantum Zeno effect, the final state after the whole sequence of measurement should be very close to the final ground state $|G(g=1)\rangle$. The standard QPE and a randomization procedure were proposed in Ref. Somma2008 to achieve the measurement approximately. Below, we use our spectral projection algorithm for the measurement in the QSA and perform classical simulations for two different Hamiltonians, and we loosely refer to this also as quantum annealing.

VI.1 Transverse-field XzY model

Here, we consider a different spin chain DegerWei2019 than the Ising model:

[TABLE]

One reason of choosing this transverse-field XzY model is because, for the qubit number $N_{q}$ being odd, there is a crossing in the lowest few energy levels when the parameter $g$ is varied; see e.g. Fig 10. But for $N_{q}$ being even, there is a small gap above the ground state (for finite $N_{q}$). Therefore, it is interesting to compare the two different cases (but in the same model) for the quantum annealing. In our simulations, we will take $r=0.5$.

We begin with the initial state either as the ground state or the highest-energy state of $H_{\rm XzY}(g=0,r=0.5)$ and run the simulations for the quantum annealing with our spectral projection algorithm as a subroutine. The projection subroutine works by performing 180 times the primitive in Sec. II, thereby approximately projecting the system to eigenstates of $H_{\rm XzY}(j\Delta g,r=0.5)$, where in this simulation $\Delta g=0.05$ and $j$ successively goes from 1 to 20, reaching $g=1$ at the end. We see, in Fig. 10a with $N_{q}=5$, that the quantum annealing does not work as there is an energy level crossing and the state of the system follows its path smoothly in the energy space crossing the lowest energy curve, and similarly for the initial highest-energy case. However, the quantum annealing indeed does work when there is a gap throughout the range of $g$ (except at the end) for the $N_{q}=6$-qubit case in Fig. 10b.

VI.2 Transverse-field Ising model

Here, we return to the transverse-field Ising model (24) and perform the quantum annealing with our spectral projection as a subroutine. The ground state at $g=0$ is unique and is given by $|0^{\otimes N_{q}}\rangle$. But the ground states at $g=1$ are doubly degenerate, and they are $|+^{\otimes N_{q}}\rangle$ and $|-^{\otimes N_{q}}\rangle$. Similar to the previous section, we examine small system sizes with $N_{q}=5$ and $N_{q}=6$, shown in Fig. 11. Given that there is small gap in both cases, the quantum annealing works.

In the above simulations we have used the $e^{-i\hat{h}\Delta t}$ without decomposing it into Trotter terms. In order to simulate larger systems, we separate the Hamiltonian into two parts: $H_{\rm e}(g)$ and $H_{\rm o}(g)$ for even and odd bonds, where terms in $H_{\rm e}$ commute with one another and similarly for the terms in $H_{\rm o}$. Thus we can apply a Trotter-Suzuki decomposition to $e^{-i(H_{\rm o}+H_{\rm e})\Delta t}\approx e^{-iH_{\rm o}\Delta t}e^{-iH_{\rm e}\Delta t}$. This simulates the scenario that in the quantum circuit one can apply simultaneously the commuting terms of the controlled version of $e^{-iH_{\rm o}\Delta t}$ and subsequently those of $e^{-iH_{\rm e}\Delta t}$. In our classical simulations, we use a 4-th order Trotter-Suzuki decomposition Sornborger1999 ; DhandSanders for $e^{-i(H_{\rm o}+H_{\rm e})\Delta t}$:

[TABLE]

where $a_{1}=(2+\sqrt{2})/4$, $a_{2}=-a_{1}$, and $a_{3}=(1+\sqrt{2})/2$.

As seen in Fig. 12 with $N_{q}=114,16,18,20,\&22$, the annealing proceeds at initializing the state at the ground state of $H_{\rm TFI}(g=0)$, which is $|00\dots 0\rangle$. Then the spectral projection is applied successively at $g=j\Delta g$ for $j=1,2,\dots,40$ and $\Delta g=0.025$ (with the primitive being run 210 times in each projection procedure), ending at $g=1$ at the end of the annealing. The final energy after the anneal is seen to be close to the final ground-state energy, which is $-N_{q}$. The accuracy in this case can be increased by making the $\Delta g$ smaller and total number of iteration steps larger. The energy variance is generally the largest around $g=0.5$, and this is expected as, in the thermodynamic limit, there is a second-order quantum phase transition at $g_{c}=0.5$, and it is known that the gap closes as ${\cal O}(1/N)$ when $g$ approaches $g_{c}$ from below. As $g$ approaches 1, the ground state becomes doubly degenerate.

VI.3 Effect of decoherence in the annealing

Here, we take into account of the decoherence effect in our spectral projection and discuss how it affects the quantum annealing. In general, our algorithm does not project against decoherence, as discussed in Sec. V, and hence the resulting quantum annealing will be worse than the noise-free case. The degree of inaccuracy depends on the error rate $\epsilon$. We use, as an illustration, the contrived scenario discussed above that the decoherence with $\epsilon=0.01$ occurs at every 31th step in our spectral projection subroutine, in which the primitive is run for 210 steps. We test this on the 5-qubit transverse-field Ising model ${H}_{\rm TFI}(g)$ and the results of three different runs for the quantum annealing are shown in Fig. 13. As opposed to the noise-free case, there is some probability (depending on the noise rate and strength) that the final state may end up far from the final ground state. But there is also some probability that the final state is close to the final ground state. Developing noise-protecting spectral projection is thus a desirable goal that can yield a noise-protecting quantum annealing algorithm.

VII Concluding remarks

We have proposed a quantum algorithm for projecting to eigenstates of any Hermitian operator, provided one can access the associated control-unitary evolution and measurement of the controlling ancilla qubit. The procedure is iterative and the distribution of the projected eigenstates obeys the Born rule. It is robust against imprecision in timing. But it has only limited resilience against decoherence; the iterative procedure takes the system towards eigenstates, even after the influence of decoherence such as a depolarizing. It has no capability of error correction or prevention. We view our method as a simpler algorithm to project the system into eigenstates of a Hermitian observable than the standard QPE and it can also be used to extract eigenvalues. We compare our spectral decomposition to the standard QPE and an iterative version in Table 1. Our algorithm can be used as a subroutine in the quantum annealing procedure by measurement Somma2008 to drive to the ground state of a final Hamiltonian. We have performed simulations that demonstrate the utility of our algorithm. We note that a previously proposed scheme of ground state cooling quantum computation also uses ancilla measurement for the cooling GSCQC . Our scheme uses ancilla measurement for the spectral projection and the way it is used in the QSA is similar to the quantum Zeno effect. It will be useful to develop a noise-protecting spectral projection. A proof-of-principle demonstration of our spectral projection algorithm on currently available quantum computers will also be desirable.

Post-selection allows projection to the ground state, but the probability for obtaining the desired post-selected outcome is exceedingly small. The algorithm that we have attempted for the imaginary time evolution suffers some problems that make it not practical; see Appendix A. The fact that we end up with a spectral projection that obeys the Born rule seems to indicate that we may need to go beyond the primitive used in this paper to achieve an imaginary-time evolution quantum algorithm, as done in Ref. Motta2019 . But whether imaginary-time evolution can be achieved without using an effective Hamiltonian is an interesting question to consider.

In order to classically simulate our spectral projection algorithm, it will generally take exponential time in the number of qubits of the system, as one needs to compute $e^{-i\hat{h}\Delta t}|\psi\rangle$. Thus, it will be interesting if such a procedure can be carried out in a quantum computer for system sizes beyond the capability of classical simulations. This might be a useful playground for demonstrating quantum advantage.

Acknowledgements.

This work was supported by National Science Foundation under grants No. PHY 1620252 and No. PHY 1915165. T.-C.W. acknowledges useful discussions with Fernando Brandão, David Poulin and Barry Sanders. We also thank an anonymous referee for his/her suggestions that help improve the original manuscript.

Appendix A A failed attempt to construct an adaptive procedure for imaginary-time evolution

We have considered a primitive similar to the Hadamard test, except using an ancillary state of $\alpha|0\rangle+\beta|1\rangle$. The controlled unitary gate is of the form $c-e^{-i\hat{h}\Delta t}$. We consider the post-measurement state of the system up to the first order in $\Delta t$,

[TABLE]

Our motivation here is to achieve the nonunitary action $e^{-\hat{h}\Delta\tau}$ on $|\psi\rangle$, which, to first order, is $[1-\hat{h}\Delta\tau]|\psi\rangle$. Let us choose to make it work for the $m=0$ outcome by requiring that

[TABLE]

then the nonunitary action is achieved, i.e., the effective action on the system is (ignoring normalization)

[TABLE]

obtaining an effective time step $\Delta\tau=\Delta t/r$ in the imaginary-time evolution. To satisfy Eq. (33), $\alpha$ and $\beta$ can be taken as

[TABLE]

and the probability for each outcome (without approximation) is

[TABLE]

The Pauli X measurement on the ancilla can be realized by first applying the Hadamard gate $H$ before measuring in the standard Z basis; see Fig. 1.

However, for the outcome ‘1’, the system will be collapsed to an undesired state, to the first order in $\Delta t$,

[TABLE]

The second term is not harmful, as by applying to the post-measurement state the ‘correcting’ unitary

[TABLE]

the system becomes

[TABLE]

to the first order in $\Delta t$. We note that this additional step is not necessary as it only modifies the relative phases of different eigen-components, but not the amplitudes.

The second term inside the bracket of Eq. (39) and Eq. (34) represents the step size of a random walk in the exponent of an action $e^{-\hat{h}\Delta\tau_{i}}$ on a quantum state $|\psi\rangle$, where $\Delta\tau_{0}=\Delta t/r$ and $\Delta\tau_{1}=-\Delta t/(r+4/r)$ for the two respective measurement outcomes ‘0’ and ‘1’; see Fig. 1b for illustration. The corresponding probabilities (36) are approximately,

[TABLE]

where $h_{{\psi}}\equiv\langle{\psi}|\hat{h}|{\psi}\rangle$ is the average energy for the state $|{\psi}\rangle$ of the system prior to this iteration. The dependence of $p$’s on the system state $|\psi\rangle$ prevents us from getting a closed-form expression for the outcomes of a long sequence of iterations.

By post-selecting the ‘0’ outcome in the primitive, and by repeating this one $n$ times we can achieve exponential decay to the ground state, via

[TABLE]

Imaginary time evolution is employed in many classical numerical methods, such as the iTEBD method for ground states Vidal2004 . However, for our quantum procedure the desired branch of having all ‘0’ outcomes occurs with an exponentially small probability, so it is not very useful in practice.

Instead of postselection, one may perform an additional operation if the undesired outcome ‘1’ occurs. We have attempted such idea but we did not succeed. What is described below is such a failed attempt.

Let us define one iteration to be the process from entangling the system with an ancilla to measuring the ancilla and possibly correcting with the unitary if needed. If the first step yields the ‘0’ outcome, then one arrives at the desired imaginary-time evolution Eq. (34). We ask what can one do if one obtains the ‘1’ outcome and arrives at a state in Eq. (39)? We can proceed with a second iteration by choosing a different parameter $r^{\prime}$. The desired outcome ‘0’ after this iteration would put the system in the state

[TABLE]

If we choose $r^{\prime}$ such that ${1}/{(r+4/r)}-{1}/{r^{\prime}}=-{1}/{r}$, i.e.,

[TABLE]

then the outcome ‘0’ leads to the desired imaginary-time evolution Eq. (34).

However, if instead the measurement still gives the undesired outcome of ‘1’, we need to correct it further by repeating the iteration until outcome ‘0’ is obtained by choosing the parameter $r_{n+1}$ in the $(n+1)$-th round via

[TABLE]

and we terminate the iteration when ‘0’ outcome is obtained. Then the desired one-step imaginary-time evolution will give

[TABLE]

This procedure is summarized in Fig. 14.

However, this procedure suffers from the occurrence of long sequences of ‘1’ outcomes, as our simulations show. As a rough estimate by dropping the first-order contribution, the probability of $n$ successive ‘1’ outcomes is

[TABLE]

which does not decay exponentially. Figure 15 shows the values of $r_{n}$ with $r_{1}=1$. One can start with a larger $r_{1}$ so as to get a smaller ratio of $r_{n}/r_{1}$, but the scaling is still not exponentially small.

Appendix B Some derivations

B.1 Energy change

Let us list the post-measurement state,

[TABLE]

and the probability that it occurs,

[TABLE]

where it is convenient to define ${\cal R}_{1}$ and a related ${\cal R}_{h}$:

[TABLE]

Thus the change in energy is

[TABLE]

Then, expanding the above expression in series of $\Delta t$ is straightforward.

B.2 Change in average energy variance

In the main text, we have the expression for the average energy variance

[TABLE]

By expanding the square and explicitly summing over $m$, we obtain

[TABLE]

Then, expanding the above expression in series of $\Delta t$ is straightforward.

B.3 Average $\hat{Q}$ action

By expanding the post-measurement state $|\psi^{\prime}_{m}\rangle$ to the second order in $\Delta t$, we have

[TABLE]

The goal is to the above equation to the exponentiated form $|\psi^{\prime}_{m}\rangle\sim e^{\hat{P}_{m}}|\psi\rangle$, correct to the second order. Naturally, $\hat{P_{m}}$ will contain the second term in the square bracket. But we also need to take into account of other contribution to the second order. So we can set

[TABLE]

Expanding $e^{\hat{P}_{m}}$, we have to the second order,

[TABLE]

which should equals

[TABLE]

Therefore, we obtain

[TABLE]

From this, it is straightforward to obtain $\hat{Q}_{m}=(\hat{P}_{m}+\hat{P}_{m}^{\dagger})/2$ and perform the sum $\sum_{m=0,1}p_{m}\hat{Q}_{m}$. In the end, we arrive at

[TABLE]

Then, expanding the above expression in series of $\Delta t$ is straightforward.