Quantum power method by a superposition of time-evolved states

We propose a quantum-classical hybrid algorithm of the power method, here dubbed as quantum power method, to evaluate $\hat{\mathcal{H}}^{n}|\psi\rangle$ with quantum computers, where $n$ is a nonnegative integer, $\hat{\mathcal{H}}$ is a time-independent Hamiltonian of interest, and $|\psi\rangle$ is a quantum state. The quantum power method is formulated on the fact that the Hamiltonian power $\hat{\mathcal{H}}^{n}$ can generally be approximated by a linear combination of time-evolution operators $\hat{U}(t)=\mathrm{e}^{-\mathrm{i} \hat{\mathcal{H}}t}$ at $n+1$ different time ($t$) variables. The formalism is based on a time-discretized form of the higher-order derivative $\hat{\mathcal{H}}^{n}={\rm i}^n \mathrm{d}^n\hat{U}(t)/\mathrm{d} t^n|_{t=0}$ incorporated with the central-finite-difference scheme and the symmetric Suzuki-Trotter decomposition, by which the approximated Hamiltonian power retains its Hermiticity under a controlled accuracy. In the quantum power method, the Suzuki-Trotter decomposition is employed to evaluate each time-evolution operator $\hat{U}(t)$ separately at different time $t$ close to zero on quantum computers, while a linear combination of these time-evolution operators is treated on classical computers. The number of gates required for approximating $\hat{\mathcal{H}}^n$ is $O(N n)$, where $N$ is the number of qubits assuming a local Hamiltonian $\hat{\mathcal{H}}$. For an application of the quantum power method, we combine this method with a multireference Krylov-subspace diagonalization scheme and show, by noiseless numerical simulations for a spin-$1/2$ Heisenberg model with the Hamiltonian power $\hat{\mathcal{H}}^{n}$ up to $n=11$ (but not limited), that the estimated ground-state energy and the ground-state fidelity over a variational-quantum-eigensolver scheme is systematically improved with increasing the power $n$.


I. INTRODUCTION
Numerically solving quantum many-body systems is one of the most useful approaches for yet challenging issues in condensed matter physics and quantum chemistry [1][2][3]. With classical computers, a repeated multiplication of a Hamil-tonianĤ of interest to a properly chosen state, i.e., the power iteration, is an essential element of various practical and advanced numerical techniques such as Krylovsubspace methods [4] including the Lanczos method [5][6][7][8][9][10], and polynomial-expansion methods [11]. Such methods allow for calculating not only ground states but also dynamics [12][13][14][15][16] of quantum many-body systems. A major obstacle in these methods is, however, the exponential growth of the dimension of the Hilbert space with its system size N. The Lanczos method has been implemented also with the variational Monte Carlo technique to systematically improve variational states towards the exact ground state [17]. While the variational Monte Carlo method allows for substantially larger N than the full-Hilbert space approaches, an affordable number of the Lanczos iterations is practically limited to a few due to the O(N n ) number of terms constitutingĤ n .
Recently, simulating quantum many-body systems with quantum computers [18][19][20][21] attracts great interest due to experimental realizations of and advances on quantum de-vices [22][23][24][25][26][27][28][29][30][31]. Quantum computers will allow for a rather more direct access to quantum states defined in a Hilbert space of potentially huge dimensions that cannot be treated with classical computers. At present, quantum computers are prone to noises and computations have to be accomplished with a small number of gates. In this regard, the variational-quantum-eigensolver (VQE) scheme has been proposed to simulate quantum many-body systems using noisy intermediate-scale quantum devices [32] and classical computers in a hybrid manner [33][34][35][36][37][38][39]. VQE calculations with noisy quantum devices are now getting affordable for fairly larger systems [40] than in the earlier studies. While the majority of VQE schemes is devoted for gate reduction at the expense of the increased number of measurements, a measureand-reuse technique has been proposed for reducing the number of qubits [41]. Such a qubit-reuse technique has been demonstrated by evaluating the ground-state energy of the one-dimensional Heisenberg model accurately only with a few trapped-ion qubits [42].
Moreover, to bypass variational parameter optimization and ansatz-state cultivation inherent in the VQE scheme, several versions of Krylov-subspace methods have been proposed. The QLanczos method [43] generates a Krylov subspace by evolving a reference state with an approximate imaginarytime evolution [44][45][46]. The multireference-selected quan-arXiv:2008.03661v1 [quant-ph] 9 Aug 2020 tum Krylov algorithm generates a set of states spanning a Krylov subspace by evolving selected reference states in real time [47]. As a related method, a quantum version of the filter-diagonalization method has been developed [48]. A version of the inverse-iteration method suitable for quantum computers [49] makes use of an integral representation of the inverse of the Hamiltonian [50]. Recently, an implementation of the exact imaginary-time evolution with the help of ancillary qubits and Grover's search algorithm has been proposed [51]. Subspace-diagonalization schemes, with subspaces not restricted to a Krylov subspace but intended to approximate a particular set of eigenspaces of the Hamiltonian of interest, have been implemented for calculating not only the ground state but also excited states of correlated quantum-chemistry systems [52][53][54][55][56][57].
In this paper, we propose a quantum power method, a version of the power method suitable for quantum-classical hybrid computing of quantum many-body systems. The method is based on a time-discretized form of the higher-order deriva-tiveĤ n = i n d nÛ (t)/dt n | t=0 of the time-evolution operator U(t) = e −iĤt , by which the Hamiltonian powerĤ n is represented as a linear combination ofÛ(t) at different time (t) variables close to t = 0. The approximated Hamiltonian power retains its Hermiticity by engaging the time-discretized formalism with a central-finite-difference scheme for the time derivatives and the symmetric Suzuki-Trotter decomposition of the time-evolution operators. Assuming that the Hamilto-nianĤ is local, the number of the gates required for approx-imatingĤ n in the quantum power method is O(Nn) where N is the system size (i.e., the number of qubits). We numerically demonstrate that the quantum power method can control the systematic errors, due to the finite-difference scheme for the time derivatives and the Suzuki-Trotter decomposition of the time-evolution operators, in approximating the Hamiltonian powerĤ n with n as large as 100 for N up to 24. We apply the quantum power method to generate a Krylov subspace and perform, using noiseless numerical simulations, the multireference Krylov-subspace diagonalization for a onedimensional spin-1/2 Heisenberg model with various reference states including those obtained by the VQE scheme. We find that the estimated ground-state energy as well as the ground-state fidelity are significantly improved with increasing the power n, thus providing a way to systematically improve the VQE scheme. We also briefly outline other applications of the quantum power method.
The rest of the paper is organized as follows. In Sec. II, we introduce the spin-1/2 Heisenberg model on a onedimensional periodic chain. In Sec. III, we formulate the quantum power method by describing the central-finitedifference scheme for the time derivatives, basic properties of the approximated Hamiltonian power, the Richardson extrapolation to systematically eliminate the systematic errors, and the Suzuki-Trotter decomposition for the time-evolution operators. In Sec. IV, we review the Krylov-subspace diagonalization scheme for an application of the quantum power method. In Sec. V, we first numerically demonstrate that the systematic errors in the quantum power method are well controlled to be essentially exact and then show numerical results of the Krylov-subspace diagonalization combined with the quantum power method for the spin-1/2 Heisenberg model. The paper is summarized with discussions in Sec. VI. Explicit forms of the higher-order symmetric Suzuki-Trotter decompositions and their error analysis are provided in Appendix A. An alternative formalism of approximating the Hamiltonian power is discussed in Appendix B. For other applications of the quantum power method, some properties of the moments and cumulants are discussed in the context of the quantum power method, and the connected-cluster expansion (CMX) for the imaginary-time evolution is demonstrated by numerical simulations in Appendix C. The Lanczos method with an emphasis on its connection to the moments is also described in Appendix D. Throughout the paper, we set = 1.
where i + 1 in the subscript should be read as 1 if i = N due to the periodic-boundary conditions. For the later use in the Suzuki-Trotter decomposition of the time-evolution operator, we divide the Hamiltonian into two parts aŝ Notice that [P 2i,2i+1 ,P 2 j,2 j+1 ] = [P 2i−1,2i ,P 2 j−1,2 j ] = 0, where [Â,B] =ÂB −BÂ is the commutator of two operatorsÂ and B.

III. FORMALISM
In this section, we formulate the quantum power method. Figure 1 illustrates an overview of the formalism based on the higher-order derivative of the time-evolution operator, which is decomposed approximately using the symmetric Suzuki-Trotter decomposition.

A. Hamiltonian power as a linear combination of unitary time-evolution operators
The time-evolution operatorÛ(t) of the time-independent HamiltonianĤ at time t is given bŷ where t is real. The nth power of the Hamiltonian,Ĥ n , is thus given by the nth derivative of the time-evolution operator at t = 0, i.e.,Ĥ n = i n d nÛ (t) dt n t=0 .
By introducing a small time interval ∆ τ , we replace the time derivative in Eq. (8) with the central finite-difference aŝ whereĤ n (∆ τ ) = n k=0 c n,kÛ n 2 − k ∆ τ (10) and A derivation of the coefficients c n,k is shown in Fig. 1(b). Equations (9) and (10) imply that the nth power of the Hamiltonian,Ĥ n , can be approximated with a controlled accuracy as a linear combination of the time-evolution operators at n+1 different time variables. From the unitarity of the time-evolution operator and its accordance with the time-reversed evolution, it follows that the approximated Hamiltonian powerĤ n (∆ τ ) is Hermitian and an even function of ∆ τ i.e., In the last equality, we have used that c n,k in Eq. (11) is an even (odd) function of ∆ τ when n is even (odd). SinceĤ n (∆ τ ) is an even function of ∆ τ , the systematic error E FD in odd powers of ∆ τ is absent in Eq. (9). Moreover, with the multiplication law of the time-evolution operatorÛ (t)Û (t ) =Û (t + t ), Eq. (10) can be written aŝ The last line in Eq. (14) indicates that the approximated Hamiltonian powerĤ n (∆ τ ) satisfies a law of exponentŝ Namely,Ĥ n (∆ τ ) is exactly the nth power ofĤ n=1 (∆ τ ) for n 0. In fact, Eq. (14) can be understood simply aŝ

B. Richardson extrapolation
As we shall discuss in Sec. V A, the systematic error E FD due to the finite-difference scheme for the time derivatives in Eq. (9) can be controlled by varying the time discretization step ∆ τ . However, it is often practically useful to reduce the systematic error E FD by not taking too small ∆ τ in the algorithmic level. The Richardson extrapolation can achieve a better error estimate by systematically eliminating lower-order errors in Eq. (9).
We can use the Richardson extrapolation recursively to further eliminate the O(∆ 4 τ ) error in Eq. (17). Namely, the rth order Richardson extrapolation can be obtained recurrently aŝ Overview of the quantum power method proposed here. (a) The Hamiltonian power H n is approximated as a linear combination of the time-evolution operators [Û(∆ τ /2)] n−2k for k = 0, 1, . . . , n, in which eachÛ(∆ τ /2) is further decomposed intoŜ (p) 2m (∆ τ /2) using the symmetric Suzuki-Trotter decomposition. Here, ∆ τ is a small time interval, and thus real positive number. E FD and E ST denote systematic errors due to the finite-difference scheme for the time derivatives and the symmetric Suzuki-Trotter decomposition of the time-evolution operators, respectively. (b) An illustration of the central-finite-difference scheme for the nth order derivative of the time-evolution operatorÛ(t) at t = 0. Pascal's triangle with an alternating sign in time t and power n provides coefficients c n,k of a linear combination of the time-evolution operators that approximates the Hamiltonian powerĤ n . The systematic error due to the finite-difference scheme is . The systematic error E ST due to the Suzuki-Trotter decomposition for approximating the Hamiltonian powerĤ n in (a) is O(∆ 2m τ ) because of the factor 1/∆ n τ in c n,k . D (p) 2m (= 2p m−1 + 1) is the circuit depth for a singleŜ (p) 2m (∆ τ ), and p is typically an O(1) integer parameter for the symmetric Suzuki-Trotter decomposition, independent of the number N of qubits. The figure refers to m = 1, p = 3, and N = 6. The rth order Richardson extrapolation improves systematically the systematic errors as ) at the expense of increasing the number (r + 1)(n + 1) of terms in the linear combination. This implies that the lowest-order symmetric Suzuki-Trotter decomposition with m = 1 is adequate to control these systematic errors consistently. The number of gates, indicated by small blue rectangles in (c), required to approximately represent the Hamiltonian power H n scales as O(Nn).
is a linear combination of n + 1 unitaries,Ĥ n (r) (∆ τ ) is a linear combination of (r + 1)(n + 1) unitaries. Note also thatĤ n (r) (∆ τ ) is no longer the nth power ofĤ n=1 . In our numerical simulations, we choose h = 2 when the Richardson extrapolation is used.
Three additional remarks are in order regarding the properties of the approximated Hamiltonian powerĤ n (∆ τ ). First, if a forward or backward, instead of central, finite-difference scheme is employed in Eq. (10), the Hermiticity and the even dependence on ∆ τ ofĤ n (∆ τ ) in Eq. (13) are both violated. Therefore, the central finite-difference scheme is a crucial choice. Second, when the time-evolution operatorÛ(∆ τ ) is approximated by a Suzuki-Trotter decomposition, the corresponding Suzuki-Trotter error E ST appears in Eqs. (10) and (14). Since the implementation of a higher-order Suzuki-Trotter decomposition on quantum computers requires many layers of gates, it is essential to control E ST with a lower order Suzuki-Trotter decomposition. Third, if a symmetric Suzuki-Trotter decomposition, which retains the equivalence between the inverse of the time evolution and the time-reversed evolution [the right-most equality in Eq. (12)], is employed to decompose the time-evolution operators in Eqs. (10) and (14), the resultingĤ n (∆ τ ) still satisfies the Hermiticity and the even dependence on ∆ τ . Therefore, it is important to adopt a symmetric Suzuki-Trotter decomposition (see Sec. III C 3 for details).

C. Suzuki-Trotter decomposition
The formalism so far is based on the exact time-evolution operatorÛ(t) in Eq. (7). However, on quantum computers, the time-evolution operator with its exponent composed of the sum of non-commuting operators usually has to be represented as a product of time-evolution operators with each exponent composed of the sum of commuting operators. For this purpose, the Suzuki-Trotter decomposition is employed to approximately decompose the time-evolution operator.
We should emphasize that one of the crucial steps for the successful quantum power method is to determine properly in which stage the time-evolution operators inĤ n (∆ τ ) should be approximated by the Suzuki-Trotter decomposition, either in Eq. (10) or in Eq. (14). Although Eqs. (10) and (14) are exactly the same if the exact time-evolution operators are used, they are no longer the same in general once the time-evolution operators are approximated. Therefore, there are at least two routes to formulate the quantum power method. As we shall discuss in details, these two approaches give us two different algorithms that scale differently in the power n. It turns out that when the power n is larger than four, the algorithm formulated on the basis of Eq. (14) with the lowest-order symmetric Suzuki-Trotter decomposition is preferable, otherwise the formalism based on Eq. (10) with the higher-order symmetric Suzuki-Trotter decompositions is favored in terms of the gate counts (see Appendix B).
To understand the difference of these two approaches, in this section, we briefly summarize a systematic construction of the higher-order symmetric Suzuki-Trotter decompositions [58][59][60] for the quantum power method.
Starting withŜ (p) 2 (∆ τ ) ≡Ŝ 2 (∆ τ ), the higher-order decom-positionŜ (p) 2m (∆ τ ) for m 2, satisfyinĝ can be constructed recursively aŝ wherek (p) , and p is an odd integer with p 3 [72]. The superscript "(p)" implies thatŜ (p) 2m consists of a product of pŜ (p) 2m−2 's. The parameter k (p) m is determined so as to eliminate the residual term x 2m−1R 2m−1 in lnŜ (p) 2m (∆ τ ) and thuŝ Namely, k (p) m is the solution of (p−1) k (p) the residual terms of even power such as x 2mR 2m are absent in the exponent ofŜ (p) 2m (∆ τ ) in Eq. (29), shown by the same argument for m = 1. Some of the higher-order symmetric Suzuki-Trotter decompositions are explicitly provided in Appendix A 1. As shown in Appendix A 2, the parameter p affects the accuracy of the decomposition for a given m. 3. Unitarity and time-reversed evolution ofŜ (p) 2m (∆ τ ) As implied in Eq. (30),Ŝ (p) 2m (∆ τ ) retains not only the unitarity but also the equivalence between the inverse and timereversed evolution, i.e., Therefore, the Hermiticity and the even dependence on ∆ τ of H n (∆ τ ) in Eq. (13) are both retained even when the exact time-evolution operators in Eqs. (10) and (14) are approximated by simply replacing them withŜ (p) 2m 's. In contrast, an asymmetric Suzuki-Trotter decomposition F(∆ τ ), such asF(∆ τ ) = e xĤ A e xĤ B , results in Thus,F(∆ τ ) retains the unitarity but the inverse is no longer equivalent to the time-reversed evolution. In this case, either the Hermiticity or the even dependence on ∆ τ ofĤ n (∆ τ ) in Eq. (13) is violated if the exact time-evolution operators in Eqs. (10) and (14) are approximated byF's. For example, if we consider an operatorĤ H ( Therefore, the symmetric Suzuki-Trotter decom-positionŜ (p) 2m (∆ τ ) is essential for the resulting Suzuki-Trotter approximatedĤ n (∆ τ ) to retain both the Hermiticity and the even dependence on ∆ τ . Note that asymmetric Suzuki-Trotter decompositions and their connection to symmetric ones have been studied in Ref. [73].

Circuit depth for a single time-evolution operator approximated by the Suzuki-Trotter decomposition
We now consider the circuit depth D (p) 2m required for a single time-evolution operatorÛ(∆ τ ) approximated by the symmetric Suzuki-Trotter decompositionŜ (p) 2m (∆ τ ), as in Eq. (27) [also see Fig. 1(c)]. We define D (p) 2m as the number of noncommuting exponentials appearing inŜ (p) 2m (∆ τ ). For example, the depth of consists of a product of pŜ (p) 2m−2 's, the depth ofŜ (p) 2m (∆ τ ) without contracting commuting exponentials is pD (p) 2m−2 . However, sinceŜ (p) 2m (∆ τ ) involves p − 1 products of two consecutiveŜ (p) 2m−2 's, between which two commuting exponentials reside, p − 1 exponentials can be contracted. We thus obtain that D (p) . By using this relation recursively, we can find that Substituting D (p) 2 = 3 in Eq. (33) yields that Recalling that p is a typically O(1) integer parameter, the depth increases exponentially with m but is independent of the number N of qubits. Therefore, the lower-order Suzuki-Trotter decomposition is highly desirable to shallow the depth of a quantum circuit.

D. Quantum power method
While the time-evolution operators satisfy the multiplication lawÛ(∆ τ )Û(∆ τ ) =Û(∆ τ + ∆ τ ), this is no longer correct when the time-evolution operators are approximated by the Suzuki-Trotter decomposition, i.e.,Ŝ (p) 2m (∆ τ )Ŝ (p) 2m (∆ τ ) Ŝ (p) 2m (∆ τ + ∆ τ ). Therefore, it is crucial to carefully consider when the time-evolution operators in the approximated Hamiltonian powerĤ n (∆ τ ) should be replaced with the symmetric Suzuki-Trotter decomposition, either in Eq. (10) or in Eq. (14), implying that there exist two different routes to formulate the quantum power method. As we shall show here and in Appendix B, these two approaches provide two different algorithms of the quantum power method that differ in the scaling of complexity but control the systematic errors E FD and E ST with essentially the same accuracy. In this section, we formulate the quantum power method based on Eq. (14) that scales much better when the power n is large. In Appendix B, we describe an alternative algorithm formulated on the basis of Eq. (10), which is favored when the power n is small (n 4).
By incorporating the symmetric Suzuki-Trotter decomposition into the approximated Hamiltonian powerĤ n (∆ τ ) in Eq. (14), the Hamiltonian powerĤ n is finally approximated asĤ where O(∆ 2 τ ) represents the systematic error E FD due to the finite-difference scheme for the time derivatives, and O(∆ 2m τ ) denotes the systematic error E ST due to the Suzuki-Trotter decomposition of the time-evolution operators.Ĥ n ST (∆ τ ) is the central quantity in the quantum power method that approximates the Hamiltonian powerĤ n aŝ Here, we have used the fact that Ŝ (p) . It should also be noticed that the order O(∆ 2m τ ) of the Suzuki-Trotter error E ST in Eq. (35) is decreased by one from the naively expected order O(∆ 2m+1 τ ) as in Eq. (27), because of the factor 1/∆ n τ in c n,k . Equation (35) already reveals a remarkable advantage in the quantum power method formulated on the basis of Eq. (14): in order to control the systematic errors E FD and E ST with the same order of accuracy, it is enough to adopt the lowest-order Suzuki-Trotter decomposition with m = 1, independently of the power n. This is in sharp contrast to the other formalism based on Eq. (10), in which the order m of the Suzuki-Trotter decomposition has to be increased with the power n and, as described in more details in Appendix B, essentially the complexity increases exponentially with the power n.
Equation (37) indicates thatĤ n ST (∆ τ ) satisfies the law of exponentsĤ Moreover,Ĥ n ST (∆ τ ) is Hermitian and an even function of ∆ τ , i.e.,Ĥ indicating that the systematic error E ST in odd powers of ∆ τ is absent in Eq. (35). Therefore, recalling that the systematic error E FD in odd powers of ∆ τ is also absent (see Sec. III B), the Richardson extrapolation can eliminate the finite-difference error E FD and the Suzuki-Trotter error E ST simultaneously aŝ whereĤ n ST(r) (∆ τ ) is the rth order Richardson extrapolation of the approximated Hamiltonian power, i.e., is a linear combination of (r + 1)(n + 1) unitariesŜ (p) 2m . Equation (40) reveals another significant feature of the quantum power method that only polynomial resources with the lowest-order symmetric Suzuki-Trotter decomposition suffice to systematically and consistently eliminate the lower order systematic errors in E FD and E ST . In Sec. V A, we will show by numerical simulations that these systematic errors in the approximated Hamiltonian power are well controlled with the time-discretization step ∆ τ for the power n as large as 100.
For the application purpose of the quantum power method, it is important that the symmetry of the HamiltonianĤ is still respected in the approximated Hamiltonian powerĤ n ST(r) (∆ τ ). This is indeed the case in the quantum power method formulated here because Therefore, the symmetry of the HamiltonianĤ is preserved in the quantum power method within the systematic error E ST due to the Suzuki-Trotter decomposition that can be well controlled. Notice that there is no contribution from the systematic error E FD due to the finite-difference scheme of the time derivatives in the right hand side of Eq. (42) because Ĥ ,Ĥ n (r) (∆ τ ) = 0. Figure 1 summarizes the quantum power method formulated here. In the quantum power method, the Hamiltonian powerĤ n is approximated toĤ n ST (∆ τ ) represented as a linear combination of the n + 1 Suzuki-Trotter decomposed 2m (∆ τ /2)] n−2k is evaluated on quantum computers. As illustrated in Fig. 1(c), a quantum circuit for a singleŜ (p) 2m (±∆ τ /2) has the circuit depth D (p) 2m = 2p m−1 + 1, and thus the circuit depth required forĤ n ST (∆ τ ) is at most O(n) with a prefactor D (p) 2m ∼ O(1). Hence, assuming thatĤ consists of O(N) local terms, the number of gates required forĤ n . When the rth-order Richardson extrapolation is employed, the number of gates required remains the same, but the number of terms in the linear combination of the Suzuki-Trotter decomposed time-evolution operators is O(rn). Therefore, for example, to evaluate the expectation value ofĤ n ST(r) (∆ τ ) with respect to a given state |ψ , the O(rn) state overlaps such as ψ|[Ŝ (p) 2m (∆ τ /2)] n−2k |ψ have to be estimated. However, these quantities can be evaluated on quantum computers separately in parallel. Considering the gate count that scales as O(Nn) for approximating the Hamiltonian powerĤ n , the quantum power method is a potentially promising application for nearterm quantum devices.
In contrast to the quantum power method, the direct evaluation of the expectation value ψ|Ĥ n |ψ requires the average of O(N n ) operators, possibly containing long strings of Pauli operators, provided thatĤ consists of O(N) local terms. Although the depth of the circuits for these averages is O(1), the number O(N n ) of averages required makes the direct evaluation of ψ|Ĥ n |ψ unfeasible as soon as the power n is large.
In classical computation, the computational complexity scales as O(N D n) for the evaluation ofĤ n |ψ , when the HamiltonianĤ is local and thus the Hamiltonian matrix is sparse. Here, N D is the dimension of the Hilbert space, i.e., N D = 2 N for the spin-1/2 Heisenberg model given in Eq. (1). This implies that the computational complexity of the classical computation scales exponentially in N. Therefore, the quantum power method introduced here, though it is approximate, would have a quantum advantage over the classical counterpart of the power method.
Finally, we should note that the form of the approximated Hamiltonian powerĤ n ST (∆ τ ) in Eq. (37) suggests a direct treatment of the linear combination of the Suzuki-Trotter decomposed time-evolution operators with a single quantum circuit [74][75][76] forming a simple recursive structure. Figure 2 shows such a circuit structure for probabilistically generating the state ∝ [Ŝ (p) 2m (∆ τ /2) −Ŝ (p) 2m (−∆ τ /2)] n |ψ , among 2 n superposed states, in the N register qubits using n ancilla qubits. However, finding the desired state in the register qubits becomes exponentially small in general if n is large. Let us define P b 1 b 2 ···b n as the probability for finding a bit string b 1 b 2 · · · b n by measuring the n ancilla qubits (b k = 0 or 1 for 1 k n). Then the probability for finding the bit string 11 · · · 1, which is relevant forĤ 2n ST (∆ τ ) [77], is given by If |ψ were an eigenstate ofŜ (p) 2m (∆ τ /2) with an eigenvalue e iλ(∆ τ ) , it oscillates as P 11···1 = [sin λ(∆ τ )] 2n , but otherwise is exponentially small. Therefore, the linear combination of the Suzuki-Trotter decomposed time-evolution operators is better treated with classical computers in the form of Eq. (36).

IV. KRYLOV-SUBSPACE DIAGONALIZATION
As an application of the quantum power method, here we consider the Krylov-subspace diagonalization. We first define a block Krylov subspace and review the subspacediagonalization scheme [78]. We then describe how the quantum power method is combined with the Krylov-subspace diagonalization. Other applications of the quantum power method are outlined in Appendix C and Appendix D.

A. Block Krylov subspace
The block Krylov subspace of the HamiltonianĤ with reference states {|q k } M B k=1 is given as where we call M B 1 the block size. We should note that the reference states {|q k } M B k=1 do not have to be orthogonal to each other but they are linearly independent.
k=1 reduces to the conventional Krylov subspace. By defining with i = k + (l − 1)M B and l = 1, 2, · · · , n, the block Krylov subspace can be written simply as

B. Rayleigh-Ritz technique
Suppose that the ground state |Ψ 0 of the HamiltonianĤ, satisfyingĤ with E 0 being the ground-state energy, should be approximated with the (non-orthonormal) basis states where {v i } nM B i=1 are the expansion coefficients to be determined. The expansion coefficients {v i } nM B i=1 can be determined by minimizing the energy expectation value Ψ KS |Ĥ|Ψ KS under the constraint Ψ KS |Ψ KS = 1. To this end, let us define the following function: where is a Lagrange multiplier, is the subspace Hamiltonian matrix, and is the subspace overlap matrix. Then the condition ∂F /∂v * i = 0 for 1 i nM B yields a generalized eigenvalue problem Since both H and S are Hermitian, the condition ∂F /∂v i = 0 for 1 i nM B yields the same equation. The lowest eigenvalue and the corresponding eigenvector v in Eq. (51) provide an approximation to the ground-state energy E 0 and the expansion coefficients {v i } nM B i=1 in Eq. (47), respectively. Note that when M B = 1, the matrices H and S correspond to the Hankel matrices M n−1 and L n−1 , respectively, defined in Eqs. (D6) and (D5).
Since S is a Hermitian matrix, it can be diagonalized by a unitary matrix V as where s is the diagonal matrix that contains the eigenvalues of S. Note that s > 0 because S is a Gram matrix and hence is positive definite. By using a matrix Eq. (51) can be transformed to a standard Hermitian eigenvalue problem of the form where and q = W −1 v. Thus, by solving the eigenvalue problem of Eq. (54), one can obtain and v = W q. The eigenvector v with the lowest eigenvalue provides the coefficients in the approximate ground state |Ψ KS [see Eq. (47)] with its energy E KS of the HamiltonianĤ in the Krylov subspace K n Ĥ , {|q k } M B k=1 . We note that if we use the Cholesky decomposition S = R † R with R being an upper-triangular matrix, instead of the eigen decomposition in Eq. (52), T reduces to the tridiagonal matrix in the Lanczos method when M B = 1 [78].
Considering the Rayleigh-Ritz technique in a quantumclassical-hybrid computation, it is suited for quantum hardware to evaluate the matrix elements of H in Eq. (49) and S in Eq. (50), because the states {|u i } nM B i=1 are defined on the Hilbert space of N D = 2 N dimensions. On the other hand, the eigenvalue problem in the nM B -dimensional block Krylov subspace given in Eq. (51) or Eq. (54) can be solved on classical computers, assuming that the Krylov subspace approximates reasonably well the eigenspace of the ground state with relatively small n and M B , despite that the dimension N D of the full Hilbert space could be much larger than nM B . This feature is shared with other quantum-classical-hybrid subspacediagonalization schemes reported previously [52][53][54][55][56][57].
We can now approximate the Hamiltonian powerĤ l−1 appearing in the Krylov-subspace basis |u i given in Eq. (45) as where with i = k + (l − 1)M B and l = 1, 2, · · · , n. Note that the systematic errors in Eq. (56) are absent when l = 1. As described in Sec. III D, to approximate the Hamiltonian power H l−1 byĤ l−1 ST(r) (∆ τ ) as in Eq. (57), the Suzuki-Trotter decomposed time-evolution operatorsŜ (p) 2m (±∆ τ /2) have to be applied at most l − 1 times to a state |q k . This implies that the circuit depth required for constructing the block Krylov subspace The circuit depth does not depend on the order r of the Richardson extrapolation.
However, for the purpose of solving the generalized eigenvalue problem in Eq. (51) or the corresponding standard eigen value problem in Eq. (54), we can evaluate the matrix elements in Eqs. (49) and (50) more directly as and with i = k + (l − 1)M B and j = k + (l − 1)M B for 1 k, k M B and 1 l, l n in the block Krylov subspace k=1 . To be more specific, the matrix elements ofH andS for r = 0, i.e., without the Richardson extrapolation, are given as The number of state overlaps required for constructing all matrix elements of bothH andS is thus O(nM 2 B ). If the rth order Richardson extrapolation is employed, the number of state overlaps to be evaluated is increased by a factor of (r + 1). The state overlaps in Eqs. (62) and (63) can be evaluated with an Hadamard-test like circuit, for example [47,[79][80][81].

V. RESULTS
In this section, we use numerical simulations to first show how the quantum power method can control the systematic errors in approximating the Hamiltonian powerĤ n and then demonstrate the application of the quantum power method combined with the multireference Krylov-subspace diagonalization for the spin-1/2 Heisenberg model.

A. Degree of approximation
We first examine quantitatively how the Hamiltonian power H n is approximated byĤ n ST(r) (∆ τ ). For this purpose, we define a distance d(Â,B) between operatorsÂ andB as where Â ,B F denotes the Frobenius inner product betweenÂ andB defined by and ||Â|| F denotes the Frobenius norm ofÂ, i.e., Note that Â , Evaluating the distance is costly as it demands matrixmatrix multiplications or diagonalizations. To avoid such costly operations, we employ a stochastic evaluation of the trace as [82][83][84][85][86] whereX ∈ Â †Â ,B †B ,Â †B and is a random-phase state with {|x } being a complete orthonormal basis set such that x|x = δ xx and φ ζ (x) being a random variable drawn uniformly from [0, 2π). Note that φ ζ |φ ζ = 2 N , i.e., the dimension N D of the Hilbert space. We choose {|x } as the orthonormal basis set that diagonalizes the local Pauli Z operators. The stochastic evaluation of the trace in Eq. (67) requires only sparse matrix-vector multiplications and a single inner-product calculation for each ζ, ifX is represented as a product of sparse matrices, which is indeed the case here. Instead of taking the limit R → ∞, we fix R = 16 for N 12 and R = 256 for N = 10 and estimate error bars.
Since φ ζ |Â †Â |φ ζ , φ ζ |B †B |φ ζ , and φ ζ |Â †B |φ ζ forÂ =Ĥ n andB =Ĥ n ST(r) (∆ τ ) are highly correlated to each other, error bars of d(Â,B) must be estimated using the corresponding 3 × 3 covariance matrix. Figure 3 shows the distance as a function of ∆ τ for n = 1, 2, . For each n, the distance with the Richardson extrapolation is an order of magnitude smaller than that without the Richardson extrapolation. The leading systematic error inĤ n , and the distance indeed scales almost linearly in ∆ 4 τ . As expected from Eq. (40), essentially no difference can be found between the results withŜ 2 andŜ (3) 4 , indicated respectively by empty and filled symbols in Fig. 3. These results clearly demonstrate that the systematic errors in approximating the Hamiltonian power H n are well controlled. Figure 4(a) shows the n dependence of the distance for N = 24 with various values of ∆ τ calculated using the lowestorder symmetric Suzuki-Trotter decompositionŜ 2 . The distance first increases with n and tends to saturate at n ∼ 100. It is remarkable to find in Fig. 4(b) that, even with the large power exponents as large as n = 100, the linear dependence of the distance on ∆ 2 τ remains in a wide range of ∆ τ (∆ τ J 0.1) and the distance is smoothly extrapolated to zero in the limit of ∆ τ → 0, clearly demonstrating the controlled accuracy of the quantum power method. Figures 4(c) and 4(d) show the same results but obtained by using the first-order Richardson extrapolation (r = 1), for which the systematic errors in approximating the Hamiltonian powerĤ n are expected to be O(∆ 4 τ ). Indeed, our numerical simulations find the linear dependence of distance on ∆ 4 τ for at least ∆ τ J 0.05 when n = 100 [see the inset in Fig. 4(d)]. Notice also that the distance itself becomes smaller by the factor of approximately 5 even for large n when the first-order Richardson extrapolation is employed.

B. Ground-state energy and fidelity
We now perform numerical simulations of the Krylovsubspace diagonalization combined with the quantum power method to calculate the ground-state energy and fidelity of the spin spin-1/2 Heisenberg model described by the Hamiltonian H in Eq. (1) on a periodic chain of N = 16 sites (i.e, qubits).
Considering the Krylov-subspace diagonalization as an application of the quantum power method on near-term quantum computers, it is crucial to reduce the circuit depth. As discussed in Sec. III D and Sec. IV C, the depth of the circuit required for constructing the block Krylov subspace k=1 scales as O(n) with a prefactor D (p) 2m . Since m and p in the symmetric Suzuki-Trotter decomposi-tionŜ (p) 2m can be set to the minimum values m = 1 and p = 3, at least for the system sizes examined in the previous section including N = 16, the primary objective here is to reduce the power n. For this purpose, we first describe the selection of the reference states, aiming that the block Krylov subspace K n Ĥ ST(r) (∆ τ ), {|q k } M B k=1 spanned by these reference states can approximate reasonably well the target subspace, which in the present case is the eigenspace of the ground state ofĤ. Then we show by numerical simulations how the selection of the reference states affects the convergence to the ground state with n.

Selection of reference states
Equation (47) suggests that the ground state |Ψ 0 can be well approximated if the reference states {|q k } M B k=1 are chosen so that these states have substantial overlap with the exact ground state. Therefore, as the reference states, we introduce the following product states for the subspace diagonalization: where is the spin-singlet state which is an eigenstate of the swap operatorP i j with eigenvalue −1 and is also known as one of the Bell states, are the eigenstates ofŶ i with eigenvalues ±1, and |0 i and |1 i are the eigenstates ofẐ i with eigenvalues ±1. |Φ A and |Φ B are the ground states ofĤ A andĤ B , respectively, while others are the Néel states that are the ground states when a mean-field theory is applied to the Hamiltonian. These product states are expected to have a sizable overlap with the exact ground state (also see Fig. 6) and, moreover, are easy to be prepared from |0 ⊗N with appropriate combinations of Pauli, Hadamard, phase, and CNOT gates.
Another relevant candidate might be a variational state that has a substantial overlap with the ground state. We thus introduce as another reference state, where |Ψ VQE is an approximate ground state prepared with a VQE scheme. Specifically, we choose |Ψ VQE as a resonating-valence-bond-type wave function without the symmetry projection operator, containing 64 optimized variational parameters for N = 16 that do not reflect the spatial symmetry of the Hamiltonian, as reported in Ref. [87]. While the exact ground-state energy is E 0 /N J = −0.196393522, our variational state |Ψ VQE has the variational energy Ψ VQE |Ĥ|Ψ VQE /N J = −0.1885 (also see Fig. 5) and the ground-state fidelity | Ψ 0 |Ψ VQE | 2 = 0.771 (also see Fig. 6).
In our previous study [87], we have shown that restoration of the spatial symmetry that is broken by a circuit ansatz greatly improves the ground-state-energy estimation as well as the ground-state fidelity. Motivated by this finding, we introduce another set of the reference states {|q k } N k=1 with whereT k is a unitary operator representing the onedimensional k-lattice-space translation withT 0 =Î, and |Ψ VQE is the same state given in Eq. (77). With this set of the reference states, the translational symmetry that is broken in the apparent circuit structure of |Ψ VQE can be restored as a linear combination of the states in the block Krylov subspace, without applying a projection operator to |Ψ VQE . For example, a simple sum of these N reference states {|q k } N k=1 , i.e., N k=1 |q k , is translationally symmetric with momentum zero. The reference states |Φ A , |Φ B , |Ψ VQE , and {|q k } N k=1 introduced above are all spin-singlet states, i.e, the total spin and the Z-component of the total spin being zero, while the X-, Y-, and Z-components of the total spin are zero for the reference states |X AFM1(2) , |Y AFM1(2) , and |Z AFM1(2) , respectively. Because the HamiltonianĤ considered here is spin SU(2) symmetric and the quantum power method preserves the Hamiltonian symmetry as shown in Eq. (42), the Krylov subspace generated from these reference states remains in the same symmetry sector of the Hilbert space as the reference states. We select these reference states because it is known that the ground state of the spin-1/2 Heisenberg model considered here is spin singlet [88]. Let us first focus on the results for n = 1, where no Hamiltonian power is incorporated in the Krylov subspace. It is not surprising to find that the energy and the fidelity are substantially improved if the reference states include the VQE state |Ψ VQE . The improvement is even more significant if we incorporate the spatially translated VQE states {|q k } N k=1 . Note that, if M B = 1, the energies indicated at n = 1 are merely the expectation values ofĤ with respect to the corresponding reference state, e.g., Φ A |Ĥ|Φ A /N J = −0.125 and Ψ VQE |Ĥ|Ψ VQE /N J = −0.1885. The multireference scheme with M B > 1 further decreases the energy and improves the fidelity without applying the Hamiltonian power to the reference states.

Ground-state energy and fidelity
With increasing the power n, the energy decreases monotonically and the fidelity keeps increasing towards one, implying that the ground state estimation can be improved system- atically over a chosen set of reference states without any parameter optimization. The nearly linear behavior of E KS − E 0 in the semilog plot shown in Fig. 5(b) suggests the exponential convergence to the exact ground-state energy as a function of n, as in the Lanczos method [10]. Notice also that the energy as well as the fidelity for M B = 16 is consistently better than those for M B 9 for every n. Moreover, the slope in the semilog plot of E KS − E 0 and also the slop of the fidelity tend to be steeper for M B > 1 than for M B = 1, implying that the convergence towards the ground state is improved more efficiently in the multireference scheme with M B > 1. Interestingly, even if |Ψ VQE is not included in a set of reference states, the multireference schemes with M B = 2 and M B = 8 surpass the scheme including only |Ψ VQE with M B = 1 at n = 5 and 3, respectively, in terms of the ground-sate energy E KS . Therefore, the multireference scheme with M B > 1 works effectively for reducing the power n and hence the number of gates in a circuit, even if simple product states with no variational parameters are chosen for the reference states. Table I summarizes the minimum dimension n per block size of the Krylov subspace and the corresponding circuit depth required for converging the ground-state energy E KS with an accuracy (E KS − E 0 )/N J 10 −4 for N = 16. Note here that the commuting exponentials in [Ŝ 2 (±∆ τ /2)] n−1 are contracted when the circuit depth is counted.

VI. CONCLUSION AND DISCUSSION
We have proposed the quantum power method that approximates the Hamiltonian powerĤ n with a linear combination of the time-evolution operators. The key ingredients of the quantum power method are the central-finite-difference scheme for the time derivatives and the symmetric Suzuki-Trotter decomposition that decomposes each time-evolution operator, both of which guarantee that the approximated Hamiltonian powerĤ n ST (∆ τ ) retains the Hermiticity and the even parity in ∆ τ , i.e.,Ĥ n ST (∆ τ ) = Ĥn ST (∆ τ ) † =Ĥ n ST (−∆ τ ), with the controlled accuracy of the finite-difference error E FD ∼ O(∆ 2 τ ) and the Suzuki-Trotter error E ST ∼ O(∆ 2m τ ). The number of gates required for approximating the Hamiltonian powerĤ n is O(Nn), where N is the number of qubits and the Hamilto-nianĤ is assumed to be local. This should be contrasted to the classical power method which scales exponentially in N.
The rth order Richardson extrapolation can be adopted to systematically improve the systematic errors as E FD ∼ O(∆ 2+2r τ ) and E ST ∼ O(∆ 2m+2r τ ) in the approximated Hamiltonian powerĤ n ST(r) (∆ τ ), without increasing the number of gates required in each quantum circuit, although the number of terms in the linear combination, which can be treated classically, increases by the factor r + 1. Thus, both with and without the Richardson extrapolation, the systematic errors E FD and E ST can be consistently treated with the lowestorder Suzuki-Trotter decomposition with m = 1, independently of the power n, which reduces significantly the circuit depth as compared with the algorithm that requires the higherorder Suzuki-Trotter decomposition with increasing the power n (see Appendix B). Therefore, the quantum power method proposed here is potentially promising for near-term quantum devices.
By numerical simulations, we have tested the quantum power method and found that the Hamiltonian powerĤ n for the spin-1/2 Heisenberg model can be well approximated bŷ H n ST(r) (∆ τ ) with the controlled accuracy to be essentially exact for the power n up to 100 and N as large as 24 qubits, corresponding to the Hilbert space dimension N D = 2 N ≈ 10 7 .
As an application of the quantum power method, we have demonstrated, with noiseless numerical simulations, the multireference Krylov-subspace diagonalization combined with the quantum power method for the spin-1/2 Heisenberg model on an N = 16 qubit ring to evaluate the ground-state energy and the ground-state fidelity. Considering the Hamiltonian powerĤ n up to n = 11, we have shown that the multireference Krylov-subspace diagonalization scheme with the block size M B > 1 greatly accelerate the convergence to the ground state, even with simple parameter-free product states for the reference states. We have also found that the Krylov-subspace diagonalization scheme with M B = 1, corresponding to a quantum version of the standard Lanczos method [78], improves the ground-state energy of the VQE state |Ψ VQE almost exponentially with increasing n. Thus, the Krylov-subspace diagonalization combined with the quantum power method, which satisfies the variational principle by definition, can provide a systematic way to further improve a VQE state that has already a reasonable overlap with an exact ground state. This is a quantum analog to the Lanczos iteration scheme in the variational Monte Carlo method on classical computers [17], but here one can treat higher powers of the Hamiltonian on quantum computers though approximately yet with a controlled accuracy.
The quantum power method proposed here can be easily generalized to higher spatial dimensions or a more complicated system described by a Hamiltonian that is decomposed intoĤ =Ĥ A +Ĥ B +Ĥ C + · · · , where generally [Ĥ Γ ,Ĥ Γ ] 0 if Γ Γ but terms within eachĤ Γ commute to each other (here, Γ, Γ = A, B, C, · · · ). Even in this case, the number of gates required in approximating the Hamiltonian powerĤ n scales similarly, except for the additional prefactor of O(N Γ ), where N Γ is the number of sub-divided Hamiltonians inĤ. This is because when the lowest-order symmetric Suzuki-Trotter decomposition is employed, e.g., , the time-evolution operator e −i∆ τĤ is approximated generally as a product of 2N Γ − 1 time-evolution operators of the subdivided systems. In the multireference Krylov-subspace diagonalization, one can choose, for example, simply the ground states |Φ A , |Φ B , |Φ C , · · · of the sub-divided Hamiltonianŝ H A ,Ĥ B ,Ĥ C , · · · for the reference states. The mean-field ground states can also be used for the reference states.
Although we have simulated only the ground-sate energy, the expectation value of other observables that commute with the HamiltonianĤ can be evaluated similarly. When an ob-servableÔ does not commute with the HamiltonianĤ, the expectation value with respect to the approximate ground state , can also be evaluated as where |ũ i =Ĥ l−1 ST(r) (∆ τ )|q k , as given in Eq. (57), and the explicit form ofĤ l−1 ST(r) (∆ τ ) with r = 0 is used in the second line.
Our numerical simulations clearly demonstrate a promising potential that the quantum power method combined with the multireference Krylov-subspace diagonalization enables us to perform systematic and optimization-free calculations for quantum many-body systems, which is suitable to nearterm quantum computers. Other applications of the quantum power method include various moment based methods, which are briefly outlined in Appendix C and Appendix D. In these appendixes, we show that the power method can evaluate Ψ|Ĥ 2 |Ψ with exactly the same amount of resource that is required for Ψ|Ĥ|Ψ , and therefore, for example, the energy variance σ 2 = Ψ|Ĥ 2 |Ψ − Ψ|Ĥ|Ψ 2 can be easily obtained. Here, |Ψ is a given quantum state. Using numerical simulations, we also demonstrate the CMX for the imaginarytime evolution. This formalism can be easily extended to other methods, e.g., the high-temperature series expansion [89].
Finally, we remark that the quantum power method proposed here can generally be applied to any sparse Hermitian operatorÂ. In this case, the nth power ofÂ is given aŝ In this Appendix, we provide a Python program that generates coefficients required for the higher-order symmetric Suzuki-Trotter decompositionsŜ (p) 2m (∆ τ ) introduced in Sec. III C 2, and examine numerically the systematic errors due to the Suzuki-Trotter decompositionsŜ (p) 2m (∆ τ ) with different parameters m and p. Note that m is an integer with m 1 and p is an odd integer with p 3

Coefficients for higher-order Suzuki-Trotter decompositions
Listing 1 shows a Python program that generates the coef- i=1 for a given set of parameters m and p in the symmetric Suzuki-Trotter decompositionsŜ (p) 2m (∆ τ ): where x = −i∆ τ and D (p) 2m = 2p m−1 + 1 derived in Eq. (34). The program includes an example for m = 2 and p = 5. In this case, the symmetric Suzuki-Trotter decomposition has a form S (5) 4 (∆ τ ) = e xs 1ĤA e xs 2ĤB e xs 3ĤA e xs 4ĤB × e xs 5ĤA e xs 6ĤB e xs 7ĤA × e xs 8ĤB e xs 9ĤA e xs 10ĤB e xs 11ĤA (A2) and the output of the program gives the 11 coefficients for other values of m and p.

Notice that the coefficients {s
i=1 are symmetric, i.e., and satisfy the following sum rule: for any m and p. Although it is sufficient to find the coeffi- for our purpose, the program can also output a cumulative sum T i of the coefficient s i defined as By plotting T i as a function of i (or i/d (p) 2m ) for several values of m with a fixed p, one can find a fractal feature appearing in the higher-order Suzuki-Trotter decompositions [58,72].

Listing 1. A Python program for generating the coefficients {s
in the symmetric Suzuki-Trotter decompositionŜ (p) 2m .

Numerical examination of a Suzuki-Trotter error
Here we numerically examine the systematic errors due to the Suzuki-Trotter decompositionsŜ (p) 2m (∆ τ ) with different parameters m and p. The Trotter formula [90][91][92] where M is an integer such that t = M∆ τ . Figure 7 shows the real part of the difference between the exact propagator and the approximated propagator i.e., with ∆ τ J = 0.07 and 0.1 for the spin-1/2 Heisenberg model on an N = 16 qubit ring described by the HamiltonianĤ in Eq. (1). Here, |Ψ 0 is the exact ground state. The exact propagator is simply given by K(t) = e −iE 0 t , where E 0 is the exact ground-state energy. As expected, when p is fixed, the error decreases by orders of magnitude with increasing m. It is also found that, when m is fixed, the error decreases by orders of magnitude with increasing p. Although we only show ReδK(t), the imaginary part of the difference, ImδK(t), behaves similarly. We should emphasize here that while the deviation of the approximated propagatorK(t) from the exact one K(t) becomes larger in the long time limit (tJ 1), the quantum power method proposed here is formulated on the basis of the time-evolution operatorsÛ(t) at time t close to zero, for which the deviation is small. Therefore, this is another advantage of the quantum power method in controlling the Suzuki-Trotter error over other quantum algorithms that require the long-time dynamics approximately described by the Suzuki-Trotter decomposed time-evolution operators. It is also found that the error decreases with increasing p for a fixed m, independently of ∆ τ . This suggests that the increase of p reduces the coefficient of the leading-error term by orders of magnitude. However, as shown in Eq. (34), p is the base of the exponential which determines the circuit depth. Thus, as far as noisy near-term quantum computers are concerned, p = 3 might be a more suitable value than p 5.
Finally, we note that several exponential-product formulas, not limited to those found by Suzuki, up to the depth 11 with an error analysis can be found in Ref. [93]. Other error analysis of the Suzuki-Trotter decomposition devoted for quantum computing can be found in Refs. [94][95][96].
Appendix B: Another formalism for approximating the Hamiltonian power As discussed in Sec. III D, we can formulate at least two different algorithms for evaluating the Hamiltonian powerĤ n , depending on in which stage the time-evolution operators in the approximated Hamiltonian powerĤ n (∆ τ ) is replaced with the symmetric Suzuki-Trotter decomposition, either in Eq. (10) or in Eq. (14). In the quantum power method described in Sec. III D, the time-evolution operators in Eq. (14) are approximated by the symmetric Suzuki-Trotter decomposition. In this Appendix, we describe the other formalism by approximating the time-evolution operators in Eq. (10) and show that the resulting algorithm scales differently from the one formulated in Sec. III D.
By incorporating the symmetric Suzuki-Trotter decomposi-tionŜ (p) 2m into the approximated Hamiltonian powerĤ n (∆ τ ) in Eq. (10), the Hamiltonian powerĤ n is now approximated aŝ O(∆ 2 τ ) represents the systematic error E FD due to the finitedifference scheme for the time derivatives, and E ST denotes the systematic error due to the Suzuki-Trotter decomposition of the time-evolution operators, the order of E ST being discussed bellow. We should emphasize here thatĤ n for ∆ τ −∆ τ , although the exact time-evolution operators satisfy the multiplication lawÛ(∆ τ )Û(∆ τ ) =Û(∆ τ +∆ τ ). Note also thatĤ 1 ST (∆ τ ) =Ĥ 1 ST (∆ τ ). We can readily confirm thatĤ n ST (∆ τ ) is Hermitian and an even function of ∆ τ , i.e., as in the case ofĤ n ST (∆ τ ) given in Eq. (39) and hence the systematic error E ST (as well as the systematic error E FD , see Sec. III A) in odd powers of ∆ τ is absent in Eq. (B1). We can also show thatĤ n ST (∆ τ ) does not satisfy the law of exponents, i.e.,Ĥ for n 2, simply because of Eq. (B3), but only satisfies it approximately within the systematic errors. This is in sharp contrast to the case ofĤ n ST (∆ τ ), which satisfies exactly the law of exponents in Eq. (38).
At first glance, one would tend to conclude thatĤ n ST (∆ τ ) in Eq. (B2) is more suitable to approximate the Hamiltonian powerĤ n thanĤ n ST (∆ τ ) in Eq. (36), because each team in H n ST (∆ τ ) contains a singleŜ (p) 2m , not a product of multipleŜ (p) 2m 's as inĤ n ST (∆ τ ), thus expecting the less number of gates in the circuit. However, the disadvantage ofĤ n ST (∆ τ ) in Eq. (B2) is that the higher-order Suzuki-Trotter decompositions are required for approximating the Hamiltonian powerĤ n with larger n.
This can be understood by recalling thatŜ (p) 2m (t) has a form of Eq. (29): Accordingly, the higher-order derivative ofŜ (p) 2m (t) at t = 0 is given by for n 2m but for n > 2m. For example, if n = 2m + 1, the derivative reads It is now important to notice that the right-hand side of Eq. (B2) corresponds to the central finite-difference approx- , i.e., In other words, the approximated Hamiltonian powerĤ n ST (∆ τ ) in Eq. (B2) is given by the higher-order derivative ofŜ (p) 2m (t) at t = 0 asĤ It is now obvious that the formalism in Eq. (B2) breaks down if n > 2m because in this case, according to Eq. (B8), Ĥn , which contradicts to Eq. (B1). Therefore, 2m n (B12) is required for approximating the Hamiltonian powerĤ n bŷ H n ST (∆ τ ) under a controlled accuracy with the systematic error This is the most important difference from the algorithm described in Sec. III D, where the lowest-order Suzuki-Trotter decomposition with m = 1 is adequate for any power n.
There are two remarks in order. First, the approximated Hamiltonian powerĤ n ST (∆ τ ) in Eq. (37) can be considered aŝ Therefore, the lowest-order symmetric Suzuki-Trotter decomposition with m = 1 is adequate to satisfy Eq. (B12) and indeed, as discussed in Sec. III D, it approximates the Hamiltonian powerĤ n with the controlled accuracy. Second, al- 4 (m = 2, p = 7, D = 15) 6 (m = 3, p = 5, D = 51) though we have emphasized that the violation of the multiplication lawŜ (p) is the essential point that distinguishes the two algorithms described here and in Sec. III D, this equation is satisfied within the systematic error. i.e., Accordingly, the two algorithms described here and in Sec. III D should be the same within the systematic error. In fact, the approximated Hamiltonian powersĤ n ST (∆ τ ) and H n ST (∆ τ ) in Eqs. (36) and (B2), respectively, are equivalent within the systematic error becausê provided that 2m + 1 > n, which is consistent with Eq. (B12).
As an example, we show in Fig. 9 the expectation values Ĥn ST (∆ τ ) and Ĥn ST (∆ τ ) with respect to the quantum states |Φ A and |Ψ VQE of the spin-1/2 Heisenberg model on an N = 16 qubit ring for n = 3. Here, a simplified notation of the expectation value · · · ≡ Ψ| · · · |Ψ (B19) with |Ψ ∈ {|Φ A , |Ψ VQE } is introduced. According to Eqs. (B7)-(B9), Ĥ3 ST (∆ τ ) in the limit of ∆ τ → 0 should converge as for m = 1. Here, the explicit form of the residual termR 3 in Eq. (B21) can be derived by using the Baker-Campbell-Hausdorff formula forŜ 2 as [59,93,97] The numerical results in Fig. 9   for m = 1. Note also that the linear convergence of these quantities to the exact values as a function of ∆ 2 τ shown in Fig. 9 corroborates the systematic errors expected for Ĥn ST (∆ τ ) in Eqs. (B10) and Ĥn ST (∆ τ ) in Eq. (35). Now we discuss the gate count for approximatingĤ n witĥ H n ST (∆ τ ). As described above, Eq. (B12) sets the order of the Suzuki-Trotter decomposition such that 2m n, i.e., the smallest order m of the Suzuki-Trotter decomposition to eval-uateĤ n being m = n/2 , where · is the ceiling function that returns the minimum integer larger than or equal to the argument. Therefore, assuming that the HamiltonianĤ is local, the number of gates required for approximatingĤ n witĥ H n ST (∆) is O(N p n/2 ) because the circuit depth D (p) 2m for the single Suzuki-Trotter decomposed time-evolution operatorŜ (p) 2m is given by Eq. (33), and thus increases exponentially in the power n. In contrast, as described in Sec. III D, the number of gates required for approximatingĤ n withĤ n ST (∆ τ ) is O(Nn) with a prefactor D (p) 2m ∼ O(1), i.e., increasing polynomially in N and n.
To apply the quantum power method formulated in this appendix to the Krylov-subspace diagonalization scheme, it is crucial to reduce the maximum power n appearing in the formalism. By defining for the basis set generated in the block Krylov subspace , the matrix elements H and S in Eqs. (49) and (50) are now approximated by replacing |u i with |ũ i as where i = k + (l − 1)M B and j = k + (l − 1)M B for 1 k, k M B and 1 l, l n. As compared with Eqs. (60) and (61), the power exponents are now distributed to the left and the right basis states.
To be more specific,H i j andS i j in terms ofŜ (p) 2m are given asH where t (l−1) Finally, we note that the systematic errors E FD and E ST in Eq. (B1) can be improved systematically, without increasing the gate count of each circuit, by adopting the Richardson extrapolation asĤ whereĤ n ST(r) (∆ τ ) is the rth order Richardson extrapolation of the approximate Hamiltonian power, i.e., is a linear combination of n + 1 unitariesŜ (p) 2m ,Ĥ n ST(r) (∆ τ ) is a linear combination of (r + 1)(n + 1) unitariesŜ (p) 2m .

Appendix C: Moment methods
In this appendix, we outline moment methods as other applications of the quantum power method to evaluate the moments and cumulants of the Hamiltonian. By using numerical simulations, we demonstrate the connected-cluster expansion (CMX) for a short-time imaginary-time evolution and estimate the ground-state energy of the spin-1/2 Heisenberg model. These numerical results are compared with those obtained by the multireference Krylov-subspace diagonalization combined with the quantum power method discussed in Sec. V B 2. We also show that the quantum power method can particularly simply evaluate the lowest order moments.

Moment and cumulant
The Feynman propagator with respect to a state |Ψ can be written as whereÛ(t) is the time-evolution operator given in Eq. (7) and is the nth Hamiltonian moment. We also define the generating function Φ(t) of the cumulants {κ n } as Thus, the nth moment µ n and cumulant κ n are given by the nth time derivative of generating functions K(t) and Φ(t), respectively, as and We note that, recently, a method making use of the expectation value of the time-evolution operator has been proposed for evaluating eigenvalues of the Hamiltonian [98]. It should be noticed [99] that the nth moment µ n can be expressed as and, equivalently, the nth cumulant κ n can be expressed as Therefore, from the moments {µ k } k n , one can obtain the cumulants {κ k } k n and vice versa. A remarkable difference between these two quantities is that the magnitude of the moment grows exponentially in n as µ n ∼ O(N n ), while the magnitude of the cumulant remains as κ n ∼ O(N) [100].
We can apply the same argument for the cumulants. Using the central-finite-difference method, the nth cumulant is evaluated as Because Φ(−∆ τ ) = Φ(∆ τ ) * , κ n (∆ τ ) for n odd and n even can be expressed, respectively, as and where Φ(0) = ln K(0) = 0 is used in Eq. (C15). Note that, if we write the propagator as K(t) = a(t)e iϕ(t) with a(t) and ϕ(t) real, then ReΦ(t) = ln a(t) and ImΦ(t) = ϕ(t), implying that the cumulants with odd order are related to the phase of K(t), while the cumulants with even order are related to the amplitude of K(t). Recently, an efficient method for estimating the overlap amplitude of two pure states has been proposed [101]. Such a method might be utilized for evaluating the cumulants with even order.

Quantum power method for moment and cumulant
As shown explicitly in the previous section, the nth moment µ n can be approximated as a linear combination of the Feynman propagator K(t), i.e., the expectation value of the timeevolution operatorÛ(t), evaluated at different time variables t (n) i = n 2 − i ∆ τ for i = 0, 1, . . . , n. Similarly, the nth cumulant κ n can be approximated as a linear combination of Φ(t), i.e., logarithm of the Feynman propagator K(t), evaluated at different time variables t (n) i . Therefore, an important quantity here is again the time-evolution operatorÛ(t).
To implement on quantum computers, the time-evolution operator is further decomposed approximately by using the symmetric Suzuki-Trotter decomposition as in Eq. (27). However, at this point, it is crucially important to recall the argument given in Sec. III D and Appendix B. Although the time evolution operatorÛ(t) evaluated at time t (n) i = n 2 − i ∆ τ satisfies thatÛ and thus the approximated nth moment µ n (∆ τ ) in Eq. (C9) is equivalent to these are no longer generally correct when the time-evolution operators are approximated by the Suzuki-Trotter decomposition, i.e.,Ŝ (p) 2m Therefore, the Feynman propagator K(t (n) i ) in Eq. (C9) can be approximated either as If the Feynman propagator K(t (n) i ) is approximated as in Eq. (C20), the nth moment µ n is given by and thus the lowest-order symmetric Suzuki-Trotter decom-positionŜ (p) 2m with m = 1 can be adopted (see Sec. III D). This approach is suitable for the calculations of higher order moments and cumulants. On the other hand, if the Feynman propagator K(t (n) i ) is approximated as in Eq. (C19), the higher-order symmetric Suzuki-Trotter decompositionŜ (p) 2m is required. As discussed in Appendix B, in order to evaluate the nth moment µ n with the controlled accuracy, the order of the symmetric Suzuki-Trotter decompositionŜ (p) 2m must be 2m n [see Eq. (B12)]. In this case, the systematic error is O(∆ 2 τ ), i.e., Therefore, this approach is not suitable for large n but is more preferable than the other approach when n 4. The same argument is applied for the cumulant κ n .

First and second moments
The first and second moments are the most fundamental quantities for many practical purposes because µ 1 = Ĥ is the average of the energy and µ 2 = Ĥ2 is related to the variance of the energy. The first moment Ĥ is directly evaluated by measuring each term of the HamiltonianĤ on quantum computers. Perhaps, Ĥ2 could also be evaluated in the same way, although terms to be measured are increased by a factor of O(N), assuming that the HamiltonianĤ is local. The quantum power method can provide an alternative approach to evaluate these quantities with the same amount of resource.
From Eqs. (C10) and (C11), we can approximate the first and second moments µ 1 and µ 2 as 10. Quantum circuit to evaluate Re Ψ|Ŝ 2 (∆ τ )|Ψ or Im Ψ|Ŝ 2 (∆ τ )|Ψ . θ in the circuit denotes the phase gates such that θ|0 = |0 andθ|1 = e iθ |1 . Since P 0 − P 1 = Re[e iθ Ψ|Ŝ 2 (∆ τ )|Ψ ], one can evaluate Re Ψ|Ŝ 2 (∆ τ )|Ψ if θ = 0 and Im Ψ|Ŝ 2 (∆ τ )|Ψ if θ = −π/2 from the difference of the probabilities P 0 and P 1 . Here, P b is the probability for finding a bit b (= 0, 1) by measuring the ancilla qubit. and µ 2 (∆ τ ) = 2 respectively. This is already remarkable because the second moment µ 2 is also estimated simply by the expectation value of a single time-evolution operator. To evaluate these quantities on quantum computers, the time-evolution operator U(∆ τ ) is approximated by the lowest-order symmetric Suzuki-Trotter decompositionŜ 2 (∆ τ ) (see Appendix C 3). Therefore, in the quantum power method, the first and second moments µ 1 and µ 2 are estimated simply by evaluating Im Ŝ 2 ∆ τ 2 and Re Ŝ 2 (∆ τ ) , i.e., and respectively. Although we have to introduce an ancilla qubit (see Fig. 10), µ 1 = Ĥ and µ 2 = Ĥ2 can be thus estimated with exactly the same amount of resource. If noise in quantum devices is not destructively serious, this approach based on the quantum power method might be more suitable than the direct approach measuring all terms inĤ andĤ 2 . Figure 11 shows the numerical results of µ 1 and µ 2 evaluated from Eqs. (C25) and (C26) for the spin-1/2 Heisenberg model with two different quantum states. We also show the results obtained by employing the first-order Richardson extrapolation, i.e., for n = 1 and 2, which expects that the systematic error scales as O(∆ 4 τ ), in stead of O(∆ 2 τ ) without the Richardson extrapolation. Our numerical simulations clearly demonstrate that the systematic errors are well controlled and the results converge smoothly to the exact values in the limit of ∆ τ → 0. The quantum power method for the first and second moments could be useful to, e.g., the energy variance minimization for optimizing a parametrized quantum circuit [102]. Here, we only consider the first and second moments, but the higher order moments can be similarly evaluated. For example, the third and fourth moments are given as and respectively. To implement these on quantum computers, U (∆ τ ) is now approximated by using the higher-order Suzuki-Trotter decompositionŜ (p) 4 (∆ τ ) with m = 2 (also see Fig. 9), which is still affordable.
(C32) Figure 12 shows the exact E(τ) and the CMX of the energy truncated at the n max th cumulant E n max (τ) = n max −1 n=0 (−τ) n n! κ n+1 (C33) for the spin-1/2 Heisenberg model, where the VQE state |Ψ VQE in Eq. (77) is selected for the quantum state |Ψ in Eq. (C30). The energy E(τ) with the exact ITE decreases monotonically in τ, because the first derivative of E(τ) is minus of the energy fluctuation [99] dE(τ) dτ where the equality satisfies if and only if |Ψ VQE (τ) is an exact eigenstate (e.g., the ground state) ofĤ. On the other hand, due to the truncation of the series at finite order, E n max (τ) at large τ diverges to −∞ for even n max 2 or to +∞ for odd n max 3. Note that E 2 (τ) = κ 1 − κ 2 τ is the tangent line of E(τ) at τ = 0. We also find that the convergence of E n max (τ) to the exact ground-state energy E 0 with respect to the power exponents n max in the cumulants required is rather slower, as compared with the Krylov-subspace diagonalization with either M B = 1 or M B = 9 discussed in Sec. V B. This is not quite surprising because the form of E(τ) in Eq. (C33) is an expansion around τ = 0, which is analogous to the high-temperature expansion.

Appendix D: Lanczos method
In this appendix, we briefly outline the Lanczos method with an emphasis on its aspect as a moment method [103,104], i.e., a potential application of the quantum power method.
It is instructive to give the explicit forms of the first few matrix elements of T n . The first three matrix elements required for constructing the 2 × 2 matrix T 2 are given by where · · · = q 1 | · · · |q 1 . Therefore, α 1 and β 2 1 are the energy expectation value and the energy variance with respect to the initial state |q 1 , respectively.

Ratio of Hankel determinants
We now describe a way to calculate recursively the ratio of the determinants appearing in Eqs. (D3) and (D4). Let us first review the determinant and the matrix-inversion formulas for general matrices. Let A n be an (n × n) matrix, b be an (n × 1) matrix, c be an (n × 1) matrix, and b be a (1 × 1) matrix, and let us consider an (n + 1) × (n + 1) matrix A n+1 of the form If we define the determinant of A n+1 is given by and the inverse A −1 n+1 is given by Now we apply the above formulas to recursively evaluate the ratios of the determinants of L n and L n−1 . Due to its particular structure, L n can be expressed in terms of L n−1 as which involves the inverse L −1 n−1 whose dimension is less than that of L −1 n by 1. The inverse matrix L −1 n can be calculated using Eq. (D13). Starting with where (L −1 n ) T = L −1 n is used. Thus, starting with the known L −1 0 and using Eqs. (D17) and (D19), one can obtain {r n } recursively as L −1 0 → r 1 → L −1 1 → r 2 → L −1 2 → r 3 → · · · . It should be noted that Eq. (D17) involves a matrix-vector multiplication and, in addition, Eq. (D19) involves a rank-1 update. Therefore, the complexity for computing the ratio of determinants in Eq. (D16) is O(n 2 ). This is more efficient when n is large because the direct calculation of a determinant from scratch, e.g., by using the LU decomposition, requires O(n 3 ) operations. Noticing that [L n ] i j = µ i+ j−2 while [M n ] i j = µ i+ j−1 , the similar recursive formula for M n can be readily derived simply by replacing the indexes for the moments in the above as {µ i } 2n i=0 → {µ i+1 } 2n i=0 .