NISQ Algorithm for Hamiltonian Simulation via Truncated Taylor Series

Simulating the dynamics of many-body quantum systems is believed to be one of the ﬁrst ﬁelds that quantum computers can show a quantum advantage over classical computers. Noisy intermediate-scale quantum (NISQ) algorithms aim at eﬀectively using the currently available quantum hardware. For quantum simulation, various types of NISQ algorithms have been proposed with individual advantages as well as challenges. In this work, we propose a new algorithm, truncated Taylor quantum simulator (TQS), that shares the advantages of existing algorithms and alleviates some of the shortcomings. Our algorithm does not have any classical-quantum feedback loop and bypasses the barren plateau problem by construction. The classical part in our hybrid quantum-classical algorithm corresponds to a quadratically constrained quadratic program (QCQP) with a single quadratic equality constraint, which admits a semideﬁnite relaxation. The QCQP based classical optimization was recently introduced as the classical step in quantum assisted eigensolver (QAE), a NISQ algorithm for the Hamiltonian ground state problem. Thus, our work provides a conceptual uniﬁcation between the NISQ algorithms for the Hamiltonian ground state problem and the Hamiltonian simulation. We recover diﬀerential equation-based NISQ algorithms for Hamiltonian simulation such as quantum assisted simulator (QAS) and variational quantum simulator (VQS) as particular cases of our algorithm. We test our algorithm on some toy examples on current cloud quantum computers. We also provide a systematic approach to improve the accuracy of our algorithm.


Introduction
Digital quantum computers have made immense progress in recent years, advancing to solving problems considered to take an unreasonable time to compute for classical computers [1,2]. In short, we are now in the Noisy Intermediate-Scale Quantum (NISQ) era [3,4], which is characterized by quantum computers with up to a few hundred noisy qubits and lacking full scale quantum error correction. Thus, noise will limit the usefulness of the computations carried out by these computers [3], preventing algorithms that offer quantum advantage for practical problems, such as Shor's algorithm for prime factorization [5], from being implemented.
However, just because such algorithms cannot be implemented on NISQ devices does not mean that quantum advantage for practical problems cannot be found with NISQ devices. There is currently great interest in the quantum computing and quantum information community to develop algorithms that can be run on NISQ devices but yet deal with problems that are practical [4,6,7]. Some of the most promising avenues deal with the problems in many-body physics and quantum chemistry. One major problem in this field is to develop algorithms capable of estimating the ground state and energy of many-body Hamiltonians. To such ends, algorithms like variational quantum eigensolver (VQE) [8,9] and quantum assisted eigensolver (QAE) [10,11] have been proposed.
The other major problem is to be able to simulate the dynamics of these many-body Hamiltonians. This task can be extremely challenging for classical computers, and Feynman proposed that this would be one of the areas where quantum computers could exhibit clear advantages over classical computers [12]. Powerful methods to simulate quantum dynamics on fault-tolerant quantum computers have been proposed, such as the concept of using truncated Taylor series by Berry et al [13].
On NISQ devices, a standard approach in simulating quantum dynamics is to utilize the Trotter-Suzuki decomposition of the unitary time evolution operator into small discrete steps. Each step is made up of efficiently implementable quantum gates, which can be run on the quantum computer [14][15][16][17][18][19][20]. However, the depth of the quantum circuit increases linearly with evolution time and the desired target accuracy. On NISQ devices, fidelity rapidly decreases after a few Trotter steps [21], implying long time scales will be unfeasible to simulate with this method. Alternatives to Trotterization have been proposed, such as VQS [22][23][24], subspace variational quantum simulator (SVQS) [25], variational fast forwarding (VFF) [26,27], fixed state variational fast forwarding (fsVFF) [28], quantum assisted simulator [29,30] and generalized quantum assisted simulator (GQAS) [31] to name a few.
Recently, Otten, Cortes and Gray have proposed the idea of restarting the dynamics after every timestep by approximating the wavefunction with a variational ansatz [32]. Building on that, Barison, Vicentini and Carleo have proposed a new algorithm [33] for simulating quantum dynamics. Their algorithm, named projected variational quantum dynamics (pVQD) combines the Trotterization and VQS approaches [22,23]. They replace the differential equation with an optimization problem, although not well characterized, and require much simpler circuits compared to VQS. However, pVQD requires a quantum-classical feedback loop and might suffer from the barren plateau problem [34] as well the optimization problem may be computationally hard [35]. Further, the feedback loop mandates that one has to wait for each computation to finish before the next computation is run, which can be a major bottleneck on cloud-based quantum computers that are accessed via a queue.
Here, we propose the truncated Taylor quantum simulator (TQS) as new algorithm to simulate quantum dynamics. Our algorithm is building on the ideas of pVQD [32,33] combined with the ansatz generation of QAS [29], which we further enhance by applying the concept of truncated Taylor series by Berry et al [13]. Our contributions and our algorithm are as following: 1. We recast the simulation of the quantum dynamics as a quadratically constrained quadratic program (QCQP). This optimization problem, unlike the optimization problem in pVQD, is well characterized and invites rigorous analysis. The QCQP in our algorithm admits a semidefinite relaxation [10]. Moreover, based on ideas from [10], one can provide a sufficient condition for a local minimum to be a global minimum, which a solver can further use as a stopping criterion. Since the classical optimization program in QAE is also a QCQP, it helps us achieve a conceptual unification of TQS with QAE.
2. The differential equations which form the classical part of QAS and VQS can be recovered in our framework. Since VQS is already a particular case of QAS [29], our approach yields both VQS and QAS as special cases of TQS.
3. We remove the need for the classical-quantum feedback loop in pVQD. The absence of the feedback loop yields our algorithm to be exceptionally faster than the feedback loop based NISQ algorithms for simulating quantum dynamics such as [22,[25][26][27][28].
4. Our algorithm avoids the trainability issues that plague other variational quantum algorithms. The choice of a problem-aware ansatz and the structure of the TQS algorithm helps bypass the barren plateau problem. It is known that in variational quantum algorithms that rely on a parametric quantum circuit, there will always be a tradeoff between trainability and expressibility, implying that a highly expressible ansatz cannot be easily trainable [36]. In our case, we do not rely on parametric quantum circuits, thus we bypass this problem. Furthermore, our algorithm provides a systematic way to obtain a more expressible ansatz, which is missing in other algorithms.

TQS Approach
Let us first assume that the Hamiltonian H is expressed as a linear combination of r tensored Pauli matrices with coefficients β i ∈ . The unitary evolution under the action of the Hamiltonian H for time ∆t is given by We do not need to implement the action of the unitary evolution in such a way. However, for purposes of describing the algorithm, we will use this power series expansion first, and talk more about alternatives later. We will now truncate the series, similar to [13]. If we choose small values of ∆t with respect to the eigen energies of H, we can approximate the unitary evolution with V (∆t) The classical evolution timestep ∆t should be chosen smaller than all relevant timescales of the Hamiltonian H to be simulated. This requires knowledge of the spectrum of H, which in general is not available. However, we can find appropriate values for ∆t in an heuristic manner. In our algorithm, the evolution with ∆t is performed on a classical computer only and thus we can choose any value for ∆t without requiring any additional quantum computational cost. Thus, we can simply evolve with a very small value for ∆t. To verify it is small enough, we can repeat the classical evolution for an even smaller value such as ∆t/2. If the results for both ∆t and ∆t/2 match, we can assume that ∆t provides sufficient accuracy. Let us next choose the ansatz at time t as linear combination of elements from cumulative K-moment states, K (refer to [29] for the formal definition). These states are defined in the same way as in [29] and will be constructed with the help of the given Hamiltonian, by essentially considering powers of the Hamiltonian. In terms of Pauli matrices, given a set of r tensored Pauli unitary matrices obtained from the unitary terms of the Hamiltonian P ≡ {P i } r i=1 and a positive integer K and some efficiently preparable quantum state |ψ〉, the K-moment states are the set of quantum states of the form for P i l ∈ P, where the indices i all run from 1 to r. We note that we only include unique states within the set {|χ〉} K . This corresponds to removing any repeated Pauli unitary in P. It should also be mentioned that the way the K-moment states are being generated is closely related to the Taylor expansion of the time evolution operator. If we consider the evolution of an arbitrary state by the time evolution operator, by observing that the Taylor expansion involves powers of the Hamiltonian H, it is clear that choosing the ansatz in such a way is suitable, as the |χ i 〉 ∈ {|χ〉} K states are essentially states in the Hilbert space of H K |ψ〉. This set is denoted by K . The cumulative K-moment states K are also defined in [29] as K ≡ ∪ K j=0 j . Now the ansatz is expressed as with some α i ∈ . For small values of ∆t, the ansatz at time t + ∆t is given by Using the ideas in [33], our goal now is to variationally approximate the time evolution of the system by adjusting our variational parameters. The crucial difference in our case is that our variational parameters α are coefficients which do not change the basis quantum states |χ i 〉. Thus, they can be solely updated via a classical computer and do not require a quantumclassical feedback loop. To evolve by time ∆t, we update the α i parameters to α i such that the following fidelity measure is maximized Using the notation |φ〉 = V (∆t) |ψ (α)〉 K , the expression for fidelity becomes Using the notation W φ ≡ |φ〉〈φ| 〈φ|φ〉 , the above expression further simplifies to The goal is to maximize the fidelity subject to the constraint that 〈ψ α |ψ α 〉 = 1. Thus, the optimization program at timestep t is given by Using the elements from K and the Hamiltonian H, we define the overlap matrices E and D as the following Because of the way the |χ n 〉 states are constructed, these values can be easily computed on a quantum computer, as they simplify to the expectation values of Pauli strings acting on the original quantum state |ψ〉. The constraint in the optimization program 12 can written in terms of α as We proceed to write the objective in the optimization program 12 in terms of the overlap matrices E and D. In first order, we can simplify the expression Further, using the notation G ≡ (E − ι∆tD) we find Using Eq.15,16,17 and the notation W α ≡ Gαα † G † α † Eα , the optimization program in 12 can be re-expressed in terms of overlap matrices as The aforementioned optimization program is a quadratically constrained quadratic program with a single equality constraint. As described in [10], this QCQP admits a direct convex SDP relaxation. Moreover, the results from [10] provide a sufficient condition for a local minimum to be a global minimum, which a solver can further use as a stopping criterion. Alternatively, the problem can be solved with the classic Rayleigh-Ritz procedure by finding the eigenvector associated with the largest eigenvalue λ of the generalized eigenvalue problem W α α = λEα . It can be shown that in the limit of small ∆t, TQS reduces to QAS (see Appendix C). This could potentially give us a way to obtain systematic higher-order corrections to the QAS matrix differential equation. Interestingly, this is a conceptual unification of the ground state problem (QAE) with the dynamics problem (QAS) in the quantum assisted framework. In QAE, finding the ground state and ground state energy of a Hamiltonian was formulated to become a QCQP. In TQS, the problem of simulating the dynamics is also given as a QCQP. This is conceptually satisfying as the problem of finding the dynamics is expressed as e −i t H |ψ〉, which is mathematically similar to using imaginary time evolution to finding the ground state via e −τH |ψ〉. The aforementioned connection is also one of the primary justifications for ansatz selection in [11]. We note that as alternative it is possible implement the unitary evolution operator U(∆t) directly instead of the Taylor series expansion of Eq.7, however this would require the usage of Hadamard tests (see Appendix D).
We want to emphasize again that the quantum computer is only required to measure the overlap matrices E and D at the start of the algorithm. No quantum-classical feedback loop for optimization is required. The only optimization steps required are performed solely on the classical computer with knowledge of the overlap matrices. The algorithm is as follows: 1. Choose an efficiently implementable initial state |ψ〉, then choose some K>0 and form the unique K-moment states |χ i 〉 to construct the ansatz.
2. With knowledge of the Hamiltonian H, calculate the overlap matrices E and D on the quantum computer. The job of the quantum computer is now done.
3. Choose a small ∆t with respect to the eigenvalues of H and evolve the state forward in time using a classical computer, by solving the optimization program 18 subject to the constraint 19.
If a higher fidelity for the simulation is desired, one can increase K to acquire an ansatz with a higher expressibility. The timestep ∆t could be increased by including higher order terms in the power series expansion of U(∆t) in our calculations (Described in Appendix E).

Results
We first use TQS to simulate a 2 qubit Heisenberg model We apply it to evolve an initial randomized 2 qubit state |ψ 2 〉. This initial state is generated by 5 layers of U 3 rotations and CNOT gates on the 2 qubits (see Appendix A). We ran the TQS algorithm on the 5-qubit quantum computer ibmq_rome, available through IBM Quantum Experience. We used error mitigation by calibrating the measurement errors and applying a filter obtained from that calibration on our data with the toolbox provided in Qiskit [37]. The results are shown in Fig.1. The evolution of the state under TQS reproduces the exact behavior very well for an ansatz with K = 1 moment states, even in the presence of the noise of a real quantum computer. Next, we apply TQS to simulate a 4 qubit XX chain model on a quantum computer Although this Hamiltonian is analytically solvable, we simulate this as a proof of principle. In Fig.2, we simulate (21) on ibmq_rome with an initial randomized 4 qubit state, generated by 5 layers of U 3 rotations and CNOT gates (see Appendix A). We run it for the K = 1 to K = 3 moment states. The evolution of the state under TQS again reproduces the exact behavior very well for the K = 3 case. We would like to point out that our algorithm can accurately simulate dynamics even for long time periods. The only errors that enter our algorithm are due to the ansatz being not expressible enough, and noise in the measurement of the matrix elements. Both type of errors affect only the initial conditions of the classical post-processing part. However, errors do not enter during the computation of the evolution itself as they are fully calculated on the classical computer. If we are able to obtain very accurate initial measurements for our matrix elements, and use an ansatz that fully captures the solution space, we believe that our algorithm in general will be able to simulate the dynamics accurately for long timescales.   Figure 3: Time evolution of TQS on a 8 qubit state, with Hamiltonian H 8 , simulated on a classical computer, with a random initial state. The initial state was generated with 3 successive layers of U 3 rotations with randomized parameters on each qubit, followed by CNOT/entangling gates. This is further described in Appendix A. a) Expectation value of 〈Z 1 〉. b) Fidelity of the state.
Next, we investigate in Fig.3 the transverse Ising model with 8 qubits by simulating TQS on a classical computer.
With an initial random state, we find that the evolution of the state reproduces the exact dynamics for the case of K = 3 moment expansion. Lastly, we compare TQS to pVQD for a 2 qubit transverse Ising model on a simulation. We consider the 2 qubit transverse Ising Hamiltonian We compare the algorithms with noisy simulators, where the noise models taken from the IBM Quantum Experience provider. The results are shown in Fig.4. While both TQS and pVQD show errors when simulating this Hamiltonian in the presence of noise, the expectation values for TQS are closer to the exact results most of the time. This is especially the case for the expectation value of 〈Z 1 〉. However, while the results might be argued to be somewhat similar, the resource requirements of both algorithms on the quantum computer are quite different. The TQS algorithm requires ≈ 30 circuits to be run, while the pVQD simulator requires well over 4000 circuits, which is a major effort to run on the IBM Quantum Experience. We note that to increase the simulation time for this example, no extra circuits are required with TQS as the algebra already has closed, whereas the number of circuits in pVQD scales linearly with simulation time. Furthermore, TQS performs well with circuit that are shallower compared to pVQD, which requires a circuit with 8 variational parameters. This behavior of TQS requiring less variational parameters to get a similar result seems to be consistent for the small models we tested, as other variational algorithms usually need an over-parameterized ansatz to perform well.
We also compare TQS to Trotterization on a noisy simulator for the same 2 qubit transverse Ising Hamiltonian. A simple Trotterization of the time evolution operator for this case is decomposed as The results are shown in Fig.5. As can be seen, even for a simple case such as this, due to the circuit lengths in Trotter increasing linearly with the time, the circuit lengths rapidly grow too long to obtain any meaningful results from the quantum computer. This is in contrast to TQS, which is able to capture the dynamics faithfully.
In Fig.6, we study our algorithm for up to thousands of qubits N . We use a Hamiltonian H = r i=1 P i that consists of r randomly chosen N -body Pauli strings P i = ⊗ N j=1 σ j , where σ j ∈ {I, X , Y, Z}. The cumulative K-moment states close at order K = r and yield the full ansatz space necessary to describe the dynamics exactly. We use the product state |0〉 N as initial state for the dynamics. This choice makes the dynamics tractable for classical computation. However, choosing an highly entangled initial state |ψ〉 would require a quantum computer to evaluate the overlaps. For such intractable states, our method provides a possible quantum advantage.

Discussion and Conclusion
The currently proposed NISQ algorithms face problems in scaling up to system sizes where classical computers cannot simulate the same systems, or in other words, to the point where we would see quantum advantage. For example, VQS/SVQS/pVQD require the use of a quantumclassical feedback loop, usually require complicated circuits, share similar problems as VQE like the barren plateau problem, and lack a systematic way to generate a parameterized ansatz. VFF and fsVFF also suffer from lacking a systematic way to generate the ansatz, usually require complicated circuits and have to run a quantum-classical feedback loop at the start. Further, the no fast-forwarding theorem suggests that not all Hamiltonians will be able to be accurately diagonalized with a reasonable amount of gates and circuit length, and the optimization step of the cost function in VFF might be too difficult to carry out efficiently. However, the barren plateau problem and ansatz state generation could be improved upon by applying various techniques [36,[38][39][40][41].
One problem that VQS and QAS share is that they require solving a differential equation which includes the pseudo-inverse of a matrix, whose elements are measured on a quantum computer. This matrix can be ill-conditioned. This procedure, via singular value decomposition, can be numerically unstable and sensitive to noise, especially as the system increases in size [42]. However, the sensitivity of these matrices has not been rigorously analyzed and more work has to be done to understand the scaling of the sensitivity.
In this work, we develop TQS for simulating quantum dynamics on digital quantum computers. TQS recasts the dynamical problem as a QCQP optimization program, which is well characterized unlike the optimization program in pVQD, allowing us to avoid the aforementioned problem in VQS and QAS.
At the same time, TQS retains the advantages of QAS, namely providing us a systematic method to select the ansatz, avoiding complicated Hadamard tests and controlled unitaries, avoiding the barren plateau problem, and only requiring usage of the quantum computer at the start, all of which are problems that are present in pVQD.
However, there are still many problems to tackle in our approach. One problem is an inherited problem from QAS. As the Hamiltonian size and complexity increase, large K values may be needed to generate enough states for a sufficiently expressible ansatz to produce accurate results. It is clear from the connection between the Taylor expansion of the time evolution operator and our K-moment states that in the general case, the further in time we want to simulate, the exponentially larger our ansatz should be and the harder the difficulty of generating that ansatz. However, this is fundamentally a complexity theoretic statement which can not be bypassed in the general case by any quantum simulation algorithm based on parametric quantum circuits (variational quantum algorithms) or linear combination of quantum states (our algorithm). This problem particularly emerges in variational algorithms for time evolution. For example, in algorithms such as VQE for finding the ground state of Hamiltonians, we know that the ground state of locally gapped Hamiltonians fulfil area laws of entanglement and thus do not need exponentially many parameters to be described. However, for the time evolution over longer times a similar statement about the complexity of the problem is not known. Though our algorithm uses a problem aware ansatz, more information from the problem such as the combination coefficients β i and symmetries of the Hamiltonian could be employed to further tame the complexity. A further discussion and analysis on the number of states needed is given in Appendix B.
As the system size increases, it may be required to reduce ∆t to preserve accuracy in the classical post-processing part of the algorithm. This will increase the computational cost of the classical computer, however it requires no additional quantum computations. The number of classical optimization steps to be carried out increases linearly with the number of discretiza-tions steps of the evolution time. Determining whether this poses a bottleneck for TQS when applied to large systems requires further studies.
Furthermore, in the presence of noise, the calculated fidelity of our states can go above one. A possible origin are small eigenvalues in the E overlap matrix, which can give the procedure of optimizing or solving the generalized eigenvalue problem numerical instability. As we scale up the system and consider more ansatz states, this issue can become more prevalent.
We expect our algorithm not to provide quantum advantage in the general case. However, we believe our algorithm is capable of providing quantum advantage over classical methods for certain cases. The conditions where we believe our algorithm will do so are: • The basis states which are used to represent the initial quantum state are highly entangled such that they cannot be stored on a classical computer. This will render the calculation of corresponding overlaps classically hard, as it boils down to a circuit sampling task. Note that the Quantum Threshold Assumption (QUATH) by Aaronson and Chen [43] says that there is no polynomial-time classical algorithm which takes as input a random circuit C and can decide with success probability at least 1 2 + Ω 1 2 n whether |〈0 n |C|0 n 〉| 2 is greater than the median of |〈0 n |C|x n 〉| 2 taken over all bit strings x n . In other words, the circuit sampling task is difficult and hence classical algorithms will not be able to compete with algorithms based on circuit sampling as system size scales. The quantum part of TQS is based on circuit sampling which is classically difficult.
• The Hamiltonian possesses a particular structure. For example, the Hamiltonian consists of a small number of unitaries, the Krylov subspace closes fast, or the Hamiltonian is a low-rank matrix. We demonstrated such an example for a Hamiltonian consisting of a limited amount of multi-body Pauli strings where our method can simulate the dynamics of thousands of qubits. These Hamiltonians would be challenging for other methods such as Trotter or variational quantum algorithms. For those algorithms, the multi-body interactions and the large number of qubits would require an extensive number of gates and circuit depth to accurately represent the evolved state. A further example where our algorithm can perform well are quantum many-body scars. This quantum many-body phenomena can arise when the Krylov subspace closes fast at a low order K [44], which is exactly the condition needed for our algorithm to perform well. The timescales that can be reasonably approximated by our algorithm is dependent on the Hamiltonian in question. Arbitrary Hamiltonians without the aforementioned conditions explore the full Hilbertspace during the evolution. Thus, it will be difficult for our ansatz to cover the whole solution space and approximate the dynamics accurately. Note that other variational quantum algorithms suffer similar problems as their ansatz is restricted to polynomial number of parameters. In the case of general Hamiltonians, our algorithm can provide systematic approximations for the quantum evolution of short time scales.
• The system size of interest and the amount of entanglement of the quantum state should be beyond the reach of classical simulation methods. Here, our algorithm can make use of the power of the quantum computer to prepare and measure classically intractable states.
In the future, the NISQ community should investigate these challenges, so that we can successfully run NISQ algorithms for larger qubit numbers.

Acknowledgments
We are grateful to the National Research Foundation and the Ministry of Education, Singapore for financial support. The authors acknowledge the use of the IBM Quantum Experience devices for this work. This work is supported by a Samsung GRC project and the UK Hub in Quantum Computing and Simulation, part of the UK National Quantum Technologies Programme with funding from UKRI EPSRC grant EP/T001062/1.

A Details on running circuits on the IBM quantum computer
For the runs on the real quantum computer, we generated an initial state with randomized parameters to evolve with the following circuit. It comprised 5 layers of successive U 3 rotation with randomized parameters on each qubit, followed by a CNOT/entangling gate (see Fig.7 and 8). We sampled from each circuit with 8192 shots.

|0〉
R Figure 7: Circuit for two qubits that generate one set of U 3 rotation with randomized parameters, followed by a CNOT gate between the 2 qubits. 5 successive layers of this circuit were used to generate the initial starting state for the 2 qubit case on the IBM quantum computer for our runs of TQS. The Θs were randomly generated.
|0〉 Figure 8: Circuit for four qubits that generate one set of U 3 rotation with randomized parameters, followed by a series of CNOT gates between the adjacent qubits. 5 successive layers of this circuit were used to generate the initial starting state for the 4 qubit case on the IBM quantum computer for our runs of TQS. The Θs were randomly generated.

B Number of basis states considered for each K, and discussion on scaling
The number of basis states that was used to construct the hybrid ansatz, for each K moment expansion, for each Hamiltonian, is given in Table 1. Given a scalar τ, an N × N matrix A and an N × 1 vector v, the action of the matrix exponential operator exp (τA) on v can be approximated as where p K−1 is a K − 1 degree polynomial. The approximation in equation 25 is an element of the Krylov subspace, Thus, the problem of approximating exp (τA) v can be recast as finding an element from K r K . Note that the approximation in equation 25 becomes exact when K −1 = r ank(A). In our case, we can identify v with the initial state |ψ〉, τ with −ι t and A with the Hamiltonian H.
In the worst case, the number of overlaps scales as O(r K ) for r terms in H. By observing the Taylor expansion of the time evolution operator exp(−iH∆t), we can see that at longer times we would struggle with finding an expressible enough ansatz in the general case, as we need to keep considering higher powers of H. This is fundamentally an expressibility problem, present in all NISQ variational algorithms, be it based on linear combination of states or those based on parametric quantum circuits. It is known that to prepare an arbitrary state on an n qubit quantum computer, we require a circuit depth of at least 1 n 2 n [45][46][47][48]. This suggests that it is very hard to produce an expressible enough Ansatz to reproduce an arbitrary quantum state in the Hilbert space.
It is known that the the Krylov subspace spans the entire space when you exponentiate the Hamiltonian H to the power of K − 1, where K − 1 = r ank(H). Thus, the number of states that we require in our Ansatz scales linearly with the rank of the Hamiltonian.
Furthermore, one of the major contributions of the TQS algorithm is that, by using this problem-aware Ansatz, it provides a systematic way to obtain a more and more expressible Ansatz. The other variational algorithms like VQS and VFF still do not have a systematic method to generate an expressible enough Ansatz, or to improve on an Ansatz in a efficient way. Also, it has been shown that if we use a hardware efficient Ansatz, we would in general expect to encounter the barren plateau problem, which makes it very hard for the algorithm to train and optimize [34,49]. Furthermore, the usual technique of using more and more layers of hardware efficient Ansatz circuits gives no guarantee that it will become more and more expressible in an efficient manner, when compared to the number of variational parameters that we are adding. There is also no guarantee that this will indeed improve the appropriateness of the Ansatz. This is especially true for larger systems. In TQS, with the way we generate the Ansatz with K moment states, we can see that at worst, we get an Ansatz with as many states as the size of the Hilbert space, which is fully expressible. This is due to the group of Pauli strings closing on itself eventually. Also, we can see that as we increase the K, we will definitely improve our Ansatz and get to a point where it is eventually expressible enough. In future, using the coefficients of the terms in the Hamiltonian, we expect to be able to slow down the growth of the number of states.
Our algorithm also relies on being able to calculate expectation values of powers of the Hamiltonian, 〈ψ|H k |ψ〉 in an efficient manner. If we look at the Pauli string level (break our Hamiltonian into linear sums of Pauli strings), the number of Pauli terms in H k grows exponentially in k. Right now, for current implementation of our algorithm on available quantum computers, this breaking into Pauli strings is necessary due to the imperfections in said quantum computers. However, if we allow more complex operations that cannot be performed very well right now, such as complex controlled unitaries, the resources needed to measure such 〈ψ|H k |ψ〉 values might scale less [50].
We would also like to mention that depending on the Lie algebra of the Pauli terms in the Hamiltonian and the rank of the Hamiltonian, the number of required overlaps can be a lot smaller compared to the upper bound. By considering specific kinds of Hamiltonians, the number of states needed will be manageable. As an example, for a system size with a multiple of 3 qubits, if we consider the Hamiltonian of the form H = X Y Z X Y Z...X Y Z + Y Z X Y Z X ...Y Z X + Z X Y Z X Y...Z X Y + X X X X X ...X X X , the set of Kmoment states is maximally size 8, implying that 8 ansatz states are sufficient to simulate the dynamics with our algorithm.

C QAS and VQS as special cases of TQS
In this appendix, we show that in the limit of choosing a very small ∆t, one obtains QAS from TQS. Since VQS is a special case of QAS [29], we get VQS also as special case of TQS. We start out with the series expansion of |ψ( α + δ α)〉 |ψ( α + δ α)〉 = |ψ( α)〉 + j ∂ ∂ α j |ψ( α)〉 δα j .
Now in the same manner as QAS, using the Mclachlan's variational principle [23,29,30,51], we demand that the variation of this fidelity is equal to 0 with respect to α j : Now we substitute in U(δt) = I − i∆t H: Now we take the derivative of this equation with respect to ∆t. Note that d d∆t δα * k = δα * k . We then discard any terms remaining that are linear in ∆t or in δα (implying we have chosen such a small ∆t that δα is also very small).
Using the above definition of the E and D matrices in equation 13 and 14, this simplifies to: