Error-resilient Monte Carlo quantum simulation of imaginary time

Computing the ground-state properties of quantum many-body systems is a promising application of near-term quantum hardware with a potential impact in many fields. The conventional algorithm quantum phase estimation uses deep circuits and requires fault-tolerant technologies. Many quantum simulation algorithms developed recently work in an inexact and variational manner to exploit shallow circuits. In this work, we combine quantum Monte Carlo with quantum computing and propose an algorithm for simulating the imaginary-time evolution and solving the ground-state problem. By sampling the real-time evolution operator with a random evolution time according to a modified Cauchy-Lorentz distribution, we can compute the expected value of an observable in imaginary-time evolution. Our algorithm approaches the exact solution given a circuit depth increasing polylogarithmically with the desired accuracy. Compared with quantum phase estimation, the Trotter step number, i.e. the circuit depth, can be thousands of times smaller to achieve the same accuracy in the ground-state energy. We verify the resilience to Trotterisation errors caused by the finite circuit depth in the numerical simulation of various models. The results show that Monte Carlo quantum simulation is promising even without a fully fault-tolerant quantum computer.


Introduction
Solving a theoretical model is basic in physics for producing useful predictions. However, many models in quantum mechanics are computationally hard. In the 1980s, Richard Feynman conceived a way to solve this problem [1], i.e. "the possibility that there is to be an exact simulation, that the computer will do exactly the same as nature." This idea is one of the main motivations for developing quantum computing technologies. In the 1990s, Seth Lloyd proposed the Trotterisation algorithm to simulate the real-time evolution on a quantum computer [2]. In comparison with real-time evolution, the ground-state problem draws more attention as it determines the properties of matter, such as nuclei [3], molecules [4] and condensate matter systems [5,6]. We can solve the ground-state problem on a quantum computer with the quantum phase estimation (QPE) algorithm [7,8].
Trotterisation and QPE are quasi-exact, in which the algorithmic error decreases at a polynomial cost in time and qubits. However, approaching the exact solution is at a cost. It is widely believed that the implementation of these two algorithms needs a fault-tolerant quantum computer [9][10][11]. Practical fault-tolerant technologies are still up to development because of the large qubit overhead required by quantum error correction [12,13]. In this situation, quantum simulation algorithms suitable for an intermediate-scale noisy quantum (NISQ) [14] computer without error correction become important [15][16][17][18][19][20]. QPE as a ground-state solver (GSS) relies on an accurate real-time simulation (RTS), which is usually called Hamiltonian simulation. In quantum mechanics, real-time evolution is represented by unitary operators. Trotterisation is a specific realisation of the time evolution operator according to the Trotter formula. For an accurate Trotterisation, the evolution time has to be sliced into many steps, resulting in a deep quantum circuit. In this paper, we propose a quantum GSS algorithm resilient to Trotterisation errors allowing one to compute the ground state with shallow Trotterisation circuits. This situation is similar to classical computing. Some classical GSS algorithms are successful in certain models, for instance, quantum Monte Carlo (QMC) [3][4][5] and density-matrix renormalisation group [6]. However, their counterparts for generic RTS are inefficient [21,22]. In quantum computing, our GSS algorithm is more feasible than RTS with the same number of Trotter steps because GSS is still accurate when RTS has significant errors.
We propose a quantum algorithm to simulate the imaginary-time evolution [23], i.e. the imaginary-time simulation (ITS). Instead of realising the non-unitary evolution by measurement or approximating it with a unitary evolution [17-19, 24, 25], we simulate the imaginary-time evolution by randomly sampling quantum circuits. Imaginary-and realtime evolution operators (i.e. e −βH and e −iHt ) can be viewed as functions of the Hamiltonian. An integral formula connects these two functions: we can express the imaginary-time evolution operator as a weighted integral of the real-time evolution operator with different time values. The weight function is a product of Lorentz and Gaussian functions. Real-time evolution operators are randomly sampled and implemented with quantum circuits. Then, we can evaluate an observable in imaginary-time evolution as an average over real-time evolution. Although based on the Monte Carlo method, our algorithm utilising quantum circuits is free of the sign problem in conventional QMC. Note that the weight function is always positive. The sign is determined by real-time evolution. We show that the sign oscillation caused by real-time evolution can be controlled by taking the proper parameters in our algorithm. Analysing the complexity, we find that the Trotter step number (i.e. circuit depth) increases polynomially with the system size and imaginary time and polylogarithmically with the desired accuracy. This polylogarithmic scaling behaviour is obtained by cancelling Trotterisation errors according to the leading-order rotation (LOR) formula [26]. Given the imaginary-time evolution, we solve the groundstate problem following the projector QMC approach [27,28]. The circuit depth required in our GSS algorithm increases polynomially with the system size and accuracy as the same as QPE [29] in the general case, and the depth increases polylogarithmically with the accuracy when an energy gap above the ground state exits.
In addition to the complexity analysis, we demonstrate that our algorithm has strong resilience to Trotterisation errors in the numerical simulation of various quantum manybody models. In all the models, ITS and GSS in our algorithm are much more accurate than RTS, given the same number of Trotter steps. Compared with QPE, the step number can be thousands of times smaller to achieve the same accuracy. Because of the error resilience, our algorithm is a promising candidate for near-term quantum computers.

Quantum-circuit Monte Carlo for imaginary-time simulation
This section describes our quantum algorithm for simulating the imaginary-time evolution. We propose that one can realise the non-unitary imaginary-time evolution operator  We take the ten-spin one-dimensional transverse-field Ising model with the parameter λ = 1.2 as an example to generate the data; see Appendix G. The Trotter step number is N t = 20. (b) Two-time correlation along the dashed line in (a). Because of the small Trotter step number, the difference between the exact correlation and the Trotterisation result is significant. (c) Expected value of the Hamiltonian in the imaginary-time simulation. The expected value converges to the ground-state energy when β is large. In comparison with the exact value, we can find that the Trotterisation result is accurate, although the Trotter step number is small. (d) Schematic of our Monte Carlo quantum simulation algorithm. The sample generator produces random times (t, t ) according to the distribution P (t, t ) = C −2 g(t)g(t ). by randomly sampling real-time evolution operators, which are unitary and can be implemented with quantum gates. The complexity analysis shows that the required circuit depth scales with the permissible error η as O (ln η) 2 .
Without loss of generality, we assume that the HamiltonianH is a traceless operator.
In the algorithm, we take H =H − E 0 1, in which E 0 is a tunable constant. This constant does not change eigenstates and the time evolution. The ground-state energy ofH is E g . ITS is to evaluate where β is the imaginary time, |Ψ(0) is the initial state, and O is an observable. Because the time evolution operator e −βH is non-unitary, we cannot realise it directly with unitary circuits [24,25]. In our algorithm, we construct e −βH as an integral of real-time evolution operators.

Formalism
We connect real-and imaginary-time evolution operators with an integral formula where is a product of Lorentz and Gaussian functions (i.e. a product of probability density functions of Cauchy-Lorentz and Gaussian distributions, up to normalisation). If we neglect the Gaussian function (i.e. take the limit τ → ∞), the integral results in the exact imaginary-time evolution operator e −βH when H is positive semi-definite (see Appendix A); This is why we choose the Lorentz function. The Gaussian function is introduced to reduce the impact of Trotterisation errors. When the real-time evolution is realised with Trotterisation, usually the error increases with |t|. Therefore, a small g(t) at large |t| is preferred. The Gaussian function decreases exponentially with |t| and, with proper parameters, only slightly modifies the exact imaginary-time evolution operator as we show next.
With the Gaussian function, the operator realised with the integral is The error due to the Gaussian function has an upper bound . Here • 2 denotes the matrix 2-norm. We can find that the error in G(H) is small when we take proper E 0 and τ .
In our algorithm, we simulate the imaginary-time evolution by taking G(H) as an approximation to e −βH . Note that the error in the approximation can be arbitrarily small and is controlled by parameters E 0 and τ . Substituting G(H) for e −βH in Eq. (1), we obtain the approximation to O (β) in the form where and The complete error analysis is given in Sec. 2.5.

Algorithm of imaginary-time simulation
In this section, we give the details of our ITS algorithm. As depicted in Fig. 1 We generate random values of the real time according to importance sampling. The probability density of (t, t ) is taken as ) is the normalisation factor. Later, we will show that C determines the variance of the Monte Carlo computing. With the distribution P (t, t ), we can rewrite an imaginarytime correlation as the expected value of the real-time correlation, i.e.
Therefore, we can estimate O G (−iβ, iβ) by computing the average over random (t, t ).
We propose to evaluate real-time correlations O (t, t ) with the Hadamard test. A general-purpose Hadamard-test circuit requires an ancillary qubit in addition to qubits for encoding the simulated system [30]. To measure O (t, t ), we need to implement a controlled-U gate, in which the ancillary qubit is the control qubit, and U = e iHt Oe −iHt (suppose O is unitary). The controlled gate complicates the circuit, which is an undesired feature for applications on a NISQ system. For certain models, such as fermion models with particle number conservation, the ancillary qubit can be removed, and the circuit can be simplified accordingly [31,32]. In this work, we analyse our algorithm focusing on the general-purpose circuit with an ancillary qubit.
For a general observable O, we decompose it as a linear combination of Hermitian unitary operators, i.e. O = j a j O j . Here, a j are real coefficients, and O j are Hermitian unitary operators, e.g. Pauli operators. We can evaluate each O j with quantum circuits given in Appendix B. Real-time correlations have real and imaginary parts, which are measured with different circuits. For each circuit shot (implementation of the circuit and measurement for one time), the measurement outcome is a number µ R or µ I (µ R , µ I = ±1) corresponding to real and imaginary parts, respectively. The correlation is the expected value of measurement outcomes from corresponding circuits, i.e. O j (t, t ) = E[µ R ]+iE[µ I ]. We remark that the correlation depends on the parameter E 0 . The overline denotes that the correlation is computed by taking E 0 = 0, i.e. O (t, t ) = Ψ(0)|e iHt Oe −iHt |Ψ(0) . Correlations with any E 0 can be derived from E 0 = 0 according to O (t, t ) = e iE 0 (t−t ) O (t, t ). Therefore, only correlations with E 0 = 0 are evaluated on the quantum computer, and E 0 is an input parameter in the Monte Carlo estimator stage after the quantum computing, as shown in Fig. 1.
With real-time correlations, we can compute imaginary-time correlations and further compute O G (β) according to Eqs. (6) and (4), respectively. The detailed pseudocode of imaginary-time simulation is given in Algorithms 1 and 2. In the algorithms, N s denotes the number of random (t, t ), and M s denotes the number of circuit shots for each of real and imaginary parts for each (t, t ). Then, the total number of circuit shots is 2N s M s for each of O G (−iβ, iβ) and 1 G (−iβ, iβ). We remark that we can take M s = 1 [i.e. two circuit shots for each (t, t )], and in this case, the estimator of O G (−iβ, iβ) still converges to its true value in the limit of large N s . In the error analysis, we will focus on M s = 1, and the total number of circuit shots is 4N s .

Variance and sign problem
A potential issue in the Monte Carlo algorithm is the statistical error. In this section, we show that the factor C determines the variance in evaluating imaginary-time correlations. Because C ≤ 1, the variance is under control. Overall, the time cost of our algorithm is a polynomial function of the system size, evolution time and accuracy as shown in the full complexity analysis in Sec. 2.5. Despite this, we also analyse the sign problem through the average phase in this section. The sign problem is the problem of numerically evaluating Algorithm 1 Imaginary-time correlation.

4:
Generate j l with the probability a −1 Implement the circuit for M s shots to evaluate the real part of O j l (t l , t l ), and record measurement outcomes {µ R,l,k | k = 1, . . . , M s }.

6:
Implement the circuit for M s shots to evaluate the imaginary part of O j l (t l , t l ), and record measurement outcomes {µ I,l,k | k = 1, . . . , M s }.
Algorithm 2 imaginary-time simulation. the integral of a highly oscillatory function. Therefore, an average phase approaching zero indicates the sign problem. We show that the average phase is always finite under the same assumption as in QPE, i.e. the initial state has a finite overlap with the true ground state.

Variance of correlation estimators
The variance depends on the factor C. We consider the variance ofÔ in Algorithm 1, which is the estimator of O G (−iβ, iβ). Because O is Hermitian, O G (−iβ, iβ) is always real. According to Eq. (6) and Algorithm 1 ( , where µ = Re e iθ (µ R + iµ I ) , and e iθ corresponds to the phase factor sgn(a j )e iE 0 (t−t ) . Measurement outcomes take µ R , µ I = ±1, see Appendix B. Therefore, |µ| ≤ √ 2. The variance of the estimator is The parameter a O is defined in Algorithm 1, and a 1 = 1 when O = 1.
We remark that the variance is amplified when using the zeroth-order leading-order rotation formula [26] to implement the real-time evolution operator (see Sec. 2.4.2), which will be analysed in Sec. 2.5.

Sign problem
The sign problem refers to the problem that the Monte Carlo summation is taken over oscillatory values. If the values are complex, this problem is also called the phase problem. When the sign problem occurs, it usually becomes difficult to achieve the desired accuracy due to statistical errors. We can directly analyse the statistical error in O (β), which will be given in Sec. 2.5. Here, we discuss the sign problem through the average phase of values in the summation. We will show that the average phase is always finite, i.e. the oscillation is not severe.
We focus on the denominator Now, we consider the case that the imaginary-time evolution operator is accurate: We take a proper E 0 and large τ such that G(H) − e −βH 2 is small. In this case, of the ground state in the initial state. Then, the average phase becomes According to this result, we prefer to take E 0 close to E g . Because G(H) − e −βH 2 decreases exponentially with (E g − E 0 ) 2 τ 2 , it is allowed to take an E 0 close to E g when τ is large. We also prefer a small C, which coincides with the analysis of the variance. With proper parameters, we can have E[e iϕ ] pg √ 2 . Therefore, as long as p g is finite, i.e. there is a finite overlap between the initial state and ground state, the average phase is finite. This assumption of finite overlap is generally required in projector QMC algorithms [27,28], in which the imaginary-time evolution operator projects the initial state onto the ground state; and it is the same in quantum algorithms that find the ground state by a projection, e.g. the QPE algorithm [7,8,33].
To understand how the sign oscillation is controlled in our algorithm, we express the initial state in the form |Ψ(0) = √ p g |ψ g + 1 − p g |ψ e . Here, |ψ g is the ground-state, and |ψ e denotes the excited-state component in the initial state. In our algorithm, we compute the imaginary-time evolution as an integral of the real-time evolution. In realtime evolution, the state is e −iHt |Ψ(0) = √ p g e −i(Eg−E 0 )t |ψ g + 1 − p g e −iHt |ψ e . When E 0 is close to E g , the phase of the ground-state term varies slowly with time. This term results in the finite phase average.
In the LOR formula, the real-time evolution operator is exact at the cost of a variance increasing with the simulated real time. In this case, we apply a truncation on the real time to control the overall variance. We focus on the LOR formula in the complexity analysis, which results in a circuit depth increasing polylogarithmically with the desired accuracy. In the Trotter formula, the real-time evolution operator is inexact, and the error increases polynomially with the real time. In this case, the truncation is unnecessary. In Sec. 4, we will show that our algorithm is resilient to the Trotterisation error in real-time evolution operator.

First-order Trotter formula and ordering operation
We consider a Hamiltonian in the formH = M j=1 H j , where H j are Hermitian operators. The first-order Trotter formula reads By cutting the evolution time into N t slices (i.e. Trotter steps), we useŨ (t) ≡ [S 1 (∆t)] Nt to approximate e −iHt , where ∆t = t/N t . The error in this approximation scales as O(t 2 /N t ).
To analyse the impact of the Trotterisation error in our ITS algorithm, we introduce the ordering operation to express the first-order Trotter formula. We define operators where N t is the number of Trotter steps: H kM +j is the H j operator in the (k + 1)th Trotter step. We use P to denote the ordering operation according to the label of operators H j (j = 1, 2, . . . , N t M ): For any product of H j operators, the ordering operation is defined as where K is the total number of operators in the product, and k j = K i=1 δ j,j i is the number of the H j operator in the product.
To define P for a general function of H j operators, we consider the spectral decomposition of each H j , i.e. H j = n j ω j,n j Π j,n j , where Π j,n j is the orthogonal projection onto eigenvector of H j with the eigenvalue ω j,n j . We can find that Therefore, for a general function, we define the ordering operation as Note that for an analytic f , we can also read P {f (H 1 , H 2 , . . . , H NtM )} as applying the ordering operation on each term in the Taylor expansion of f (H 1 , H 2 , . . . , H NtM ).
Using the ordering operation, we can reexpress the first-order Trotter formula as where

Zeroth-order leading-order-rotation formula
Instead of approximating the real-time evolution operator with a product of unitary operators, we can also approximate it with a summation of unitary operators [36,37]. We can even reproduce the exact real-time evolution operator with a summation formula including infinite terms, which can be implemented with the Monte Carlo method. Here, we take the zeroth-order LOR formula [26] as an example.
We suppose that each term in the Hamiltonian is a Pauli operator, i.e.H = M j=1 h j σ j . Here, h j are real parameters, and σ j are Pauli operators. The zeroth-order LOR formula reads where We use s to denote terms in the summation formula. Each term is a rotation operator or Pauli operator. We define unitary operators We define weights and phases as where weights w s (t) are positive, and phases θ s (t) are real. Then the formula can be rewritten as The zeroth-order LOR formula is exact and free of the Trotterisation error. However, when we realise such a formula using the Monte Carlo method, the variance is amplified by a factor of C A (t) 2 , and C A (t) = s w s (t) = 1 + C 2 To reduce the factor, we divide the time t into N t Trotter steps, and the formula becomes , i.e. we can reduce the variance by taking a large N t .
We can find the connection between the zeroth-order LOR formula and Trotterisation by neglecting Pauli-operator terms in the formula. We have With only rotation terms, the formula is a summation of stochastic Trotterisation products. The role of Pauli operators in the LOR formula is correcting errors in the summation of stochastic Trotterisation products, which is at the cost of increased variance.
Algorithm 3 Imaginary-time correlation with the zeroth-order leading-order-rotation formula.
Implement the circuit for M s shots to evaluate the real part of

6:
Implement the circuit for M s shots to evaluate the imaginary part of

Complexity analysis of imaginary-time simulation
In this section, we present a complete complexity analysis of our algorithm, which is possible because the imaginary-time evolution operator is constructed according to an explicit integral formula. In comparison, there are algorithms that approximate the imaginarytime evolution operator with a unitary operator utilising variational circuits [17][18][19]. With the explicit formula, we can avoid an algorithmic accuracy depending on the variational ansatz. Compared with the quantum-classical Monte Carlo algorithm, in which one simulates the imaginary-time evolution with the classical auxiliary-field Monte Carlo [28] incorporating a quantum trial state [20], our algorithm is free of the sign problem even without introducing the phaseless approximation and can approach the unbiased solution in a systematic way. As summarised in Theorem 1 (also see Theorem 3), the conclusion is that our algorithm requires resources scaling polynomially or polylogarithmically with the system size (assuming local interactions), imaginary time and desired accuracy.
In Algorithm 1, we have assumed that one can implement the exact real-time evolution operator directly by unitary gates, which corresponds to Trotterisation in the N t → ∞ limit. To optimise the performance of our algorithm in the complexity analysis, in this section we focus on the zeroth-order LOR formula. Accordingly, the algorithm has to be adapted, see Algorithm 3. We will explicitly give C T in the algorithm after introducing the truncation on the real time. The circuit for measuring

Truncation error
If we take the LOR formula to realise the real-time evolution operator, the variance is amplified by a factor of C A (∆t) 2Nt , which increases with |t|. In order to bound the variance, we apply a truncation at t = ±T in the integral. Given N t , the integral formula with truncation becomes Replacing where and P S, The truncation error decreases exponentially with T because of the Gaussian function in g(t). Using properties of the Gaussian function and Lorentz function, we have and is the upper bound of the truncation error. The truncation is a universal way to bound the impact of errors in the real-time evolution. Suppose a method, e.g. the first-order Trotter formula, realises the unitary operator U (t) as an approximation to e −iHt , and (t) is the upper bound of is not larger than one. Therefore, as long as the circuit complexity for controlling the error (T ) in the real-time evolution is a polynomial function of T , the circuit complexity for imaginary-time evolution is also polynomial. Specifically for the LOR formula, the realtime evolution operator is exact, therefore the contribution of the (T ) term is zero; The cost is an increased variance, which we will discuss next.

Statistical error
According to Algorithm 3 (M s = 1), the variance ofÔ has the upper bound The derivation of the upper bound is similar to evaluating Eq. (6) according to Algorithm 2, see Sec. 2.3.1. Now we work out an upper bound of C T . Because C A (t) increases monotonically with |t|, we have Here, we have used that the integral of g(t) is not larger than one. The upper bound of C T approaches one when N t is large, therefore, we can control the variance by taking a large N t . The statistical error in the correlation estimator is . According to Chebyshev's inequality, the probability that the error e O is not smaller than a O δ has an upper bound In the rigorous complexity analysis, we will use Chebyshev's inequality. However, it is worth noting that if we approximate the distribution ofÔ with the normal distribution according to the central limit theorem, the error e O is not smaller than a O δ with the probability The computation fails when the statistical error in the denominator of O (β) is not smaller than δ or the statistical error in the numerator of O (β) is not smaller than a O δ. The failure probability has the upper bound This probability determines the sampling cost of our ITS algorithm.

Circuit depth and sampling cost
There are three error sources in our ITS algorithm: the error due to the Gaussian function γ G , the truncation error γ T and the statistical error determined by C T . The bound γ G holds under the condition E g − E 0 ≥ β τ 2 , therefore, we need to take an E 0 smaller than the actual ground-state energy. Because the energy is evaluated according to Eq. (1), all errors are amplified when the denominator Ψ(0)|e −2βH |Ψ(0) is small. The denominator has the lower bound Ψ(0)|e −2βH |Ψ(0) ≤ p g e −2(Eg−E 0 )β , therefore, we prefer an E 0 close to E g . In summary, we first need an estimate of the ground-state energy and then take E 0 accordingly. There are two cases. First, we have a sufficiently accurate estimate of E g obtained from other algorithms, such as the Hartree-Fock method or QMC on a classical computer. In this case, we can directly use it in our ITS algorithm. Second, if we do not have a sufficiently accurate estimate of E g , we can start with any preliminary estimate and iteratively improve the accuracy of the ground-state energy with our ITS algorithm, which will be given in Sec. 3.2. Here, we simply assume that there is an estimate of the ground-state energyÊ g with the uncertainty δE, i.e. E g ∈ [Ê g + δE,Ê g − δE]. We also assume that we have a lower bound of the ground-state probability p b , i.e. p g ≥ p b .
Let η be the total error. In the algorithm, we take τ ∼ ln 1 η because γ G decreases exponentially with τ 2 . Similarly, we take T ∼ ln 1 η because γ T decreases exponentially with T 2 /τ 2 . To control the variance, we need to take N t ∼ T 2 ∼ ln 1 η 2 . The circuit depth is proportional to the number of Trotter steps N t , therefore, the circuit depth increases polylogarithmically with the desired accuracy. The costs of our ITS algorithm are summarised in Theorem 1. The proof is in Appendix C. See Theorem 3 for the case without a prior estimate of the ground-state energy. Theorem 1. Let η and κ be any positive numbers. Suppose thatÊ g , δE and p b satisfy conditions E g ∈ [Ê g + δE,Ê g − δE] and p g ≥ p b . The result of O G (β) from the ITS algorithm with the zeroth-order LOR formula isÔ 1 . If we take proper parameters in the algorithm, the inequality holds with a probability higher than 1 − κ, with the Trotter step number and sample size The circuit depth and dependence on the system size can be derived from the theorem. To evaluate the zero-order LOR formula, the gate of each Trotter step is a controlled U s (∆t) gate, where U s (∆t) is a rotation gate or a Pauli gate as defined in Eq. (17). Such a controlled gate can be decomposed into O(n) controlled-NOT and single-qubit gates [26], where n is the number of qubits for encoding the simulated system (i.e. the system size). Then the total number of gates for evaluating the zero-order LOR formula (the two controlled U S (t) gates, see Appendix B) is O(nN t ). Some additional gates (e.g. single-qubit gates on the ancillary qubit) are used in the circuit, which however does not change the polynomial scaling with N t and n. Note that gates for preparing the initial state are not taken into account, and we have to assume that the number of these gates scales polynomially with n. In addition to this direct dependence on n, h tot also depends on the system size. Usually, h tot scales with n polynomially in Hamiltonians that only involve local interactions (see Appendix G for example Hamiltonians). Therefore, the overall dependence on the system size is polynomial. Besides the circuit depth and sample size, the qubit cost in our algorithm is up to one ancillary qubit in addition to the n qubits representing the system.
Later, we will consider the first-order Trotter formula in addition to the zero-order LOR formula. For the Trotter formula, the circuit depth is also proportional to N t . When we use the first-order Trotter formula to approximate the real-time evolution operator, each Trotter step corresponds to a controlled S 1 (∆t) gate, see Eq. (10). If each term H j is a Pauli operator, S 1 (∆t) is a product of M (the number of terms in the Hamiltonian) rotation gates. Therefore, we can realise the controlled S 1 (∆t) gate with O(M n) controlled-NOT and single-qubit gates, and the total gate number is O(M nN t ). Similar to the parameter h tot in the LOR formula, the term number M usually depends on the system size and scales with n polynomially in Hamiltonians that only involve local interactions.

Monte Carlo quantum ground-state solver
We can compute the ground-state energy by simulating the imaginary-time evolution. Specifically, we take O =H in Eq. (1) and evaluate the energy in the imaginary-time evolution. When the time β increases, the energy approaches the ground-state energy, i.e. lim β→∞ H (β) = E g , under the assumption that the probability of the ground state in the initial state |Ψ(0) (i.e. p g ) is finite.
In this section, we first analyse the projection error, i.e. the error due to a finite β. Then, we present an iterative algorithm for computing the ground-state energy and analyse its complexity. Our ITS algorithm requires a preliminary estimate of the ground-state energy. In the iterative algorithm, we start with an estimate with a large uncertainty, then we let the uncertainty decrease exponentially with the number of iterations. In addition to the iterative method, we will introduce two other approaches in Sec. 5: the optimisation of E 0 and quantum subspace diagonalisation.

Projection error
In the projector method, a projection onto the ground state is applied on the initial state to compute the ground-state energy. We approximate the projection with e −βH . The final state reads e −βH |Ψ(0) = D n=1 √ p n e −β(Eg+En−E 0 ) |φ n , where D is the dimension of the Hilbert space, p n is the initial probability in the eigenstate |φ n with the eigenenergy E g + E n , i.e. E n ≥ 0. We suppose that n = 1 corresponds to the ground state, i.e. p 1 = p g , |φ 1 = |Ψ g and E 1 = 0. Then, the expected value of energy reads .
where α = pg 1−pg , and we have used that D n=2 p n = 1 − p g . The second term in Eq. (35) is the error due to the imperfect projection. Note that H (β) ≥ E g .
From Eq. (35), we can find that the expected value converges to the true ground-state energy in the limit β → ∞, as long as p g is finite. The speed of convergence depends on the energy gap above the ground state, ∆ = min n=2,...,D E n . Without the gap, the projection error is inversely proportional to the evolution time and has the upper bound With a finite energy gap, the projection error decreases exponentially with time. When 2β∆ ≥ 1 + ln(1 + e −1 α −1 ), the upper bound becomes See Appendix D for proofs of the upper bounds.

Iterative ground-state solver and its complexity
In the iterative ground-state solver, we need an initial estimate of the ground-state energy. There is a simple and universal initial estimate that always works. Because H 2 ≤ h tot , the ground-state energy is in the interval [−h tot , h tot ]. Therefore,Ê g = 0 and δE = h tot is an estimate of the ground-state energy satisfying E g ∈ [Ê g − δE,Ê g + δE] (as required in Theorem 1); We take this as the initial estimate.
With the initial estimate, in each round of iteration, we choose parameters in ITS such that the uncertainty δE is reduced by a fixed factor. We take the factor of 3 4 as an example, see Algorithm 4. Details of choosing parameters in the algorithm are given in Appendix E.
In the general case (without assuming a finite energy gap), the upper bound of the projection error decreases linearly with 1 β . To reduce the error in the ground-state energy to ξ, we need to take β ∼ 1 ξ . In this case, the circuit depth N t ∼ β 2 ∼ 1 ξ 2 scales polynomially with the desired accuracy. If there is a finite energy gap above the ground state, the upper bound of the projection error decreases exponentially with β. Then, β ∼ ln 1 ξ , and the circuit depth N t ∼ ln 1 ξ 2 scales polylogarithmically with the desired accuracy. The costs of our iterative GSS algorithm are summarised in Theorem 2. The proof is in Appendix E. Choose β such that the projection error [Eq. (36) in the general case or Eq. (37) with a finite energy gap] is smaller than δE 2 .

5:
Take parameters in ITS such that the ITS error h tot η is smaller than δE 4 . See Appendix E.

6:
Implement ITS according to Algorithm 2. Theorem 2. Let ξ and κ be any positive numbers. Suppose that p b satisfies p g ≥ p b . The result of ground-state energy from the iterative GSS algorithm isÊ g . If we take proper parameters in the algorithm, the inequality |Ê g − E g | < h tot ξ holds with a probability higher than 1 − κ, with the largest Trotter step number N t,max and total sample size The largest Trotter step number depends on the energy gap: (i) In the general case, (ii) If there is a finite energy gap ∆ above the ground state, where ∆ b is an input parameter satisfying ∆ ≥ ∆ b . Now, we reconsider our ITS algorithm for computing any observable O. If we do not assume prior knowledge about the ground-state energy, we can work out it with our iterative GSS. Because our GSS algorithm includes ITS as a subroutine, let's call ITS for computing O the task ITS for clarity. Given the evolution time β of the task ITS, the accuracy of the ground-state energy at the level of δE ∼ 1 β is required. Before implementing the task ITS, we can run the iterative GSS algorithm taking ξ ∼ 1 htotβ . We can find that the cost in iterative GSS is independent of the desired accuracy η in the task ITS. Therefore, incorporating iterative GSS does not change the dependence on η. The proof of the following theorem is given in Appendix F.

Theorem 3. ITS incorporating iterative GSS.
Let η and κ be any positive numbers. Suppose that p b satisfies p g ≥ p b . If we take proper parameters in the algorithm, the inequality holds with a probability higher than 1 − κ, with the largest Trotter step number and total sample size where ζ = min η, 1 htotβ .

Error resilience -Suppressed error distribution in the frequency space
Our first result about the resilience to errors is the complexity analysis, which shows that the circuit depth (which is proportional to the Trotter step number N t ) scales polylogarithmically with permissible errors η and ξ (depending on the gap in GSS), in contrast to the polynomial scaling [18,29]. Instead of the zeroth-order LOR formula considered in the complexity analysis, in the following, we show that our algorithms are also resilient to Trotterisation errors in the first-order Trotter formula, and the error resilience reduces the circuit depth for implementing Trotterisation. In this section, we give an explanation of this resilience, and in the next section, we demonstrate it numerically with various models.

Imaginary-time evolution operators with Trotterisation errors
Before explaining the error resilience, we first work out the imaginary-time evolution operator with Trotterisation errors according to the integral formula. The real-time evolution operator realised with the Trotter formula in the spectral decomposition is j=1 ω j,n j ∆t Π NtM,n N t M · · · Π 2,n 2 Π 1,n 1 . According to the integral formula in Eq. (2) (without truncation), the imaginary-time evolution operator with Trotterisation errors reads In addition to the integral formula, we can also realise the imaginary-time evolution operator with a product of non-unitary operators according to the Trotter formula [18,19].
To realise e −βH , the non-unitary product is e E 0 β S 1 −i β Nt Nt , where N t is the number of Trotter steps, and S 1 (t) is defined in Eq. (10). When t is imaginary, S 1 (t) is non-unitary.

Three operators for comparison
We show the error resilience in the comparison between three operators: the real-time evolution operator V R = e −iHt realised with the Trotter formula, the imaginary-time evolution operator V I = e −βH realised with the non-unitary product and the (approximate) imaginary-time evolution operator V G = G(H) realised with the integral. To simplify expressions, we define three functions f R (ω) = e −iωt , f I (ω) = e −βω and f G (ω) = G(ω). Then we can express error-free operators as V α = f α (H), where α = R, I, G.

The three operators with Trotterisation errors areṼ
Using the spectral decomposition, we can also express these operators in a unified form,Ṽ α = P {f α (H sum − E 0 1)}.
The Trotterisation error in each operator isṼ α − V α .
In Fig. 2, we plot magnitudes of errors in three operators taking a one-qubit Hamiltonian as an example. The figure shows that the error in V R decreases rapidly with N t especially when N t is small, which opens a large gap from V R and V I , i.e. the integral formula of imaginary-time evolution operator is much more robust to Trotterisation errors.

Error distribution in the frequency space
To explain the error resilience, we consider the Fourier decomposition of the Trotterisation error. For the real-time evolution, Here, l is the label of the frequency ω l , and v l are operator-valued coefficients. The Fourier decomposition of V R is V R = n e −i(Eg+En−E 0 )t |φ n φ n |, where E g + E n are eigenvalues ofH. The Fourier decomposition ofṼ R is given by Eq. (44), notice thatṼ R = e iE 0 tŨ (t). Therefore, Nt NtM j=1 ω j,n j }. To obtain a similar decomposition of the Trotterisation error in V I , we can simply replace t with −iβ. We havẽ Therefore, for all three operators, there is a unified expression of Trotterisation errors Note that the operator-valued coefficients v l are the same for three operators, and the Trotterisation error is determined by the function f α . The Fourier decomposition of Trotterisation errors explains the error resilience. We plot three functions f α in Fig. 3. The absolute value of f R (ω) is always one. The function f I (ω) is smaller than one when ω > 0 and larger than one when ω < 0. In comparison, f G (ω) is always smaller than one: The function takes its maximum value at ω = 0 and Algorithm 5 Ground-state solver by optimising E 0 .

3:
Take E 0 = 0, computeĤ according to Algorithm 2, and save all the data.
It is similar for Algorithm 3.
decreases exponentially to zero as |ω| increases. The function f G (ω) suppresses the impact of v l . As a result, the error in V G is much smaller than in V R , but the error in V I can be even larger than in V R , as shown in Fig. 2. This observation suggests that our algorithm based on the integral formula is more robust to Trotterisation errors than ITS algorithms based on the non-unitary Trotter formula [18,19].
To verify this explanation, we plot magnitudes of errors in three operators as functions of E 0 , see Fig. 4. Because f G (ω) is smaller than one for all ω, the impact of all v l is suppressed regardless of E 0 . Therefore, we expect to observe the error suppression (compared with V R ) for all E 0 , which coincides with the numerical result in Fig. 4.

Numerical demonstration of the error resilience
In this section, we present numerical results about error resilience. First, we compare three computation tasks, RTS, ITS and GSS. We show that the impact of Trotterisation errors is smaller in ITS and GSS compared with RTS. Then, we study the accuracy of GSS when the system size increases. We find that a Trotter step number scales linearly with the system size is sufficient for computing the ground-state energy of the one-dimensional transverse-field Ising model. Finally, we compare our GSS algorithm to QPE. For the water molecule, our algorithm reduces the Trotter step number from 6 × 10 5 in QPE to only four. Similar results are also obtained with the one-dimensional transverse-field Ising model.
Instead of the iterative method introduced for rigorous complexity analysis, we propose an alternative method to determine E 0 via the variational principle and optimisation. This method may be more relevant to practical implementation than the iterative method. Our algorithm effectively prepares the state G(H)|Ψ(0) , and its energy E(β, τ, E 0 ) = H G (β) is a function of parameters β, τ and E 0 . According to the variational principle E(β, τ, E 0 ) ≥ E g , we can solve the ground state by minimising the energy, i.e. we take E min = min E 0 E(β, τ, E 0 ). In this paper, we focus on varying E 0 , and in general, we can also vary β and τ to minimise the energy.
We have the following remarks on the optimisation method. First, G(H)|Ψ(0) as a trial wavefunction is capable of expressing the ground state: When the Trotter step number N t increases, E min always approaches E g . The reason is that the operator G(H) is a projection onto the ground state. As shown in Sec. 3 Figure 5: Errors R , I and G in quantum simulation. Results of the ten-spin one-dimensional transverse-field Ising (1D TFI), 3 × 3 two-dimensional transverse-field Ising (2D TFI), ten-spin antiferromagnetic Heisenberg (AFH) and six-site (12-qubit) Fermi-Hubbard (FH) models are obtained in numerical simulation. The first-order Trotter formula is implemented for all the models, and the second-order Trotter formula is implemented for AFH and FH models. In numerical simulation, we take N t = 20, T = 4 for the FH model and T = 3 for other models, τ = 2T , E 0 = E g in real-and imaginary-time simulations and β = T in ground-state solver. In quantum subspace diagonalisation (QSD), we take d = 8 and β a = aT /d. approach the ground state when N t scales polynomially (or polylogarithmically if there is a gap). Second, the one-parameter optimisation is technically trivial by a grid search. Third, we can implement the quantum computing only for E 0 = 0 and then generate results for all E 0 , see Algorithm 5. Therefore, the optimisation is entirely on the classical computer without iteratively querying the quantum computer, and we can take a high grid resolution without increasing the quantum-computing cost. Note that variational quantum algorithms [15][16][17][18][19] usually involve a feedback loop between quantum and classical computers to update parameterised quantum circuits. Our optimisation method only uses one-way communications from the sample generator to the quantum computer to the Monte Carlo estimator, see Fig. 1. This one-way feature removes the feedback loop and potential latency.
In the numerical study, we focus on the optimisation method and Trotter formula instead of the iterative method and LOR formula. Quantum subspace diagonalisation (QSD) or quantum Lanczos [18,[43][44][45][46] is a way to reduce the error in GSS, and we also demonstrate this method in the numerical study. Following Ref. [18], we choose a set of imaginary times {β a | a = 1, 2, . . . , d} to generate a subspace. We evaluate matrices A a,b = 1 (−iβ b , iβ a )

Comparison of three tasks
To demonstrate the error resilience in our ITS and GSS algorithms, we compare the impacts of Trotterisation errors in three tasks: RTS, ITS and GSS. Specifically, we compute the two-time correlation H (t, t ) in RTS, the expected value of energy H (β) in ITS and the ground-state energy E g in GSS. Using the first-order Trotter formula, results with Trotterisation errors are denoted by H (t, t ), H (β) and E g , respectively: and E g depends on the method. Two methods are considered. The first method is the op- The second method is QSD incorporating the optimisation of E 0 , in which imaginarytime correlations A a,b are B a,b are replaced by their Trotterisation approximations. In numerical calculations, we neglect implementation errors in quantum computing, i.e. errors due to imperfect quantum gates and statistical errors. We have analysed statistical errors in Secs. 2.3.1. We demonstrate the error resilience with various quantum many-body models, including the one-and two-dimensional transverse-field Ising models, anti-ferromagnetic Heisenberg model, Fermi-Hubbard model and randomly generated Hamiltonians. Several different parameter values (denoted by λ) are taken in each model. See Appendix G for details of these models and numerical calculations.
First, we illustrate the comparison by taking the ten-spin one-dimensional transversefield Ising model as an example. We take the model parameter λ = 1.2 (coupling strengths are on the order of one), the maximum evolution time T = 3, τ = 2T and the Trotter step number N t = 20. The two-time correlation evaluated using the first-order Trotter formula is shown in Fig. 1(a). Data on the diagonal line, i.e. H (t, −t), are redrawn in Fig. 1(b) for a comparison between correlations with and without Trotterisation errors. We can find that the error in H (t, −t) is already comparable to the variation of the function itself, i.e. the Trotter step number is inadequate for even approximate RTS. Using these inaccurate RTS data, we can compute the expected value of energy in the imaginary-time evolution, and the result (assuming E 0 = E g ) is shown in Fig. 1(c). We can find that the ITS is still relatively accurate. As shown in Fig. 1(e), the minimum expected value of energy (for β = T ) is close to the ground-state energy.
For a quantitative comparison, we consider three quantities: the average error in RTS, Numerical results show that with the same Trotter step number, the error in GSS after the optimisation of E 0 is smaller than the average error in ITS, and the average error in ITS is smaller than the average error in RTS. For the ten-spin one-dimensional transversefield Ising model with λ = 1.2, we have R = 1.4, I = 0.062 and G = 0.0021. We can find that I is much smaller than R , and G is much smaller than R . This observation holds in various models as shown in Fig. 5. We can find that I is smaller than R by a factor of 10 to 200 in almost all cases, except the Fermi-Hubbard model with λ = 0.2, and G is smaller than R by a factor of 30 to 10 5 . We also show that QSD incorporating the optimisation of E 0 can significantly reduce the error. Similar results are obtained with randomly generated Hamiltonians as shown in Fig. 6. Intuitively, the error in GSS should be comparable to the error in ITS because the ground-state energy is computed with ITS. However, the error in GSS is actually smaller. There are two reasons for this result. First, the error in ITS is large when β is small, which causes a relatively large average error. See the inset of Fig. 1(b). There are two sources of the error: the integral formula [i.e. the error G(H) − e −βH ] and Trotterisation. Usually, the Trotterisation error increases with time. On the contrary, the formula error decreases with β. The operator G(H) is a projection onto the ground state when β is large (assume E 0 ≤ E g ). Therefore, when β is large, the difference between G(H) and e −βH is small because they are both projections. The error at a small β is mainly due to the integral formula. In a finite interval, this error decreases with β; If β is too large, the error increases with β because of Trotterisation. Second, taking the optimal E 0 in GSS reduces the error.

Scaling with the system size
We use the one-dimensional transverse-field Ising model to test the scaling behaviour of our GSS algorithm and infer the impact of Trotterisation errors in hundred-qubit quantum computing. We increase the number of spins n spin and take the number of Trotter steps N t = rn spin . The numerical result shows a trend that the error in the ground-state energy decreases with n spin when r is a fixed constant, see Fig. 7. Accordingly, considering N t = 50 for n spin = 100, we can infer that the error is smaller than 0.0025, which is much smaller than the energy gap ∼ 0.8 between the ground and first-excited states (given the model parameter λ = 1.2). Therefore, we can solve a hundred-qubit ground-state problem and achieves a relatively small error with tens of Trotter steps. Here, E 0 plays the role of ω: Thinking of the limiting case τ → ∞, G(H) becomes e −β|H−E 0 1| , therefore, G(H)|Ψ(0) 2 reaches a maximum when E 0 takes an eigenvalue. In the limit of a large Trotter step number N t , all these algorithms can produce an accurate result at a polynomial cost. Next, we show numerical evidence that our algorithm is much more robust than QPE when N t is small. To demonstrate robustness, we consider the water molecule encoded into fourteen qubits. We plot the energy H (β) and its normalisation factor G(H)|Ψ(0) 2 in Fig. 8. We can find that G(H)|Ψ(0) 2 has a peak around the exact E g . If we take E 0 at the maximum of the peak as the ground-state energy, the error is about 0.09 hartree. The energy H (β) has a minimum. If we take the minimum value as the ground-state energy (i.e. the E 0 optimisation method), the error is about 0.00026 hartree. The error in the minimum-energy approach is below the chemical accuracy (about 1 millihartree) and smaller than the error in the maximum-peak approach for two orders of magnitude. In the calculation, we use the second-order Trotter formula and take N t = 4 Trotter steps. See Appendix G for results of the first-order Trotter formula and other step numbers. For a direct comparison, QPE achieves the chemical accuracy with the step number N t ≈ 6×10 5 for the same molecule [9].

Comparison to quantum phase estimation
In addition to the water molecule, we also observe the robustness to Trotterisation errors in computing the ground state of the transverse-field Ising model. In Table 1 Table 1: Error G and corresponding Trotter step numbers in our quantum-circuit Monte Carlo (QCMC) algorithm [without and with quantum subspace diagonalisation (QSD)] and the quantum phase estimation (QPE) algorithm. Here, we list the result of the one-dimensional transverse-field Ising model with λ = 1.2 and n spin = 20. We use the first-order Trotter formula for both algorithms.
list the error in the ground-state energy and corresponding Trotter step numbers for the model with n spin = 20 spins. With N t = 40 Trotter steps in our algorithm, errors in energy are 9.1 × 10 −4 and 2.8 × 10 −5 depending on whether QSD is used. To achieve the same energy resolutions using QPE, Trotter step numbers 2.5 × 10 5 and 4.6 × 10 7 are required, respectively. Therefore, the circuit depth of our algorithm is thousands to a million times shallower.

Conclusions
This paper proposes a quantum algorithm for simulating the imaginary-time evolution by sampling random quantum circuits. By analysing how the finite circuit depth impacts the accuracy, we find that our algorithm is resilient to Trotterisation errors caused by the finite circuit depth. The error resilience is demonstrated in two ways, complexity analysis and numerical simulation. In the complexity analysis, we find that the circuit depth scales polylogarithmically with the desired accuracy. We have this superior scaling behaviour owing to the LOR formula and Gaussian function in the integral. In the LOR formula, Trotterisation errors are corrected by random Pauli operators leading to an exact real-time evolution operator. With the Gaussian function, the truncation time scales polylogarithmically with the desired accuracy. These two factors together result in the polylogarithmically-scaling circuit depth. Based on the simulation of imaginary time, our algorithm for solving the ground-state problem inherits the resilience to Trotterisation errors. If there is a finite energy gap above the ground state, the circuit depth for solving the ground state also scales polylogarithmically with the desired accuracy. This energy gap is a finite energy difference between the ground state and first-excited state, which exists in many finite-size Hamiltonians, rather than a finite gap in the limit of large system size (a stronger condition). In the numerical simulation, we directly compare our algorithm to the QPE algorithm and find that the circuit depth is thousands of times smaller in our algorithm. This reduction in the circuit depth can be explained by analysing the Trotterisation error in the frequency space. Optimising the algorithm parameter E 0 and utilising QSD can further reduce the impact of Trotterisation errors.
In this paper, we focus on applying our algorithm for computing the ground state. We can also use imaginary-time evolution to study finite-temperature properties [52,53]. Constructing the imaginary-time evolution operator according to the Monte Carlo method, our algorithm can be used as a subroutine and combined with conventional projector QMC algorithms, such as Green's function Monte Carlo and ancillary-field Monte Carlo [27,28]. In this way, we may further reduce the circuit depth by using more classical computing techniques and resources. Our algorithm can be completely explicit as in the iterative approach or includes a minimised variational computing, i.e. the optimisation of E 0 . Compared with quantum variational algorithms [15][16][17][18][19], our optimisation is in a one-parameter space and implemented entirely on the classical computer without involving quantum computing. It is worth noting that variational principles are efficient tools for maximising the power of shallow circuits. We can think of a variational Monte Carlo algorithm in which we optimise the distribution of circuits g(t) rather than circuit parameters. The distribution function in the Monte Carlo quantum simulation provides a new dimension to explore to develop efficient quantum algorithms in the NISQ era. As we show in this paper, the problem caused by shallow circuits can be solved to a large extent by using Monte Carlo methods in quantum computing.

A Integral formula
The integral formula reads First, we apply the spectral decomposition to the Hamiltonian, and we get where {ω} are eigenvalues of the Hamiltonian, which are real, and {|ω } are orthonormal vectors. Then, the real-time evolution operator reads The integral formula becomes and According to the matrix 2-norm, the error in the integral formula is Here, we have used that For the Lorentz-Gaussian function, we have where Next, we use the residue theorem to evaluate this integral. We consider the contour in the Here, we have used properties of the Faddeeva function. When β + ηωτ 2 = 0, we have When β + ηωτ 2 < 0, we have Here, we have used that erfc(x) + erfc(−x) = 2. Therefore, for all cases, we have Because g(t) ≥ 0, the normalisation factor is C = G(0) = erfc( β √ 2τ ). To derive the error in the integral formula, we consider β − ωτ 2 < 0. Using erfc(x) ≤ e −x 2 when x ≥ 0, we have When ∆E ≥ β τ 2 , we have β − ωτ 2 < 0 for all ω. Then, where γ G ≡ e − ∆E 2 τ 2 2 is the upper bound of the error due to the finite τ .

B Circuit
We can evaluate O (t, t ) with or without an ancillary qubit. The protocol with an ancillary qubit works for the general case, and the protocol without the ancillary qubit only works under certain conditions. We present both protocols in this section, however, we focus on the general-case protocol in the complexity analysis.
The circuit for the general-case protocol is shown in Fig. 9, which is adapted from Refs. [30]. We assume that O is a unitary operator, e.g. a Pauli operator. For a general operator, we can express it as a linear combination of unitary operators and measure each term. In the circuit, the gate U i prepares the initial state, i.e. |Ψ(0) = U i |0 ⊗n . The top qubit is the ancillary qubit. The gate B is for adjusting the measurement basis: B † ZB = X and B † ZB = Y to measure X and Y Pauli operators of the ancillary qubit, respectively. Let X and Y be expected values of ancillary-qubit Pauli operators evaluated using the circuit, then O (t, t ) = X + i Y . Therefore, the measurement outcome of X is µ R , and the measurement outcome of Y is µ I .

C Proof of Theorem 1
First, we analyse the total errors in estimators of O (−iβ, iβ) and O (β), respectively. Then, we give the proof of Theorem 1.
hold with the probability 1 − P , where P ≤ Proof. G is an integral over unitary operators, therefore, G 2 ≤ C ≤ 1. In the proof, we will also use O 2 ≤ a O . Note that |Ψ(0) is a normalised state. There are three error sources. First, we use G(H) to approximate e −βH . The corresponding error is where we have used the condition E g − E 0 ≥ β τ 2 . Second, we use G T (H) to approximate G(H). The corresponding error is Third, the statistical error in the estimateÔ is and the inequality holds with the probability 1 − P (e O ≥ a O δ). The total error inÔ is Therefore, the inequality in the lemma holds with a probability larger than 1 − P (e O ≥ a O δ). The bound of the probability is given by Eq. (29). The lemma has been proved.
holds with the probability 1 − P , where P ≤ Proof. The error inÔ 1 is Here, we have used that If e 1 < δ and e O < a O δ, we have d 1 < and d O < a O . Therefore, the inequality in the theorem holds with a probability larger than 1 − P δ . The bound of the probability is given by Eq. (31). The lemma has been proved.
The following is the proof of Theorem 1, which contains a protocol for choosing parameters in our ITS algorithm. The protocol is up to optimisation but sufficient for working out the scaling behaviour of our algorithm.

Proof.
Step-1 -We take Step-2 -We solve the equation With the solution, we have Step-3 -In the error budget, three error sources contribute equally. We take δ = 1 3 and solve the equation to work out Later, we will choose parameters such that γ G , γ T ≤ x.
Step-4 -We take Under this condition, the upper bound γ G holds. We have Step-5 -We take then Step-6 -We choose a value of C max , which is larger and close to 1, e.g. C max = 1.1, and solve the equation Step-7 -We take

D Bounds of the projection error
To work out an upper bound of the projection error, we consider the functional of the weight function w(x), where f (x) ≡ e −x x α+e −x . Taking w(x) = D n=2 p n (α + e −x ) δ(x − 2βE n ), we can express the error as H (β) − E g = (2β) −1 y. Because y ≤ max x f (x), we have The derivative of the function is The maximum value of the function is at x = x 0 , which is the solution of the equation α(1−x)+e −x = 0. We can calculate the solution via the product logarithm. Given x 0 , the maximum value is max Instead of using the exact solution, we consider x 1 = 1 + ln(1 + e −1 α −1 ), and we have α Replacing max x f (x) with ln(1 + e −1 α −1 ) in Eq. (90), we can obtain the upper bound in Eq. (36) With the gap, the upper bound is given by max x∈[2β∆,∞) f (x). We assume that 2β∆ ≥ The upper bound with a finite gap is which is the same bound as in Eq. (37) according to the definition of function f .

E Proof of Theorem 2 E.1 The general case
Each iteration. First, we give details of how to choose parameters in ITS in each round of iteration. LetÊ g and δE be outputs of the previous round. We take parameters as follows.
Step-2 -We take η = δE 4htot and choose parameters in ITS according to the protocol in the proof of Theorem 1. Note that a O = h tot when O =H. In each round of iteration, we set the failure probability upper bound as κ N i instead of κ, where N i is the number of iterations.
According to Theorem 1, the error in ITS is smaller than δE 4 , i.e.
with a probability higher than 1 − κ N i . Then the total error is smaller than 3δE 4 with the same probability lower bound 1 − κ N i . Taking β according to δE, we have Substitute Eq. (94) into Eqs. (33) and (34), we obtain and Total cost. The initial estimate isÊ g = 0 and δE = h tot . To reduce the uncertainty δE to the desired accuracy h tot ξ, we take the number of iterations N i = ln ξ ln 3

4
. Because for each iteration the failure probability has the upper bound κ N i , then the total failure probability is lower than κ.
The cost of each iteration increases with δE −1 , therefore, the last step has the largest cost. Substituting δE = O(ξh tot ) into Eq. (95), we obtain N t of the last step, which is N t,max . The sample size of the last step is Note that the factor ln 1 ξ is due to N i . The total sample size N s,tot is smaller than N i N s,max .

E.2 The case with a finite gap
Given a lower bound ∆ b of the energy gap, the imaginary-time evolution with is sufficient to reduce the projection error to the desired level. However, before implementing ITS with the imaginary time β ∆ , we have to work out a preliminary estimate of the ground-state energy with sufficient accuracy. This can be achieved following the general-case approach.
With the finite energy gap, the algorithm has two stages. In the first stage, we follow the approach in Appendix E.1 to reduce the uncertainty δE to instead of the ultimate desired accuracy h tot ξ. We use to denote this intermediate desired accuracy. Then, the cost in the first stage is (102) In the following, we assume that ξ > ξ to work out the scaling with ξ.
In the second stage, we take β = β ∆ and η = ξ 2 . We choose parameters in ITS according to the protocol in the proof of Theorem 1.
Because p b is the lower bound of p g , we can always take p b = 3 3+e if the input lower bound to the algorithm is higher than 3 3+e (for the simplicity of the proof). Then, p b ≤ 3 3+e always holds. Under this condition, ln 4 α b ≥ x 1,b , where x 1,b = 1 + ln(1 + e −1 α −1 b ). Because α b ≤ α, we have x 1,b ≥ x 1 and ln 4 α b ≥ x 1 . When ξ is a small number, ln 4 α b ξ ≥ ln 4 α b and 2β∆ ≥ x 1 . Then, we can apply Eq. (37), and where we have used that ∆ ≤ 2h tot . Note that H (β) − E g is always positive. With η = ξ 2 , the total error is smaller than h tot ξ.

F Proof of Theorem 3
First, we take ξ = 1 βhtot in iterative GSS (without the energy gap assumption). Then the final uncertainty of the ground-state energy is δE = 1 β . We choose the sample size in iterative GSS such that it fails with a probability lower than κ/2. Substituting ξ into Eqs. (39) and (38), we can work out the cost of iterative GSS.
Second, with the ground-state energy and the uncertainty δE = β −1 , we implement the task ITS. The cost of task ITS is given by Eqs. (33) and (34). Note that = O e −4 η . We also choose the sample size in task ITS such that it fails with a probability lower than κ/2. Then, the overall failure probability is lower than κ.
Consider both iterative GSS and task ITS, the largest Trotter step number is the Note that ξ = 1 βhtot .

G Details of the numerical simulation
Three models are considered in the numerical study. The Hamiltonian of the transversefield Ising model isH The Hamiltonian of the anti-ferromagnetic Heisenberg model is h is plotted, and the nuclear repulsion energy is not taken into account. The first-order formula with 16 Trotter steps and the second-order formula with 4 Trotter steps are sufficient for achieving the chemical accuracy. Note that we find these adequate step numbers by doubling the step number each time, i.e. they are not necessarily the minimum step numbers for the chemical accuracy. Quantum subspace diagonalisation is not used in the calculation.
For the Fermi-Hubbard model, we take H 1 = H X , H 2 = H Y and H 3 = H Z . For randomly generated models, each H j is a Pauli-operator term in the Hamiltonian.
We take the initial state as follows. For the transverse-field Ising model, |Ψ(0) = where |V ac = |0 ⊗n spin denotes the vacuum state. For randomly generated models, |Ψ(0) = |0 ⊗n spin . For the water molecule, we compute the ground-state energy in a minimal STO-3G basis of 10 electrons in 14 spin orbitals as the same as in Ref. [9]. The Hamiltonian of electrons in the water molecule at bound length 0.9584Å and bound angle 104.45 • is generated and encoded into qubits using qiskit nature [56]. The unit of energy is hartree (E h ). In our algorithm, we take β = 3 E −1 h and τ = 2β. For the initial state, we remove two-particle terms from the original Hamiltonian and calculate the ground state of the Hamiltonian with only one-particle terms, and we take the ground state of the one-particle Hamiltonian as the initial state. The first-and second-order Trotter formulas are used in our simulation, and results are shown in Fig. 10.
In the numerical simulation, we neglect quantum-machine errors (e.g. decoherence) and statistical errors, and we aim at an 'exact' computation of the integral over time t. In the numerical integration, we use the simplest midpoint rule with the step size δt = T /20. The numerical integration is truncated at t = ±10T .
To obtain a stable inverse in the numerical calculation, we apply a truncation on eigenvalues of A [46]. We suppose eigenvalues are λ 1 , λ 2 , . . . , λ d in descending order, and the number of eigenvalues greater than t = 10 −10 is d t . Then, the inverse matrix √ Λ −1 is replaced by the d t × d matrixΛ t with elementsΛ t;i,j = δ i,j λ −1 i . Accordingly, H ef f is a d t × d t matrix.
In QPE, eigenvalues of a Hamiltonian are estimated by measuring the phase e −iEnt due  Figure 11: (a) Error in the ground-state energy of the one-dimensional transverse-field Ising model with n spin spins. The ground-state energy is measured using quantum phase estimation, and the time evolution is implemented using Trotterisation with the step size δt. The black curves are obtained by fitting the data using = aδt 2 . (b) The factor a as a function of n spin . The black curve is obtained by fitting the data using a = un spin + v. Note that the data point of n spin = 2 is not used in fitting.
to real-time evolution. Here, E n is an eigenvalue, and t is the evolution time. To achieve the energy resolution , the required evolution time is t ∼ π −1 [9], and it is similar for time series analysis [48]. We follow the approach in Ref. [9] to analyse the impact of Trotterisation errors in QPE. In Trotterisation, the operator S 1 (δt) is implemented to approximate the exact time evolution operator e −iHδt . Therefore, the spectrum of the effective HamiltonianH = i δt ln S 1 (δt) is measured in phase estimation. Then, the error is the difference between ground-state energies of H andH. For the one-dimensional transverse-field Ising model, we obtain the error for varies n spin and δt by numerically simulating Trotterisation, and results are plotted in Fig. 11. By fitting the data, we find that the error scales with n spin and δt in the form = aδt 2 (118) and a = 0.228605n spin + 0.132962.
Therefore, given , we take δt = /a. The Trotter step number is N t = t/δt = π √ a √ . Now, we can estimate Trotter step numbers in QPE. The results are shown in Table 1. Note that we take = G to compare Trotter step numbers in QPE and our algorithm.