Perturbation theory with quantum signal processing

Perturbation theory is an important technique for reducing computational cost and providing physical insights in simulating quantum systems with classical computers. Here, we provide a quantum algorithm to obtain perturbative energies on quantum computers. The benefit of using quantum computers is that we can start the perturbation from a Hamiltonian that is classically hard to solve. The proposed algorithm uses quantum signal processing (QSP) to achieve this goal. Along with the perturbation theory, we construct a technique for ground state preparation with detailed computational cost analysis, which can be of independent interest. We also estimate a rough computational cost of the algorithm for simple chemical systems such as water clusters and polyacene molecules. To the best of our knowledge, this is the first of such estimates for practical applications of QSP. Unfortunately, we find that the proposed algorithm, at least in its current form, does not exhibit practical numbers despite of the efficiency of QSP compared to conventional quantum algorithms. However, perturbation theory itself is an attractive direction to explore because of its physical interpretability; it provides us insights about what interaction gives an important contribution to the properties of systems. This is in sharp contrast to the conventional approaches based on the quantum phase estimation algorithm, where we can only obtain values of energy. From this aspect, this work is a first step towards ``explainable'' quantum simulation on fault-tolerant quantum computers.


Introduction
Perturbation theory is one of the most important techniques to understand quantum systems. It solves problems by separating them into easy parts and difficult Kosuke Mitarai: mitarai.kosuke.es@osaka-u.ac.jp ones, and gradually taking the effect of difficult parts into account. For weakly correlated systems, it usually gives sufficiently accurate physics of the systems. A benefit of using perturbative methods is its computational efficiency compared to the cost of exactly solving quantum systems. The computational cost to obtain an exact solution of an n-body quantum system is generally exponential to n on classical computers, while that of perturbative methods is only polynomial. Another important aspect of perturbation is its physical interpretability. It provides us insights into what effect a specific interaction of the system has on its overall physical properties.
In this work, we provide a method to implement perturbation theory on quantum computers and discuss if the benefits described above can also be obtained. Our method allows one to use strongly-interacting Hamiltonians that are only solvable with quantum computers as a starting point of the perturbation. More specifically, our algorithm first constructs a ground state of an unperturbed Hamiltonian via quantum signal processing (QSP) [1][2][3] and fixed-point amplitude amplification [4]. Then, we generate a perturbative state by applying the inverse of an unperturbed Hamiltonian with quantum signal processing (QSP), and obtain an expectation value of a perturbation operator via robust amplitude estimation (RAE) [5][6][7].
For an unperturbed Hamiltonian H and perturbation V that can be written as H = h σ and V = v σ where σ are Pauli operators, the complexity of the proposed algorithm isÕ( h 1 v 2/3 /(∆δ)) and O( h 1 v 2 2/3 /(∆ 2 δ)) respectively for the first-order and second-order perturbation, where ∆ is a spectral gap of H, δ is error in the estimated perturbation energy and a p = ( a p ) 1/p . We also perform a concrete resource analysis of the algorithm for simple chemical systems such as water clusters and polyacenes to discuss its practicality. This is, to the best of our knowledge, the first such analysis of a practical application of the QSP and the QSPbased matrix inversion technique. Despite of efficiency of QSP compared to conventional techniques, it is found that our algorithm gives impractical numbers as computational cost; for example, we estimate over 10 31 calls of block-encodings would be required to perform perturbation on a pentacene molecule. This is much larger than the cost required for naively performing the phase estimation of the total Hamiltonian, which only requires 10 10 calls of block-encodings. While we could not achieve a reduction of computational cost like the classical perturbation theory, the other benefit of perturbation, that is, the interpretability of the result, is still an important point. Conventional techniques of quantum simulations based on phase estimation can give us energy and its eigenstates, but cannot provide us insights into why the energy is the obtained value. We, therefore, believe this work is a first step toward an "explainable" quantum simulation on fault-tolerant quantum computers.

Perturbation theory
We have an n-qubit Hamiltonian H total = H + V . We consider the Hamiltonian H and the perturbation V that can be decomposed into Pauli operators σ as, Note that H and V can, without loss of generality, be defined as not having I ⊗n term. Let the ground state of H total with eigenvalue E 0 be |E 0 . Also, let an eigenstate H with an eigenvalue i be | i . We assume that the eigenvalues are ordered in ascending order 0 < 1 ≤ · · · ≤ 2 n −1 . It is well known [8] that |E 0 can be approximated as where Π = I − | 0 0 |, to the first order in V if i is not degenerate. The corresponding eigenvalue E 0 can be approximated as, where (1) and (2)

Block-encoding of a Hamiltonian
First, we introduce the block-encoding [3,9]. We say a unitary U A block-encodes a matrix A when it has the following form: More formally, we say (n + l) Babbush et al. [10] showed that we can perform phase estimation of a unitary e i arccos A/α to estimate an eigenvalue of A to precision δ by using times. The phase estimation algorithm [11] takes in a state |ψ and outputs eigenvalue e iφ of a unitary U with a probability p(φ) = | φ|ψ | 2 where |φ is the eigenstate of U corresponding to the eigenvalue e iφ . If we wish to obtain a specific eigenvalue such as the ground state energy, |ψ must therefore have non-negligible overlap with the corresponding eigenstate.
An explicit block-encoding of an n-qubit operator A which can be represented as a sum of Pauli operator as A = L =1 a σ can be constructed via PREPARE a and SELECT operations introduced in Ref. [10]. PREPARE operation acts on l = log 2 L qubits as where a 1 = |a |. SELECT acts on n + l qubits and defined as Then,

Eigenvalue transformation via quantum signal processing
Given a block-encoding U A of A, we can construct a block encoding of P (A) for certain polynomials P (x) [1][2][3]. In this work, we are only interested in real polynomials P (x). The seminal work [3,Corollary 10] shows that, for a degree-d real polynomial P (x) such that we can obtain a unitary P 0 and P 1 such that whereP (x) is a degree-d polynomial such that Re[P (x)] = P (x), * denotes complex conjugate, and |g is a "garbage" state which is orthogonal to |0 l [P (A) |ψ ]. We can construct P 0 and P 1 with d calls of U A . Moreover, a unitary defined as which uses another ancilla qubit, satisfies as shown in [3,Corollary 18]. Note that P can be constructed by only using d-calls of U A rather than 2d-calls, since P 1 is obtained by inverting phase sequences of P 0 [3, Corollary 18]. After the application of P, we can post-select on the ancilla being |0 l+1 to obtain a state proportional to P (A) |ψ . The probability of success is P (A) |ψ 2 . This procedure is called quantum signal processing (QSP) [1][2][3]. In this work, we say a polynomial is QSP-implementable if it satisfies the above two conditions (13) and (14).

Amplitude estimation
Amplitude estimation refers to various techniques to obtain values of amplitudes of a quantum state |ψ within error of δ using O(1/δ) calls of a state preparation unitary U ψ such that |ψ = U ψ |0 . The original algorithm can be found in [12]. Recent developments have made the procedure significantly more efficient. For example, a state-of-the-art method called the robust amplitude estimation (RAE) [5][6][7] can empirically estimate ψ|σ|ψ for a Pauli operator σ within a mean squared error of δ 2 by using calls of U ψ . Other recent works such as iterative quantum amplitude estimation [13,14] have shown similar performance. In this work, we employ RAE to estimate the expectation values.
3 Perturbation theory on quantum computers

High-level overview
Our approach for performing perturbation on a quantum computer is as follows: Step 1 of the algorithm can be replaced with more sophisticated techniques provided in Refs. [15,16] from the naive phase estimation. Algorithms for step 2 are also proposed in e.g. [15,17,18], but previous works discuss the cost in terms of O-notations. In the following subsections, we describe the details of steps 2, 3, and 4, without resorting to O-notations. Some assumptions we make are in order. First, we assume that one can prepare a state |ψ that has non-zero overlap p with | 0 in steps 1 and 2. We do not discuss how to choose and prepare such a state in detail because they strongly depend on the target system. For our numerical resource estimation in Sec. 4, we employ Hartree-Fock states. The second assumption is the knowledge of lower bounds to the overlap p and the spectral gap ∆ = 1 − 0 of the unperturbed Hamiltonian H. This assumption is required for constructing quantum circuits for steps 2 and 4. In the following, we formulate the cost by using p and ∆ but they can always be replaced by their lower bounds. We note that Ref. [19] has considered a similar task that involves the inversion of Hamiltonians using QSP. Their main target is, however, to calculate Green's function, which is widely used to obtain response properties with respect to external fields, and not to construct general perturbation theory, which can evaluate the energy of a many-body system, like this work. From the technical viewpoint, for the calculation of Green's function, only the operator in the form of (H − z) −1 where z is a complex number is required and one does not need to construct Π(H − 0 ) −1 Π as we do here. Ref. [20] also considers the perturbative simulation of quantum systems. Their goal is to implement real-time dynamics on a quantum computer and does not focus on obtaining the energies of Hamiltonians. It is achieved by first discretizing the time evolution to small time slices and sampling the perturbative part of each of the short-time dynamics with respect to a quasi-probability distribution. Their method hence requires exponential cost with respect to the evolution time, and would not be suitable for extracting energy eigenvalues with high precision.

Preparation of a reference state
In this work, we utilize the QSP-approximation of rectangular functions introduced by Low and Chuang [21] and reviewed in Ref. [2] to prepare the reference state | 0 from which we start the perturbation. | 0 can be generated via phase estimation. However, it is more efficient to use QSP to filter | 0 when we know the value of the corresponding eigenvalue 0 . The construction of rectangular functions closely follows that of [2,21] but we give more detailed cost, that is, the degree of polynomial needed for their approximation, than the previous works. In Appendix A, we show that there exists a QSP-implementable polynomial P filter ε,κ,x th (x) such that, where where x th = x th + κ/2 and W (x) is the Lambert W function. We plot the values of n filter (ε, κ, x th ) as Fig. 1 with various parameter settings. Using an inequality W (x) > log(x) − log(log(x)), we can roughly upperbound n filter (ε, κ, x th ) by, which has O(log(1/ε)/κ) scaling. Let us assume that we know an estimateˆ 0 of 0 such Note that, since we assume that H does not have I ⊗n term, H has L H + 1 terms. Also, note that h 1 = h 1 + |ˆ 0 |. Let l = log(L H + 1) , the number of ancilla qubits needed to block-encode H . We  Tables 2 and 3.
and the second largest energy is larger than (∆−δ 0 )/ h 1 , we can set The error parameter ε controls the fidelity of | 0 . Let us assume an initial state |ψ fed to the QSP satisfies, where | ⊥ 0 is a state orthogonal to | 0 . To obtain | 0 , we first apply P filter ε,κ,x th (H / h 1 ) via QSP and get a state in the form of, The post-selection on the ancilla being |0 results in a state which is used as the reference state for the perturbation.
We obtain the perturbation energies as expectation values of observables O = V and V Π(H − 0 )ΠV with respect to |˜ 0 . We, therefore, wish to make to obtain the perturbation energy with an accuracy of δ. In Appendix B, it is shown that, Assuming O(ε 2 ) term is negligible, we can roughly upper-bound δ prep for O = V and V Π(H − 0 )ΠV respectively as, and We can remove the need for post-selection in the preparation of |˜ 0 with the fixed-point amplitude amplification algorithm [4]. Deterministic state preparation is essential to employ advanced expectation value estimation techniques such as RAE [5][6][7]. When applied to our setting, it allows us to deterministically prepare a state where |0˜ ⊥ 0 is a state orthogonal to |0 l +1 |˜ 0 , for r > 0 with approximately . It can roughly be bounded from above by 2r V and 2r V 2 /∆. We need to take sufficiently small r to make this negligible.
Finally, let us discuss the related previous works about ground state preparation. Ge, Tura, and Cirac [17] also consider how to realize filtering operations but with a so-called linear combination of unitaries approach which originates in [22]. The approach requires log 2 (d) ancilla qubits to implement degree-d polynomial. Refs. [15,18] use QSP to implement a reflection operator I − 2 | 0 0 |, and uses the fixed-point amplitude amplification. They do not perform a detailed error analysis of the protocol as we do in this paper; we expect that it would be more involved and that such analysis would yield a comparable performance to the method presented in this section. All of the above works give the cost in O-notation and do not give detailed cost estimates as we do in this work.

First-order perturbation energy
The first-order energy correction (1) 0 can be obtained by naive measurement of the operator V , that is, we perform the VQE-like measurement where we estimate each of Pauli operator σ appearing in Eq. (2) by the RAE.
The technique demands us to deterministically prepare | 0 . To this end, we employ the fixed-point amplitude amplification algorithm [4]. Let the unitary that prepares |˜ 0 in Eq. (33) be W prep . Then, the RAE can empirically estimate ˜ 0 |σ |˜ 0 with a mean squared error δ 2 σ using calls of W prep in total [6, Eq. (8)] (converted here to be a noise-free case). Employing optimal distribution of M to minimize the overall error, we can obtain an estimate of ˜ 0 |V |˜ 0 with a mean squared error δ 2 1 using calls of W prep in total. For completeness, we derive Eq. (36) in Appendix C following Ref. [6] where they have derived the total cost for noisy quantum computers.
Here, we combine the above discussion to state a formal result in the O-notation as follows: Let | 0 be the non-degenerate ground state of H. Assume that we have an estimateˆ 0 of 0 , the ground state energy of H, such that |ˆ 0 − 0 | < δ 0 < ∆. Moreover, assume that we can preprare a state |ψ such that | 0 |ψ | 2 = p. Then, we can estimate the first-order perturbation energy 0 |V | 0 within an additive error of δ 1 by using Proof. This is obtained by multiplying Eq.
We will estimate the concrete number of calls to the block-encoding of H in Sec. 4.

Second-order perturbation energy
We approximate Π(H − 0 ) −1 Π by the method presented in [2]. Ref. [2] gives a construction of a QSPimplementable polynomial to realize matrix inversion.
However, as we show in Appendix D, it can also be used for approximating Π(H − 0 ) −1 Π. More concretely, there exists a QSP-implementable polynomial P ptb ε,w,w0 (x) that satisfies the following conditions: Such P ptb ε,w,w0 (x) can be constructed by appropriately tuning the error parameter of a matrix inversion polynomial presented in [2], as shown in Appendix D. Its degree is where n sign (ε, κ, c) We plot the values of n ptb (ε, w, w 0 ) as Fig. 2 with various parameter settings. Let P ptb ε,w,w0 be a unitary that implements P ptb ε,w,w0 (H / h 1 ), that is, for a general state |0 l +1 |ψ , it acts as, (2) 0 can be approximated by the expectation value of the operator P ptb ε,w,w0 as, 0 ≈˜  Table 3.
where the expectation is taken with respect to |0 l +1 (V | 0 ). In Appendix E, we show that, by taking it is guaranteed that Let us define If we wish to obtain the overall perturbation energy within an error of δ, we have to make Eq. (50) sufficiently smaller than δ by taking small ε and δ 0 . As for ε, we can for example take ε = ( h 1 / V 2 )(δ/10) without increasing the implementation cost so much to guarantee that δ 2 is negligible with respect to δ. Note that 1/ε only contributes logarithmically to n ptb . As for δ 0 , although its effect to n ptb is logarithmic, the phase calls of U H . We therefore wish to take δ 0 as large as possible while maintaining δ0 (∆−δ0) (2) 0 δ. Now, we discuss how to obtain an estimate of P ptb ε,w,w0 and therefore˜ (2) 0 . One strategy to measure P ptb ε,w,w0 is to generate V | 0 via the block-encoding of V . This, however, would prevent us from analyzing the contribution of each term to the perturbative energy. We take a more naive strategy to avoid this issue.
Since V can be expressed as Eq. (2), P ptb ε,w,w0 can be rewritten as, We can estimate P ptb ε,w,w0 using Hadamard test involving a controlled-(σ P ptb ε,w,w0 σ ) gate. Note that controlled-P ptb ε,w,w0 does not increase the T cost by much; it only adds a single control qubit to the controlled-U H gates already used to realize P ptb ε,w,w0 . We first deterministically generate |0 l +1 |˜ 0 using filtering described in Sec. 3.2. Then, adding an ancilla qubit, we prepare a state , which is the quantity we wish to estimate. This allows us to employ the RAE to estimate the expectation values P ptb ε,w,w0 . Denoting the unitary to prepare the state in Eq. (55) by W 2 , the optimal number of calls of W 2 to estimate P ptb ε,w,w0 with a mean squared error δ 2 P is, which follows from exactly the same discussion provided in Appendix C. Since˜ is obtained by multiplying 2/(w h 1 ) to the estimated value of P ptb ε,w,w0 , we need to take δ P = w h 1 2 δ 2 to obtain˜ (2) 0 with a mean squared error δ 2 2 . Therefore, we conclude that, calls of W 2 are needed to estimate˜ (2) 0 . Combining the above discussion, we obtain the following result: ∆ , H and p be defined as in Theorem 1. Additionally, let Π = I − | 0 0 |. Then, we can estimate the secondorder perturbation energy − 0 | V Π(H − 0 ) −1 ΠV | 0 within an additive error of δ 2 by using calls of a block-encoding of H .
Proof. This is obtained by multiplying the following equations.
We will calculate the concrete values in the next section.

Resource analysis with realistic parameters
Here, we assess the overall cost for obtaining the ground state energy E 0 within the accuracy of δ by the perturbative method presented in the previous section. We analyze a rough T-cost of our approach with an assumption that the second-order perturbation is accurate enough to produce E 0 . Let the T-cost of implementing U H be T (H). After formulating the overall cost of the proposed approach, we plug in some realistic parameters to the formula and discuss its feasibility.

Formulating the overall cost
Let us first discuss how to distribute the overall error δ to δ 0 , δ 1 , and δ 2 . Importantly, δ 2 depends on δ 0 .
The contribution to the T-cost needed for each step of the presented perturbation method can be summarized as follows:

Estimation of 0
Phase estimation of H to obtain 0 within an accuracy of δ 0 takes since we can obtain the ground state energy with probability p if we use an initial state (26). Note that if we expect

Estimation of
(2) 0 It takes M (2) (Eq. (58)) calls of W 2 which makes use of single calls of W prep and P ptb ε,w,w0 . The cost of W prep is 2 log e (2/r) √ p n filter (ε filter , κ, x th )T (H).

Application to molecular systems
We first perform resource estimation for clusters of water molecule (H 2 O) m for m = 2, 3, 4, 5, 6. Geometries of systems are taken from Ref. [23]. We use the STO-3G minimal basis set to represent the Hamiltonians, which are computed with PySCF [24,25] and Open-Fermion [26]. Hamiltonians are expressed Löwdin localized orbital [27], which allows us to separate them into intra-molecule and inter-molecule interactions. Here, we estimate the cost to implement the perturbation theory presented in the previous section treating the inter-molecule interactions as the perturbative term V . Since we expect (2) 0 ∆ in this case, we set δ 0 = δ 1 = δ 2 = δ chem /3 and δ chem = 10 −3 Hartree.
For a cluster of m molecules like the ones considered here, a special treatment for generating the unperturbed ground state | 0 can be made. The unperturbed Hamiltonian H can be written in the form of, where H i only acts on qubits corresponding to localized orbitals in the i-th molecule. Note that the above form is always possible with appropriate indexing of orbitals and Jordan-Wigner transformation [28]. Then, the ground state | 0 of H is a tensor product of ground states | 0,i of H i . Let us now assume that we can efficiently generate |ψ i such that as an initial state to create | 0,i . A naive application of the discussion in the previous section would yield an inefficient protocol that scales exponentially with respect to m since the overlap between m i=1 |ψ i and | 0 = m i=1 | 0,i decays as i p i . However, instead of directly generating m i=1 | 0,i by applying global amplitude amplification to m i=1 |ψ i , we can independently prepare | 0,i using log(2/r i )/ √ p i applications of     filtering polynomial with error parameter ε filter,i . This replaces log(2/r) √ p factor with i log(2/ri) √ pi , and does not require the exponential cost. Note that we should take r i = r/m and ε filter,i = ε filter /m to maintain sufficient accuracy for | 0 . In the following numerical results, we take this approach.
We also perform resource estimation for m-acene molecules for m = 4, 5, 6. Geometries of the molecules are taken from Ref. [29]. We use STO-3G minimal basis set, and obtain an active space consisting of πorbitals using PiOS [30]. The sizes of active spaces are 36, 44, and 52 spin-orbitals for m = 4, 5, 6, respectively. After the Jordan-Wigner transformation [28], the total Hamiltonian H total is partitioned into H and V in such a way that, if σ contained in H total has any Pauli-X or Y operators acting on inactive orbitals, σ is grouped into V , and otherwise, σ is taken into H. Note that this choice of perturbation operator V makes the firstorder energy correction 0 |V | 0 = 0. We, therefore, do not estimate the cost for first-order perturbation in this case. Since openfermion [26] is not capable of handling all spin-orbitals of relatively large molecules like polyacene, we derive the expression of molecular Hamiltonians in terms of Majorana operators, each of which corresponds to Pauli operators, in Appendix F and use it for calculations. Since we expect (2) 0 ∆ in this case, we set δ 1 = δ 2 = δ 2 = δ chem /3 and δ 0 = (∆/ MP2 )δ 2 , where MP2 is the correlation energy obtained with the second-order Møller-Plesset perturbation theory using restricted Hartree-Fock state as a reference state.
The properties of Hamiltonians are summarized in Table 1. For water clusters, p and ∆ are calculated as the overlap and energy difference between the Hartree-Fock state and the exact ground state of a single water molecule isolated from the other ones. For polyacene, ∆ is taken from Ref. [29] where the authors calculated the energy gap using density matrix renormalization group. p is taken to be 0.7 assuming that we use Hartree-Fock states as input. This value of p is determined by performing complete-active-space configuration interaction (CASCI) calculation with an active space consisting of 32 spin orbitals and 14 electrons for tetracene and pentacene. The calculation is performed with OpenMolcas version 22.02 [31,32] and we find the overlap between the CASCI and Hartree-Fock states are 0.713 and 0.722 respectively for tetracene and pentacene.
Using the numbers in Table 1, we calculate the resource according to Sec. 4.1 and show the results as Table 2 and 3. The polynomial degrees are also plotted in Figs. 1 and 2. The overall cost for the first-order perturbation ranges from the order of 10 10 to 10 18 calls of U H . That for the second-order is much higher and needs 10 18 to 10 32 calls of U H . These numbers are not practical; even for the smallest system that we considered, (H 2 O) 2 , the implementation of U H needs over 10 4 T gates at least since we find L H = 3 × 10 3 (see Sec. 2.2). This means the overall T-count of the algorithm for the first-order perturbation is over 10 14 . Even if we can implement a T gate in 1 µs, the algorithm would take 10 8 seconds. Our analysis shows the importance of concrete resource analysis of quantum algorithms beyond O notation. The algorithm, at the first sight, seems to be efficient in the sense that it hasÕ( h 1 v 2/3 /(∆δ)) andÕ( h 1 v 2 2/3 /(∆ 2 δ)) cost, but practicality cannot be ensured without such an analysis.
Finally, let us compare this cost with a naive approach that uses phase estimation of H total . We can since T-cost is roughly proportional to the number of Pauli terms [33]. To obtain the ground state energy with accuracy of δ by the phase estimation of H total , we need √ 2π( h 1 + v 1 )/(2δ) calls of U H total [10]. If a reference state used in phase estimation has overlap p with the true ground state, the correct ground state energy is obtained with probability p. The T-cost for the whole process is therefore well approximated as Taking the numbers from Table 1, naive phase estimation of H total would only need a cost equivalent to 10 10 calls of U H (calculated roughly through Eq. (64)) even for the largest molecule, hexacene, considered in this work. Unfortunately, we conclude that the perturbation theory, at least in its present form, does not reduce the computational cost of solving chemical systems.

Conclusion
In this work, we provided a quantum algorithm to obtain perturbative energies and analyzed its rough computational cost for simple molecular systems. It is, to the best of our knowledge, the first concrete resource analysis of the practical application of QSP. Several remarks are in order. First, the estimated numbers are rather pessimistic; the algorithm needs over 10 14 T gates for the simplest system considered here. However, it should be remarked that the large contribution to the overall cost comes from the expectation value estimation of the perturbation operator V . This is the same problem faced by the variational quantum eigensolvers [34] where the energy expectation values have to be determined with high accuracy. We might be able to reduce the measurement counts, M (1) and M (2) , by neglecting unimportant parts of V . For example, when using the active space approximation, we empirically know that energy corrections due to interactions between the core and virtual orbitals are small, and hence might be able to neglect them. Also, the use of other amplitude estimation techniques such as the ones presented in Ref. [14] may reduce the runtime by a constant factor.
Second, although the overall cost seems to be impractical, the polynomial degrees are on the order of only 10 8 even for the largest system we considered. Using more sophisticated techniques introduced in Ref. [35] can further reduce the requirement by a factor of 2-3. Hence, we might be able to perform the generation of the perturbed state (Eq. (4)) in a practical time scale. The proposed algorithm might therefore become useful when we are interested in measuring very few observables of the perturbed state. However, it should be noted that we need to obtain QSP phase sequences for 10 8 -degree polynomials. The state-of-the-art technique for the phase sequence finding [35] is verified up to 10 4degree polynomials but it is not clear if the method still works for such high-order polynomials.
Third, the required numbers for ε and r, which correspond to the infidelity of the prepared states, are very small and become comparable to the gate fidelities of fault-tolerant quantum computers for difficult problems considered in this work. This is also caused by large V and small δ, and the requirements are likely to be relaxed for other more simple operators. However, if we wish to explore the direction considered in this work, the values of ε and r are likely to remain at the same level. In this case, gate error rates should also be taken into account for an accurate estimation of required computational resources even in the fault-tolerant quantum computer regime.
Finally, it should be stressed again that the perturbative approach allows us to interpret the physical meaning of the results although the proposed algorithm requires more computational resources for computing total energy than the naive phase estimation approach. We believe that, while the values of energy are indeed an important quantity, the interpretability of the results is key to the practical applications of quantum simulation algorithms. This work is only a first step toward this goal, which remains to be reached in the future.
Program code to reproduce the numerical results of this paper is available at https://github.com/ kosukemtr/perturbation-resource-estimate.

A Polynomial approximation of sign functions and filtering functions
Let Low and Chuang [21] first approximates sign(x) with erf(kx) = 2 √ π kx 0 e −y 2 dy, and then approximates erf(kx) with a polynomial. We follow exactly the same strategy and give a detailed cost.
It is easy to verify its properties (20) and degree (21).

B Error analysis for the state preparation
By Eq. (20), which leads to Eq. (30).

C Optimal cost of robust amplitude estimation
We consider an optimal distribution of M to achieve an overall mean squared error δ 2 in estimating ˜ 0 |V |˜ 0 = v ˜ 0 |σ |˜ 0 . For a given distribution of M , the overall error δ 2 {M } can be written as, since estimates of ˜ 0 |σ |˜ 0 are independent from each other. We wish to minimize M = M while achieving δ 2 {M } = δ 2 . To do this, we define a Lagrangian and solve ∂L ∂M = 0. This leads to, where α = 5 √ 2 2 e 2 e−1 . Plugging this into the constraint δ 2 {M } = δ 2 , we get, which is equivalent to We, therefore, conclude the optimal distribution of M is and total calls M of W is,

D Polynomial for perturbation
First, we review the method presented in Ref.