qSWIFT: High-order randomized compiler for

Hamiltonian simulation is known to be one of the fundamental building blocks of a variety of quantum algorithms such as its most immediate application, that of simulating many-body systems to extract their physical properties. In this work, we present qSWIFT, a high-order randomized al-gorithm for Hamiltonian simulation. In qSWIFT, the required number of gates for a given precision is independent of the number of terms in Hamiltonian, while the systematic error is exponentially reduced with regards to the order parameter. In this respect, our qSWIFT is a higher-order counter-part of the previously proposed quantum stochastic drift protocol (qDRIFT), whose number of gates scales linearly with the inverse of the precision required. We construct the qSWIFT channel and establish a rigorous bound for the systematic error by using the diamond norm. qSWIFT provides an algorithm to estimate given physical quantities by using a system with one ancilla qubit, which is as simple as other product-formula-based approaches such as regular Trotter-Suzuki decompositions and qDRIFT. Our numerical experiment reveals that the required number of gates in qSWIFT is signiﬁcantly reduced compared to qDRIFT. Particularly, the advantage is signiﬁcant for problems where high precision is required; for example, to achieve a systematic relative propagation error of 10 − 6 , the required number of gates in third-order qSWIFT is 1000 times smaller than that of qDRIFT.


I. INTRODUCTION
Hamiltonian simulation is a key subroutine of quantum algorithms for simulating quantum systems.Given a Hamiltonian H = L =1 h H , where h ≥ 0 and L is the number of terms, the task in Hamiltonian simulation is to construct a quantum circuit that approximately emulate time evolution U (t) := exp(−iHt) of the system for time t.Several approaches have been established for this task.The conventional approach uses the Trotter-Suzuki decompositions that provide a deterministic way for Hamiltonian simulation [1][2][3].The gate count of this approach scales at least linearly with the number of terms L in H [3]; we note that the gate count in [2] scales at least quadratically with L. Although this scaling is formally efficient but is impractical for many applications of interest, particularly for the electronic-structure problem in quantum chemistry, where the number of terms in a Hamiltonian is prohibitively large.An alternative approach is randomly permuting the order of terms in the Trotter-Suzuki decompositions [4].This randomized compilation provides a slightly better scaling for gate count over Ref. [2], but the gate count still depends on the number of Hamiltonian terms quadratically.
The quantum stochastic drift protocol (qDRIFT) [5] is another randomized Hamiltonian simulation approach but is independent of the number of terms.In qDRIFT, gates of the form exp(−iH τ ) with small interval τ are applied randomly with a probability proportional to the strength h of the corresponding term in the Hamiltonian.qDRIFT improves upon the Trotter-Suzuki approach in that its gate count is independent of L and Λ := max h (magnitude of the strongest term in the Hamiltonian), and instead depends on λ := L =1 h .However, qDRIFT has poor scaling with respect to the precision ε in contrast to that in the Trotter-Suzuki approach.We note that there are other approaches to Hamiltonian simulation with asymptotically better performance as a function of various parameters [6][7][8][9][10][11].Still, the approaches based on product formulae, e.g., Trotter-Suzuki decompositions and qDRIFT, are preferred for their superior performance in practice [12] and predominant usage in experimental implementations [13][14][15] due to their simplicity and the fact that they do not require any ancilla qubits.From this perspective, we focus on an approach based on product formulae while having better gate scaling than the previous methods.
In this paper, we propose the quantum swift protocol (qSWIFT), a high-order randomized algorithm having (i) better scaling with respect to the precision ε and (ii) the same scaling with respect to λ compared to qDRIFT.Specifically, the gate count of qSWIFT scales as O((λt) 2 /ε 1 K ) with K as the order-parameter while that of qDRIFT scales as O((λt) 2 /ε).For example, with respect to the precision ε, the gate count of qSWIFT scales as O(1/ √ ε) for the second-order and as O(1/ε 1 3 ) for the third-order.Our qSWIFT algorithm shares its simplicity with the other approaches based on product formulae.It works in the system with one ancilla qubit (we refer to the qubits other than the ancilla qubit simply as the system qubits).We can construct all gate operations with exp(iH τ ) and the swift operators Sb := S (b) ρS (b) † where S (b) is a unitary transformation; the swift operators can be constructed if we can efficiently implement controlled-H gates, as shown in Fig. 1.In the case of qDRIFT, the entire time evolution is divided into segments, and a sampled time evolution exp(iH τ ) is performed in each segment (see Fig. 2a).In qSWIFT, we utilize the swift circuit in addition to the circuit for qDRIFT.The swift circuit also has segments; in most segments, a sampled time evolution exp(iH τ ) is performed to the system qubits, but in the other segments, a sequence of the swift operators is performed (see Fig. 2b).The number of swift operators is upper bounded by about twice the order parameter.Therefore, the qSWIFT algorithm can be performed with almost no additional resources compared to the qDRIFT.
FIG. 1: Quantum circuits for implementing the swift operators S(b ).The top line corresponds to an ancilla qubit, and the bottom line corresponds to the qubits in the system.FIG.2: Examples of quantum circuits used for qDRIFT and qSWIFT.
We will now describe in more detail how qSWIFT is carried out.First, we build the qSWIFT channel that simulates the ideal time evolution.We then establish a bound for the distance between the qSWIFT channel and the ideal channel, quantified with the diamond norm, which exponentially decreases by increasing the order parameter.The established bound yields the desired scaling for the gate count of qSWIFT.It should be noted that the qSWIFT channel itself is not physical in the sense that it is not a completely-positive and tracepreserving (CPTP) map.Nevertheless, we can employ the qSWIFT channel to develop a procedure for measuring a physical quantity of interest, i.e., computing the expectation value of some given observable that is exponentially more precise than the original qDRIFT with respect to the order parameter.
Our numerical analysis also reveals the advantage of our qSWIFT algorithm.We show the asymptotic behavior of qSWIFT by using electronic molecular Hamiltonians with the number of qubits ∼ 50 and compare the performance with the other approaches based on the product formulae.Specifically, we compute the required number of gates to approximate the time evolution with the molecule Hamiltonians with a given systematic error ε.We show that the number of gates in the third-order version of qSWIFT is 10 times smaller than that of qDRIFT when ε = 0.001 for every time region.A significant reduction of the number of gates is observed when ε = 10 −6 ; the required number of gates in third-order (sixth-order) qSWIFT is 1, 000 (10, 000) times smaller than that of qDRIFT.We also simulate our qSWIFT algorithm and the other product formulae-based algorithms by using a quantum circuit simulator with the small-size (eight-qubits) molecular Hamiltonian.Its result is consistent with the result of the asymptotic behavior analysis.
The rest of the paper is organized as follows.In Section II, we briefly review approaches for the Hamiltonian simulation based on the product formula.Section III and Section IV are dedicated to proposing and analyzing our qSWIFT algorithm.In Section III, we introduce the way of constructing the second-order qSWIFT algorithm.
Then we generalize the algorithm to the higher-order in Section IV.In Section V, we validate our algorithm by numerical experiments.Finally, in Section VI, we conclude with some discussions.

II. BACKGROUND
This section covers the key background pertinent to the following sections.We begin with a brief description of Hamiltonian simulation and Trotter-Suzuki formulae in Section II A. Then we review the qDRIFT algorithm for Hamiltonian simulation in Section II B.

A. Hamiltonian simulation by Trotter-Suzuki formulae
We begin with a brief description of the Hamiltonian simulation.For a given time-independent Hamiltonian of the form H = L =1 h H , where h > 0 and H are Hermitian operators with H = 1, the task in Hamiltonian simulation is to find a good approximation of the transformation with U (t) := e iHt and t as a real parameter.We assume we can efficiently implement each e iH t by quantum gates with t as a real number.For example, if we decompose H to the sum of the tensor products of the Pauli operators, we can efficiently implement each e iH t .The conventional approach for the Hamiltonian simulation is the Trotter-Suzuki decomposition [1,2].In the first-order Trotter-Suzuki decomposition for Hamiltonian simulation, the entire simulation for time t is divided into r segments of simulations for time t/r as U (t) = (U (t/r)) r and U (t/r) is approximated as U (t/r) ≈ U which yields the approximation U (t) ≈ (U ts (t/r)) r for the entire simulation.The second-order Trotter-Suzuki decomposition is given by U (t) ≈ (U which serves as the base case for the recursive formula for the 2k th -order decomposition, where p k := 1/(4 − 4 1/(2k−1) ).
Let us discuss the Trotter-Suzuki decomposition in the channel representation.The 2k th -order Trotter-Suzuki channel U (2k) ts (t) is used to approximate the channel U(ρ) as U(ρ) ≈ (U (2k) ts (t/r)) r .For a given channel C, we denote by C r the r -repetition of C. Previous analytic work [2,5] shows that where || • || is the diamond norm.We note that α here is defined so that rαL is the number of gates used in the 2k th -order decomposition.Hence the gate count for the 2k th -order Trotter-Suzuki decomposition, denoted by G ts , is Notice that the gate count approaches to O(L 2 Λt) by increasing the order parameter 2k, but the prefactor scales exponentially with 2k.Because of this rapidly growing prefactor, Trotter-Suzuki decompositions of finite orders, typically second (k = 1) or fourth order (k = 2), are used in practice [5].We note that [3] improves the gate-count scaling to linear with respect to L though it depends at least linearly on λ = L =1 h .However, we use the gate-count scaling in Eq. ( 6) to compare qSWIFT against qDRIFT [5], particularly for comparing the numerical experiments.

B. Hamiltonian simulation by qDRIFT
Developed by Campbell [5], qDRIFT is an algorithm for Hamiltonians simulation using a randomized procedure.While the procedure is randomized, with many repetitions the evolution stochastically drifts towards the target unitary.Specifically, the exact time evolution is approximated by N repetitions of the qDRIFT channel E N as U ≈ E N N , with the qDRIFT channel defined as where is the unitary channel that we call the time operator, and are three variables used through the paper.To realize the qDRIFT channel E N , the index is sampled according to the probability p and the the quantum state ρ is evolved through the channel associated with the operator e iH τ .For evaluating the systematic error of the approximation, they define the exact short time evolution U N as Then, they show where the diamond distance is defined as They utilize the diamond distance d U, E N N as the measure of the systematic error.By using the subadditive feature of the diamond distance, they obtain the bound for the diamond distance as: In other words, to reduce the systematic error within ε, we need to set N ∈ O((λt) 2 /ε).
In most of the applications of the Hamiltonian simulation, what we have interests is computing the expectation value of an observable after applying the time evolution operator U. Let us write the expectation value as where Q is an observable and ρ init is an input quantum state.By using the qDRIFT algorithm, we can approximately compute the value of Q as where the systematic error is bounded as

III. SECOND ORDER QSWIFT
In this section, we describe our second-order qSWIFT as a preparation for introducing the general high-order qSWIFT in Section IV.To elucidate our algorithm, we use a "mixture function" in our second and higher-order qSWIFT.We begin by describing this function in Section III A. Next, we construct the second-order qSWIFT channel and discuss its error bound in Section III B. Finally, in Section III C, we explain how to apply the constructed qSWIFT channel for computing physical quantities.

A. Mixture function
In constructing our qSWIFT channels, we make use of a mixture function.As a preparation, let us first define the following sorting function.
Definition III.1 (Sorting function).Let S N be the permutation group.For positive integers k, N with k < N , let A := (A 1 , . . ., A k ) and B := (B 1 , . . ., B N −k ).We define the sorting function as where σ ∈ S N and The mixture function is defined by using the sorting function as follows.
Definition III.2 (Mixture function).For positive integers k, N with k < N , let A := (A 1 , . . ., A k ) and B := (B 1 , . . ., B N −k ).Then we define the mixture function as where the set S sub N,k is the subgroup of the permutation group S N comprised of all elements σ ∈ S N that satisfies the following condition: if both X i , X j ∈ A or ∈ B then σ(i) < σ(j) for any i < j.We remark that the number of elements in S sub N,k is N k .
For simplicity, if elements of B are identical, we denote the sorting and mixture functions as and also elements of B, are identical.In this case, X j = A for j ≤ k and X j = B for j ≥ k + 1.
Note that the sorting function in Eq. ( 17) and mixture function in Eq. ( 19) are bilinear functions.For example, if the th element of A is a linear combination of elements of another vector F, i.e., if A = n c n F n for c n ∈ C, then we have

B. Second-order qSWIFT channel
To construct the qSWIFT channel, let us define where the variables λ, p and τ are defined in Eq. ( 9).We then have for the ideal time-evolution channel in Eq. ( 10) and the qDRIFT channel in Eq. ( 7), where Using the definition of ∆ k and Eqs. ( 24) and (25), we have where we used ∆ 2 = (τ 2 /2)L (2) + ∆ 3 and linearity of M 1,N −1 to obtain the last equality.Let us denote the first two terms as We refer to E (2) as the second-order qSWIFT channel.In the following lemma, we provide a bound for the error in approximating the ideal channel U in Eq. ( 1) by the second-order qSWIFT channel E (2) , where the error is quantified as the diamond norm of their difference.
Lemma III.3.Let U be the ideal channel in Eq. (1) and let E (2) be the second-order qSWIFT channel in Eq. (32).
Then, in the region λt ≥ 1, We provide the proof in Appendix A 1. By this lemma, as far as we consider the reasonable parameter region 2) ≤ ε.This provides a quadratic improvement over the original qDRIFT with respect to ε.

C. Implementation of the second-order qSWIFT channel
The second-order qSWIFT channel we constructed is not a physical channel as it is not a CPTP map.Therefore, we cannot directly implement the qSWIFT channel itself.However, in most applications of the Hamiltonian simulation, our interest is in computing physical quantities, i.e., the expectation value of an observable after applying the time evolution operator as described in Section II B. Thus, in the following, we focus on how to compute q (2) := Tr(QE (2) (ρ init )), for a given observable Q and the input ρ init .By using q (2) , the systematic error is reduced to where we use (33), which is the quadratic improvement in terms of (λt) 2 /N from the one in qDRIFT (16).
Here we provide a way of computing q (2) by quantum circuits.We can expand q (2) as where For computing q (1) , we just need to apply the original qDRIFT channel and compute the expectation value.Thus, we focus on how to compute δq in the following.More specifically, we will show that δq is computable by using the swift circuits, composed of the time evolution exp (iH τ ) and the swift operators shown in Fig. 1.By using which can be derived from the definition (19), we obtain where . Now let us move on to the evaluation of δq r .We transform δq r as where we use (28) in the second equality.Let two probability distributions be 1 ( ) = where the input is a vector with two elements ( 1 , 2 ).Also, let Then it holds where we define a channel K(r, ) as To evaluate each term of (43), we utilize a system with one ancilla qubit.Let us write the density matrix for a given system with one ancilla qubit as the matrix form: We define the operation of a quantum channel K(r, ) that transforms the initial state ρinit = |+ +| ⊗ ρ init into a final state as where the dot [.] in the diagonal element denotes a matrix we do not have interest in now.Then it holds and consequently, where with X as the Pauli-X observable for the ancilla qubit.Therefore, we can evaluate δq r if we implement the channel K(r, ) by using quantum circuits.
We can specify the channel K(r, ) as follows: where we define ẼN := 1 ⊗ E N and = ( 1 , 2 ), and where the operation of the channel L is specified for the input having the identical non-diagonal elements as follows: The channel L can be written as the sum of two swift operators S(0) and S(1) introduced in Fig. 1 as which is easily checked by using It can be pedagogically shown that the channel (50) reproduces the transformation in (46); with ρ init = ρ init /2, it holds By substituting (52) to (50), we obtain Finally, from (38), (48), and (58), where in the second equality, we define We can evaluate (61) using Monte Carlo sampling.To this end, let be the probability distribution for product of multiple qDRIFT probability distributions, where n ∈ Z + specifies the number of qDRIFT distributions and k j , for each j, goes from 1 to L. Then where Tn ( k) is the unitary channel defined as the sequence of the time operators: We obtain an unbiased estimator of δq(s, b 1 , b 2 ) as follows: 1. Sample uniformly from {0, 1, 2. With probability P ) to (61), we obtain the estimate of δq as δ q.It should be noted that the swift circuit for evaluating each term of (65) has the structure that two swift operators are tucked in between N − 1 time operators as in Fig. 2b.Therefore, the number of gates in the swift circuit is almost the same as the original qDRIFT that requires N time operators.

IV. HIGHER-ORDER QSWIFT
In this section, we generalize the second-order qSWIFT channel introduced in Section III B to an arbitrary highorder channel.First we construct the higher-order qSWIFT channel and discuss the error bound in Section IV A. Then, in Section IV B, we construct an algorithm to apply the qSWIFT channel for computing physical quantities.

A. Higher-order qSWIFT channel
To construct a high-order qSWIFT channel, we retain higher orders of τ in the right-hand-side of where M k,N −k is the mixture function defined in Eq (19).To this end, we note that ∆ 2 is a linear combination of L (n) as per Eq. ( 28).Thus by Eq. ( 21) we obtain where δ[i, j] here is the Kronecker δ.To show the second equality, we use the identity which holds for a fixed set of integers {n j } k j=1 with n j ≥ 2. In the last equality, we use the fact that the Kronecker δ is zero if any of n 1 , n 2 • • • n k is larger than ξ.Truncating the upper limit of ξ yields a high-order qSWIFT channel.Specifically, we define the high-order qDRIFT as follows.
Definition IV.1 (Higher-order qSWIFT).We define K-th order qSWIFT channel (K ≤ N ) as Here we note that the upper limit N in the second summation can be replaced with K for K ≤ N because the terms with k > K become zero by the Kronecker δ.We remark that setting K = 2 yields the second-order qSWIFT channel in Eq. (32).Also, by setting K = 1, we reproduce the qDRIFT channel.
We now provide a bound for the error in approximating the ideal channel U in Eq. ( 1) by the high-order qSWIFT channel E (K) in the following lemma.
Lemma IV.2.Let U be the ideal channel in and let E (K) be the K-th order qSWIFT channel.Then, in the region λt ≥ 1, as far as The proof is given in Appendix A 2. As in the case of the second order, as far as we consider the reasonable

B. Implementation of the higher-order qSWIFT channel
As we discuss in Section III C, the qSWIFT channel is not a physical channel, but we can apply it to computing the physical quantities.Specifically, we discuss how to compute q (K) := Tr QE (K) . (73) Then, the systematic error is bounded as which can be exponentially small with the order parameter.
Since q (1) is again evaluable by using the original qDRIFT, we focus on the evaluation of δq (k) ( n) in the following.Similar to the second-order case, we will transform δq (k) ( n) as a sum of the terms evaluable with quantum circuits.The mixture function M k,N −k (•) included in (77) is written as the summation of N k terms as where with f σ,k,N −k as the sorting function defined in (17).Now we move on to the calculation of δq To this end, let We can write L (n) as a linear combination of D s : as per (28).Let two probability distributions: where n specifies the number of vectors in and we write the ath element of as a .We see that for n = 2, (84) is consistent with (41).Then it holds with which is consistent with (42) for n = 2.By using the bilinearity of the sorting function, we obtain where we use (83) in the first equality and (85) in the second equality.Substituting above into (81), we obtain where As in Section III C, we compute the right hand side of (91) by using the system with one ancilla qubit.Let By repeatedly operating (51) with j = 1 , • • • , n , we obtain the operation of Ln ( ) as for a given density operator ρ.Recall the definition of the sorting function (17): with Then in the system with one ancilla qubit, it holds with Since it holds that the operation of Xσ(1) Therefore, the estimation value of the observable Q = X⊗Q with the final state in (99) gives δq Next, we discuss the way of evaluating the right hand side of (100).By substituting (52) to (100), we obtain with b ∈ {0, 1} ⊗n , where Then with b j ∈ {0, where we use (101) in the first equality and we use and in the second equality.By substituting (105) to (100), we obtain Finally, combining (80), (81), and ( 90) with (108), where and we define the coefficient as We can get an unbiased estimator of (111) using Monte Carlo sampling.More specifically, we repeat the following procedure and compute the average of the output; the p-th operation works as follows: 1. Sample σ uniformly from all elements of S sub N,k .

Sample
with N shot measurements.We set the result to δ q(k) As in the second-order case, we repeat the above process N sample times and compute the average of δ q(k) we obtain the estimate of δq (k) ( n) as We write the above algorithm to estimate δq (k) ( n) as Evalcorrection and summarize it in Algorithm 1.By using the Evalcorrection, we can construct the algorithm to compute q (K) according to (75).We write the algorithm as qSWIFT and summarize it in Algorithm 2. Since we can tune N sample and N shot depending on δ q(k) ( n) we parameterize them as N sample ( n) and N shot ( n).We also write the number of circuits sampled and the number of running each circuit for calculating q (1) as N 0 sample and N 0 shot .

6:
With N shot measurements, estimate the following by the quantum circuit and set it to δ q(k) end for end for 10: end for ) by sampling N 0 sample circuits and measuring each circuit N 0 shot times; set the result to q(1) .2: Set delta = 0. for n ∈ G2(k, ξ) do

6:
Add the result of Evalcorrection (k, n, N, N sample ( n), N shot ( n)) to delta.end for 9: end for 10: Set q(K) = q(1) + delta.Output: q(K) It should be noted that there is a statistical error in q(K) due to the use of Monte Carlo sampling and the limited number of measurements.To reduce the statistical error, we need to increase N sample ( n), N shot ( n), N 0 sample , and N 0 shot .We show how the total number of quantum circuits runs scale with the order in Appendix B even though our focus in this paper is discussing the reduction of systematic error and not the statistical error.In the discussion, we show that the dominant source of the quantum circuit runs comes from the estimation of ∆q (k) ( n) with limited n and the cost for ∆q (k) ( n) with other n is negligible when N is large; consequently, the total number of quantum circuit runs scale only less than quadratically with K.

V. NUMERICAL EXPERIMENTS
In this section, we present two numerical simulations of qSWIFT.In Section V A, we show the asymptotic behavior of qSWIFT and compare it with Trotter-Suzuki decomposition and qDRIFT.In Section V B, we perform a numerical experiment with the hydrogen molecule Hamiltonian and compare its performance with those of the previous algorithms.

A. Asymptotic behaviour
We compute the required number of gates N to achieve a given systematic error for each evolution time and each algorithm: qSWIFT, Trotter-Suzuki decomposition, and qDRIFT.To clarify the difference from the original qDRIFT, we use the three molecules used in the numerical experiment of the original paper [5]: propane, ethane, and carbon dioxide.
For the qSWIFT algorithm, we utilize the third-order (K = 3) qSWIFT.We also use the sixth-order qSWIFT (K = 6) as a reference.For the Trotter-Suzuki decomposition, we use both the deterministic and the randomized methods of the first, second, and fourth order.For calculating the systematic error of qSWIFT, we use the bound (A22) derived in Appendix A 2. For the error of qDRIFT and the Trotter-Suzuki decompositions, we utilize the same upper-bound formulae as in the literature [5] (See Appendix B and Appendix C of [5]).FIG.3: The asymptotic behaviors of N for each time for each molecule to achieve the systematic error ε = 0.001.Three subfigures correspond to each molecule: propane with STO-3G basis (left), ethane with 6-31G basis (center), and carbon dioxide with the 6-31G basis (right).We show the third-order qSWIFT by the pink line (qSWIFT-3), the sixth-order qSWIFT by the pink dotted line (qSWIFT-6), and the qDRIFT by the black line (qDRIFT).For the Trotter-Suzuki decomposition, the best of the deterministic method among the first, second, and fourth order is shown by the gray dotted line (TS (Best)), and the best of the randomized method is shown by the green dotted line (RTS (Best)).FIG.4: The asymptotic behavior to achieve ε = 10 −6 .Other settings are the same as Fig. 3.
Fig 3 show the asymptotic behaviors of N for each time to achieve the systematic error ε = 0.001, which includes three subfigures corresponding to each molecule: propane with STO-3G basis (left), ethane with 6-31G basis (center), and carbon dioxide with 6-31G basis (right).To generate each molecule Hamiltonian, we use OpenFermion [16].For each figure, we show the third-order qSWIFT by the pink line (qSWIFT-3), the sixthorder qSWIFT by the pink dotted line (qSWIFT-6), and the qDRIFT by the black line (qDRIFT).For the Trotter-Suzuki decomposition, the best of the deterministic method among the first, second, and fourth order is shown by the gray dotted line (TS (Best)), and the best of the randomized method is shown by the green dotted line (RTS (Best)).
We see that the qSWIFT algorithms outperform the qDRIFT for all t in the sense that the required number of gates N is more than 10 times smaller in the qSWIFT algorithms than in the qDRIFT.While the original qDRIFT reduces the number of gates in the region t 10 8 from the Trotter-Suzuki decompositions, the qSWIFT algorithms realize a further reduction of the gates.Note that the improvement from the third-order qSWIFT to the sixth-order qSWIFT is not as large as that from the qDRIFT to the third-order qSWIFT.Since there is a trade-off between the order and the number of quantum circuits run as we discuss in Appendix B, it may be better to utilize the third-order qSWIFT rather than the sixth-order qSWIFT though it depends on the features of quantum devices.Fig. 4 is the same figure as Fig. 3 other than that, the required systematic error is ε = 10 −6 .In this case, where more precise time evolution is necessary, the merit of using qSWIFT is much clearer.In qSWIFT-3 (qSWIFT-6), the required number of gates for each time is almost 1,000 (10,000) times smaller than that of qDRIFT.The region where qDRIFT has an advantage over the Trotter-Suzuki is very limited (t 10 5 ∼ 10 6 ) due to the bad scaling of qDRIFT in terms of ε.In contrast, the region where qSWIFT has the advantage over Trotter-Suzuki decompositions does not change much from the case of ε = 0.001 (t 10 9 ∼ 10 10 ).This result shows the merit of our algorithm; the case of qDRIFT, we need to increase the number of gates to reduce the systematic error, but in the case of our qSWIFT, it can be reduced just by increasing the order parameter of the algorithm.FIG.5: The estimation error of Tr(QU(ρ init )) for each number of gates N for each method with Q = ZIIIIIII and ρ init = |+ ⊗8 .The time evolution is performed by the hydrogen molecule Hamiltonian with 6-31g basis transformed by the Bravyi-Kitaev transformation.The evolution time is set to be t = 1.In plotting each point, we perform six trials and show the mean value and the standard deviation of the mean.For qSWIFT, we show the result of the second-order (qSWIFT-2) with the purple line and the third-order (qSWIFT-3) with the pink line.The result of the qDRIFTalgorithm (qDRIFT) is plotted with the black line.For the deterministic Trotter-Suzuki decomposition, the first-order result (TS-1st) is plotted with the dotted gray line, and the second-order result (TS-2nd) is plotted with the dotted green line.For the randomized Trotter-Suzuki decomposition, the first-order result (RTS-1st) is plotted with the dotted yellow line, and the second-order result (RST-2nd) is plotted with the dotted blue line.
We estimate Tr(QU(ρ init )) by using qSWIFT algorithm described in Algorithm 2 with an observable Q, time evolution U, and an input state ρ init .For U, we implement the time evolution with the hydrogen molecule Hamiltonian with the 6-31G basis and t = 1.Again we use OpenFermion [16] to generate the molecule Hamiltonian.We utilize the Bravyi-Kitaev transformation [17] for transforming the Fermionic operators to Pauli operators.The number of terms in the Hamiltonian L is 184.The generated Hamiltonian has eight qubits, and therefore, we use a system with nine qubits (including one ancilla qubit).As for the observable Q, we choose Q = ZIIIIIII, and as for the input state, we choose ρ init = |+ ⊗8 .For comparison, we also estimate Tr(QU(ρ init )) by using the qDRIFT and the Trotter-Suzuki decomposition.As the input parameters of the qSWIFT, we set N 0 shot = N shot ( n) = 100, N 0 sample = 400, 000, and N sample ( n) = C( n) × N 0 sample .For qDRIFT and the randomized Trotter-Suzuki decomposition, we sample N 0 sample quantum circuits and perform N 0 shot measurements for each circuit.For the deterministic Trotter-Suzuki decomposition, we perform N 0 shot × N 0 sample measurements.For the quantum circuit simulation, we use Qulacs [18].Fig. 5 shows the estimation error of Tr(QU(ρ init )) for each number of gates N for each method.In plotting each point, we perform six trials and show the mean value and the standard deviation of the mean.For qSWIFT, we show the result of the second-order (qSWIFT-2) with the purple line and the third-order (qSWIFT-3) with the pink line.The result of the qDRIFT (qDRIFT) is plotted with the black line.For the Trotter-Suzuki decomposition, we show the first and second-order results.The minimum N for the first and second Trotter-Suzuki decomposition is 184 and 368 (1840 for the fourth order).For the deterministic Trotter-Suzuki decomposition, the first-order result (TS-1st) is plotted with the dotted gray line, and the second-order result (TS-2nd) is plotted with the dotted green line.For the randomized Trotter-Suzuki decomposition, the first-order result (RTS-1st) is plotted with the dotted yellow line, and the second-order result (RTS-2nd) is plotted with the dotted blue line.
We see that qSWIFT algorithms outperform the other methods.Particularly, the required N for the third-order qSWIFT for achieving ε ∼ 0.001 is almost 10 times smaller than that for qDRIFT, which is consistent with the asymptotic behavior shown in Fig. 3. Consequently, even in the region where Trotter-Suzuki decomposition works better than qDRIFT, qSWIFT algorithms outperform Trotter-Suzuki decompositions.

VI. CONCLUSION AND DISCUSSION
Hamiltonian simulation is a crucial subroutine of various quantum algorithms.Approaches based on the product formulae are practically favored due to their simplicity and ancilla-free nature.There are two representative product formulae-based methods: Trotter-Suzuki decompositions and qDRIFT [5].Trotter-Suzuki decompositions have the issue that the number of gates depends on the number of the terms in the Hamiltonian, at least linearly.In contrast, qDRIFT avoids the dependency on the number of terms but has the issue that the number of gates is dependent on the systematic error ε as O(1/ε).
In this paper, we propose qSWIFT, a high-order randomized algorithm having both the advantage of the Trotter Suzuki decompositions and qDRIFT, in the sense that its gate count is independent of the number of terms and that the gate count is asymptotically optimal with respect to the precision.We construct the qSWIFT channel and bound the systematic error by the diamond norm.We prove that the qSWIFT channel satisfies the required precision that decreases exponentially with the order parameter in terms of the diamond norm.Then we create the algorithm applying the qSWIFT channel to estimate given physical quantities.The algorithm requires a system as simple as qDRIFT; it requires just one ancilla qubit, and only the time evolution operators and the swift operators constructible with the gate in Fig. 2b are necessary to construct the quantum circuits.Our numerical demonstration shows that qSWIFT outperforms qDRIFT with respect to the number of gates for required precision.Particularly when high precision is required, there is a significant advantage of using qSWIFT; the number of gates in the third-order (sixth-order) qSWIFT is 1000 (10000) times smaller than qDRIFT to achieve the systematic error of ε = 10 −6 .
As a future direction, it is beneficial to perform case studies to investigate the performance of qSWIFT in specific problems.Particularly, the literature [19] points out that the qDRIFT does not perform well in the phase estimation problems, unlike originally expected [5] due to the relatively large systematic error for a given number of gates.In contrast, since qSWIFT can successfully reduce systematic error by increasing the order parameter, we expect that the performance in the phase estimation is significantly improved with qSWIFT.Also, as we note in Section I, there are algorithms using many ancilla qubits [6][7][8][9][10][11] that achieve a better asymptotic gate scaling with respect to λ and t though its constant prefactor is relatively large compared to qDRIFT and qSWIFT.Comparing qSWIFT with those algorithms and discussing the advantages of each algorithm in specific problems is a promising direction for future work.

e
a) A quantum circuit for qDRIFT with N segments.iH 1 τ • • • e iH r τ e iH r+2 τ • • • e iH N τ(b) A swift circuit with N segments, where two swift operators are performed in the (r + 1)-th segment.

( 2 ) 2 S(b1) 1
s ( ) sample = ( 1 , 2 ).With probability P mdrift ( k, r) and P mdrift ( k , N − 1 − r), sample k and k .3. Estimate Tr Q TN−1−r ( k ) S(b2) Tr ( k)(ρ init ) (65) with N shot measurements and set the resulting value to δ q(s, b 1 , b 2 ), which is an unbiased estimator of δq(s, b 1 , b 2 ).The value (65) can be evaluated by applying a swift circuit composed of the time evolution and the swift operators and estimating the expectation value of Q by measurements.We repeat the above process N sample times, and the estimate of δq(s, b 1 , b 2 ) is computed as the sample average.By substituting each estimate of δq(s, b 1 , b 2