Advances in time-independent Hamiltonian simulation algorithms

Abstract

Quantum simulation is a rapidly advancing field poised to revolutionize our understanding of complex quantum systems by harnessing the unique capabilities of quantum computers. In this review, we present a concise overview of key developments in quantum simulation algorithms, with a focus on time-independent digital Hamiltonian simulation of quantum dynamics. We review seminal methods—from traditional Trotter formulas to techniques such as truncated Taylor series combined with the linear combination of unitaries, and quantum signal processing. Additionally, we examine error analyses and recent innovations that enhance the efficiency and accuracy of quantum simulations.

Graphical abstract

Keywords

hamiltonian simulation / digital quantum simulation / product formulas / quantum signal processing / quantum singular value transformation

1 Introduction

The concept of quantum simulation was first proposed by Richard Feynman in 1982, who suggested that quantum systems could be efficiently simulated using quantum computers [1]. As Feynman put it:

“Nature isn’t classical, dammit, and if you want to make a simulation of nature, you’d better make it quantum mechanical, and by golly it’s a wonderful problem, because it doesn’t look so easy.”

This idea laid the groundwork for the development of various quantum algorithms designed to simulate the dynamics of quantum systems. Quantum dynamical simulation problems, such as Hamiltonian simulation, can be formulated as follows:

Problem 1 (Time-independent Hamiltonian simulation). Given a Hamiltonian H=\sum _{l=1}^L{H}_l that is independent of evolution time t (i.e., H(t)=H), Hamiltonian simulation refers to approximating the evolution operator U=e^{-\mathrm{i}Ht} by an implemented procedure \mathcal{A}. If \mathcal{A} is a deterministic unitary circuit implementing U^{\prime }, we quantify the simulation error by the operator norm \Vert U-U^{\prime }\Vert \le ϵ. If \mathcal{A} realizes an effective quantum channel \mathcal{E} of evolution, we quantify the error at the channel level, typically by the diamond norm \Vert \mathcal{E}-\mathcal{U}{\Vert }_{\diamond }, where \mathcal{U}(\rho )=U\rho U^{†}. Specifically, the goal is to determine an upper bound on the number of quantum gates required for this simulation.

The first practical quantum algorithm for Hamiltonian simulation was introduced by Seth Lloyd in 1996 [2]. Lloyd demonstrated that any local quantum system could be simulated efficiently within the quantum circuit model. Following Lloyd’s work, significant progress has been made in enhancing the efficiency and accuracy of Hamiltonian simulations. One major development was the introduction of higher-order product formulas [3], such as the Suzuki-Trotter decomposition [4], which provided improved precision by reducing the error scaling with the number of time steps. These methods allowed for more accurate simulations with fewer computational resources.

In addition to product formulas, the development of algorithms based on quantum walks and linear combinations of unitaries (LCU) [5,6] has expanded the toolkit for post-Trotter Hamiltonian simulation methods [7] (see Fig. 1(b)). One of the prominent post-Trotter methods is the truncated Taylor series (TS) method [5], which employs the LCU technique to implement Hamiltonian simulation. This approach approximates the time evolution operator by truncating the series expansion and has proven particularly useful for simulating systems with complex interactions and long-range correlations. Furthermore, quantum signal processing (QSP) and quantum singular value transformation (QSVT) represent significant advances in post-Trotter Hamiltonian simulation [8,9]. QSP leverages polynomial approximations and quantum control techniques to simulate Hamiltonian dynamics with optimal parameters. QSP has been recognized for its ability to unify many quantum algorithms through eigenvalue transformation, making it a promising tool for large-scale quantum computing [10]. QSVT takes a further step by enabling singular value transformation on non-Hermitian and non-square matrices.

The advancements in Hamiltonian simulations have not only broadened the scope of quantum simulation algorithms but also paved the way for practical applications in various scientific domains. From quantum chemistry and materials science to condensed matter physics [11–13], the ability to simulate intricate quantum systems with high fidelity is advancing our understanding and exploration of quantum phenomena.

In this review, we survey recent advancements in time-independent Hamiltonian simulation algorithms, covering product formulas (PF), truncated Taylor series (TS) methods, quantum signal processing (QSP) and qubitization, and quantum singular value transformation (QSVT). We first establish the basic problem formulation, access models, and complexity measures in Section 2. Section 3 reviews the main advances in PF-based simulation, followed by a dedicated discussion of PF variants in Section 4. Section 5 is devoted to truncated TS methods. In Section 6, we review Hamiltonian simulation based on QSP and qubitization, and in Section 7 we present the QSVT-based framework and its implications for Hamiltonian simulation. In Section 8, we discuss rigorous theoretical advantages and fundamental limitations of Hamiltonian simulation. Finally, in Section 9, we conclude with discussion and outlook.

2 Preliminaries

2.1 Notations

In the following context, the operator norm of the operator A is denoted as \Vert A\Vert. For A=\sum _j{\alpha }_j{A}_j, the 1-norm is defined as {\Vert A\Vert }_1=\mathrm{T}\mathrm{r}\sqrt{A^{†}A}. The normalized Frobenius norm is defined as {\Vert A\Vert }_F=\sqrt{\mathrm{T}\mathrm{r}(A^{†}A)/d}, where d is the Hilbert-space dimension. For a vector |v⟩, the l_2 norm is {\Vert |v⟩\Vert }_{l_2}=\sqrt{⟨v|v⟩}. For a bipartite pure state {|\psi ⟩}_{AB}, the von Neumann (entanglement) entropy is S({|\psi ⟩}_{AB})= -\mathrm{T}\mathrm{r}({\rho }_A\ln {\rho }_A)=-\mathrm{T}\mathrm{r}({\rho }_B\ln {\rho }_B), where {\rho }_{A/B} is the reduced density matrix of subsystem A/B.

We also use standard channel distances when discussing randomized or non-deterministic simulation procedures [14]. For a linear map \mathrm{\Phi }, the diamond norm is defined as

\Vert \mathrm{\Phi }{\Vert }_{\diamond }:=\underset{k\ge 1}{sup}\underset{\rho :\Vert \rho {\Vert }_1=1}{sup}\Vert (\mathrm{\Phi }\otimes {\mathbb{I}}_k)(\rho ){\Vert }_1,

where {\mathbb{I}}_k is the identity channel acting on an auxiliary space of dimension k.

For block-encodings and post-selected constructions, it is also convenient to measure the partial operator norm of an operator A. Let |G⟩{⟨G|}_a\otimes \mathbb{I} denote the success subspace projector. For a unitary A on ancilla plus system intended to encode a system operator U on the success block, we define the partial operator norm as

\Vert A-U{\Vert }_{\mathrm{p}\mathrm{a}\mathrm{r}}:=\Vert ({⟨G|}_a\otimes \mathbb{I})A({|G⟩}_a\otimes \mathbb{I})-U\Vert .

The operator (spectral) norm \Vert U-U^{\prime }\Vert is the standard worst-case metric for unitary Hamiltonian simulation, as it directly upper-bounds the state-vector error for every input state, i.e., \Vert (U-U^{\prime })|\psi ⟩{\Vert }_{{\ell }_2}\le \Vert U-U^{\prime }\Vert for all normalized |\psi ⟩. However, several influential Hamiltonian-simulation frameworks are inherently channel-valued: randomized compilations (e.g., randomized product formulas) output a classical mixture of unitaries and therefore implement a CPTP map \mathcal{E}, for which the natural worst-case metric is the diamond norm \Vert \mathcal{E}-\mathcal{U}{\Vert }_{\diamond } or the 1-norm distance between the generated mixed state \Vert \mathcal{E}(\rho )-\mathcal{U}(\rho ){\Vert }_1. Likewise, block-encoding based methods may control errors in a partial operator norm on the success subspace of an ancilla measurement; when one conditions on successful post-selection, the induced operation on the system can again be compared to U=e^{-\mathrm{i}Ht} in the operator norm.

Since we will use the product of operators extensively, we denote the products with different directions by the arrows. For any operators \{A_l{\}}_{l=1}^m, we have

\overleftarrow{\prod _{l=1}^m}{A}_l=A_m\cdots A_2\cdot A_1,\overrightarrow{\prod _{l=1}^m}{A}_l=A_1\cdot A_2\cdots A_m.

To facilitate a rigorous comparison across the diverse Hamiltonian simulation algorithms discussed in this review, we establish the operator norm \Vert U-U^{\prime }\Vert as the fundamental worst-case error metric. While algorithms are often analyzed in their native error metrics depending on their structure (e.g., block-encodings, randomized channels, or observable tracking), these metrics can all be systematically related to operator-norm quantities. Table 1 summarizes these quantitative relationships.

2.2 Input models

To rigorously analyze the complexity of Hamiltonian simulation algorithms, we must explicitly define how the target Hamiltonian H is accessed. Different algorithms operate under different input assumptions.

2.2.1 Input model for Trotter product formula

We first define the input model most commonly used for product formulas and randomized compilation. We assume the target Hamiltonian H is provided as a decomposition into L summands, H=\sum _{l=1}^L{\alpha }_l{H}_l, where each H_l is a Hermitian operator with unit spectral norm. In the context of standard digital simulation, these terms are typically weighted Pauli strings P_l\in \{I,X,Y,Z{\}}^{\otimes n}.

The fundamental algorithmic resource in this model is the elementary exponential e^{-i\theta H_l}. For deterministic product formulas, complexity is quantified by the total number of such exponentials, referred to as the Trotter step count N, required to bound the simulation error by ϵ. To translate this algorithmic count into physical circuit depth, one assumes a compilation model where each Pauli exponential e^{-i\theta P_l} is synthesized using a sequence of Clifford and T gates.

Randomized compilation schemes, such as qDRIFT [14], operate under a similar input model but utilize the exponentials differently. The resource cost is defined by the number of sampled terms, each with an arbitrary evolution time, required to approximate the ideal quantum channel.

2.2.2 Sparse matrix oracle model

For algorithms designed for general sparse Hamiltonians (where the Pauli decomposition may be inefficient) and for establishing asymptotic query lower bounds, we use the Sparse Matrix Oracle model [19]. We assume H is a d-sparse matrix, meaning there are at most d non-zero entries in any row or column and each entry is a m-bit binary number. Access is provided via two unitary oracles:

● O_H: For row index j and column index k, it computes the matrix element H_{jk}:

O_H|j,k⟩|z⟩=|j,k⟩|z\oplus H_{jk}⟩.

● O_F: To efficiently locate non-zero elements, this oracle returns the column index of the l-th non-zero element in row j:

O_F|j,l⟩=|j,f(j,l)⟩,

where l\in \{1,\dots ,d\} iterates over the non-zero entries.

Constructing a block-encoding for a sparse Hamiltonian is typically achieved via two state-preparation unitaries, U_{row} and U_{col}, to prepare superpositions of row and column indices weighted by the Hamiltonian entries.

2.2.3 Block encoding via linear combination of unitaries

Following Berry et al. [5], when the Hamiltonian is provided as a linear combination of unitaries, H=\sum _l{\alpha }_l{H}_l, where {\alpha }_l>0 and the H_l are unitary operators, such as Pauli strings, the input access is defined by the ability to perform two specific oracle operations: PREPARE (G) and SELECT.

To implement the LCU algorithm or construct a block-encoding for such Hamiltonians, we define the two input oracles as follows:

\begin{array}{rl}& G|0⟩:=\frac{1}{\sqrt{\alpha }}\sum _l\sqrt{{\alpha }_l}|l⟩,\\ & \mathrm{select}(U):=\sum _l|l⟩⟨l|\otimes H_l,\end{array}

where \alpha =\sum _l{\alpha }_l is the normalization constant. In this model, complexity is quantified by the number of queries to G and \mathrm{select}(U). These oracles are combined to construct the walk operator W=(G^{†}\otimes \mathbb{I})\mathrm{select}(U)(G\otimes \mathbb{I}), which block-encodes H/\alpha.

Regardless of whether the input is accessed via LCU oracles or sparse-matrix oracles, the objective is to construct a unitary U in the standard form of a block-encoding:

U=\left(\begin{array}{cc}H/\alpha & \cdot \\ \cdot & \cdot \end{array}\right),

such that (⟨0{|}_a\otimes \mathbb{I})U(|0{⟩}_a\otimes \mathbb{I})=H/\alpha. Once H is encoded in this standard form, we can proceed with algorithm-agnostic techniques (such as QSP) to implement functions of H, independent of the underlying input structure. In these frameworks, the algorithmic complexity is quantified by the number of queries to U and its inverse/controlled versions.

2.2.4 Projected unitary encoding (PUE)

Following Gilyén et al. [9], many modern singular-value-based techniques are most naturally stated in a model that is slightly more general than the standard ancilla-|0⟩ block-encoding. A projected unitary encoding (PUE) of an operator A is specified by a triple (U,\tilde{\mathrm{\Pi }},\mathrm{\Pi }), where U is a unitary on an enlarged Hilbert space and \mathrm{\Pi },\tilde{\mathrm{\Pi }} are orthogonal projectors onto possibly different subspaces.

A:=\tilde{\mathrm{\Pi }}U\mathrm{\Pi },

which should be viewed as the action of U from the \mathrm{\Pi }-subspace into the \tilde{\mathrm{\Pi }}-subspace. In general, A can be non-Hermitian and even rectangular when \mathrm{r}\mathrm{a}\mathrm{n}\mathrm{k}(\mathrm{\Pi })\ne \mathrm{r}\mathrm{a}\mathrm{n}\mathrm{k}(\tilde{\mathrm{\Pi }}); moreover, A is automatically a contraction in the operator norm because it is obtained by projecting a unitary.

The key conceptual distinction from a standard block-encoding is that block-encoding fixes the projectors to a particular ancilla state, whereas PUE allows \mathrm{\Pi } and \tilde{\mathrm{\Pi }} to be arbitrary efficiently checkable subspaces. Concretely, an (\alpha ,a,\epsilon )-block-encoding of an operator A is recovered as the special case of PUE obtained by choosing \mathrm{\Pi }=\tilde{\mathrm{\Pi }}=|0⟩{⟨0|}^{\otimes a}\otimes \mathbb{I}, and taking A\approx \alpha ({⟨0|}^{\otimes a}\otimes \mathbb{I})U({|0⟩}^{\otimes a}\otimes \mathbb{I}) up to error \epsilon. Thus, projected unitary encoding strictly generalizes block-encoding.

Operationally, the PUE model assumes coherent access not only to U (and typically U^{†}), but also to controlled operations that discriminate the success subspaces defined by \mathrm{\Pi } and \tilde{\mathrm{\Pi }}. A convenient primitive is the controlled-\mathrm{\Pi } NOT,

C_{\mathrm{\Pi }}\mathrm{N}\mathrm{O}\mathrm{T}:=X\otimes \mathrm{\Pi }+\mathbb{I}\otimes (\mathbb{I}-\mathrm{\Pi }),

which flips an ancilla qubit if and only if the data register lies in \mathrm{i}\mathrm{m}\mathrm{g}(\mathrm{\Pi }) (and analogously for C_{\tilde{\mathrm{\Pi }}}\mathrm{N}\mathrm{O}\mathrm{T}). Using this primitive, one can implement phase rotations conditioned on the success subspace, e.g.,

e^{i\varphi (2\mathrm{\Pi }-\mathbb{I})}=C_{\mathrm{\Pi }}\mathrm{N}\mathrm{O}\mathrm{T}(\mathbb{I}\otimes e^{-i\varphi Z})C_{\mathrm{\Pi }}\mathrm{N}\mathrm{O}\mathrm{T},

and similarly for \tilde{\mathrm{\Pi }}. These conditional phases, interleaved with applications of U and U^{†}, form the alternating sequences used by singular value transformation: they manipulate the singular values of A=\tilde{\mathrm{\Pi }}U\mathrm{\Pi } while preserving its left and right singular subspaces.

From a complexity standpoint, PUE therefore measures cost in terms of queries to U (and U^{†}, and controlled versions when needed) together with queries to the subspace-checking primitives for \mathrm{\Pi } and \tilde{\mathrm{\Pi }}.

While product formula algorithms natively utilize the direct Hamiltonian summand input model without requiring any overarching encoding, all post-Trotter approaches can be conceptually and mathematically unified under the block-encoding framework. By treating block-encoding as the standard intermediate representation, techniques like Quantum Signal Processing (QSP) can be generically applied to any Hamiltonian regardless of its original physical access model. Table 2 outlines how various input models are translated into a standard block-encoding.

3 Product formula

Product-formula-based approaches have been fundamental to quantum simulation since Lloyd’s pioneering work [2], which first demonstrated explicit digital quantum circuits for simulating quantum dynamics. Lloyd’s method introduced the concept of decomposing a Hamiltonian into simpler sub-Hamiltonians and approximating the overall evolution by concatenating their individual evolutions, a construction now known as a product formula (see Fig. 1(b)). While initially restricted to constant local Hamiltonians, this approach was significantly expanded by Aharonov and Ta-Shma [20] to encompass arbitrary row-sparse Hamiltonians, broadening its practical applicability.

The development of product formula techniques has been marked by several theoretical breakthroughs. Childs [21] and Berry et al. [22] extended the approach beyond first-order Lie-Trotter approximations to incorporate higher-order Trotter-Suzuki decompositions [23], demonstrating nearly linear gate complexity with respect to simulation time. Somma [24] further improved these bounds by analyzing the role of commutation relations among sub-Hamiltonians. This commutation-based framework ultimately led to Childs et al.’s work [3,25], which established near-optimal scaling for geometrically local Hamiltonians and provided comprehensive error analyses across common Hamiltonian families.

In this section, we begin by presenting the standard product formula through Trotter-Suzuki decomposition in a self-contained manner. The general expression of the product formula enables us to extend the intuitive insights gained from low-order cases to more complicated scenarios. We will also discuss efforts aimed at accurately analyzing simulation errors associated with the Trotter-Suzuki formula. Additionally, we review prominent studies that have advanced the analysis of simulation errors by incorporating additional structural information. These advancements establish product-formula simulation as the most accessible simulation approach in the near term.

3.1 Input model and complexity metrics

To rigorously analyze the complexity of product formulas, we first define the input model and resource metrics. We assume the target Hamiltonian H is specified by a decomposition into L summands, H=\sum _{l=1}^L{H}_l. In the standard digital simulation model, each term H_l is assumed to be an efficiently implementable Hermitian operator, typically a weighted Pauli string h_l{P}_l where P_l\in \{I,X,Y,Z{\}}^{\otimes n}.

The fundamental resource in this framework is the elementary exponential e^{-\mathrm{i}{P}_l\tau } (a Pauli rotation). We quantify the algorithmic complexity by the total number of such Pauli evolutions—often referred to as the Trotter step count or query complexity—required to bound the simulation error by ϵ for a total evolution time t. This metric treats each term-wise evolution as a primitive operation.

To translate this algorithmic count into physical circuit depth, one must explicitly model the synthesis of each exponential. For a k-local Pauli term, the unitary e^{-\mathrm{i}{P}_l\tau } is synthesized using a sequence of single-qubit Clifford gates and CNOTs, followed by a single arbitrary-angle rotation. While the specific depth of this synthesis depends on the Pauli weight and hardware connectivity, it appears as a multiplicative constant in asymptotic scaling.

3.2 Product formula with Trotter-Suzuki decomposition

Given a Hamiltonian H=\sum _{l=1}^L{H}_l with L terms, one of the most straightforward ways to realize its time evolution is to interleave the individual sub-evolutions. Specifically, when the evolution time t is small, we obtain the following Lie-Trotter approximation:

(1)e^{-\mathrm{i}Ht}\approx e^{-\mathrm{i}{H}_{L}t}\cdots e^{-\mathrm{i}{H}_{2}t}{e}^{-\mathrm{i}{H}_{1}t}=\overleftarrow{\prod _l}{e}^{-\mathrm{i}{H}_{l}t},

with a second-order error:

\Vert \overleftarrow{\prod _l}{e}^{-\mathrm{i}{H}_{l}t}-e^{-\mathrm{i}Ht}\Vert =\mathcal{O}\left(\Vert H{\Vert }^2{t}^2\right).

Based on this, a natural extension is to interleave the sub-evolutions even more finely:

\begin{array}{rl}{e}^{-\mathrm{i}Ht}& \approx e^{-\mathrm{i}{H}_{1}t/2}\cdots e^{-\mathrm{i}{H}_{L}t/2}{e}^{-\mathrm{i}{H}_{L}t/2}\cdots e^{-\mathrm{i}{H}_{1}t/2}\\ & =\overrightarrow{\prod _l}{e}^{-\mathrm{i}{H}_{l}t/2}\overleftarrow{\prod _l}{e}^{-\mathrm{i}{H}_{l}t/2}.\end{array}

This back-and-forth structure improves the simulation error to third order, \mathcal{O}\left(\Vert H{\Vert }^3{t}^3\right).

This idea of using fine-grained mixing to improve the simulation accuracy is valid for even higher-order approximations. To illustrate the specific construction, we first introduce the product formula, which describes the most general concatenation of sub-Hamiltonian terms in product form [3]:

\mathcal{U}(t)=\overleftarrow{\prod _{y=1}^{\mathrm{{\rm Y}}}}\overleftarrow{\prod _{l=1}^L}{e}^{-\mathrm{i}t{a}_{(y,l)}{H}_{{\pi }_y(l)}},

where we have \mathrm{{\rm Y}} rounds of stages in the formula, \pi denotes the permutation of the sub-Hamiltonian terms, and the coefficients \{a_{(y,l)}\} control the corresponding evolutions. For example, the first-order approximation in Eq. (1) uses the constant coefficient 1, the identity permutation, and \mathrm{{\rm Y}}=1. The second-order approximation uses constant coefficient {\displaystyle \frac{1}{2}} with a pair of back-and-forth permutations, and \mathrm{{\rm Y}}=2. We denote these two formulas by {\mathcal{U}}_1(t) and {\mathcal{U}}_2(t), respectively.

To further improve the simulation accuracy, we specialize the generic formula to the Trotter-Suzuki form [23]. Typically, given a positive even integer p, the corresponding pth-order Trotter-Suzuki formula {\mathcal{U}}_p(t) with leading error \mathcal{O}\left((\Vert H\Vert t{)}^{p+1}\right) is constructed as

(2){\mathcal{U}}_p(t)={\mathcal{U}}_{p-2}(u_{p}t{)}^2{\mathcal{U}}_{p-2}(1-4u_{p}t){\mathcal{U}}_{p-2}(u_{p}t{)}^2,

where u_p=1/(4-4^{1/(p-1)}). In the context of the product formula, this pth-order Suzuki formula requires \mathrm{{\rm Y}}=2\times 5^{p/2-1} with an arbitrary pair of back-and-forth permutations along the round index y. Its coefficients \{a_{(y,l)}\} are fully determined by the iterative relation in Eq. (2).

Note that in the previous illustration, we focus on the short-time regime so that the lowest-order term in the error is the leading term. Nevertheless, this is not always the case, especially since long-time simulation is often of greater interest. To guarantee the accuracy of the simulation for long evolution times, researchers introduced a partition of the total evolution time. For an arbitrarily long time t, we can uniformly divide it into r steps, each t/r long. Therefore, the simulation error in each step is \mathcal{O}\left((\Vert H\Vert t/r{)}^{p+1}\right). The overall error, consisting of r single-step errors, can be bounded by the triangle inequality as

\begin{array}{rl}\Vert {\mathcal{U}}_p(t/r{)}^r-e^{-\mathrm{i}Ht}\Vert & \le \sum _{i=0}^{r-1}\Vert {\mathcal{U}}_p(t/r)-e^{-\mathrm{i}Ht/r}\Vert \\ & =\mathcal{O}\left(\frac{(\Vert H\Vert t{)}^{p+1}}{r^p}\right).\end{array}

Since the order of the formula is always p\ge 1, we can suppress the error arbitrarily by increasing r.

So far, we have introduced the general construction of an arbitrary high-order formula for simulation, along with its deviation from the ideal evolution. Nevertheless, this analysis of error is far from tight, and researchers eventually came up with better error analyses. The key observation is that we have used the norm \Vert H\Vert to bound all kinds of commutators in the previous error analysis. In [24], Somma returned to a commutator-based analysis and explored the possible improvement arising from this expression. Subsequent work by Childs et al. also proved a more compact and closed-form expression of the simulation error:

Theorem 3.1 (Trotter error with commutator scaling [3]). Let H=\sum _{l=1}^L{H}_l be a Hermitian operator consisting of L summands, and let 0\le \Vert H\Vert t\le 1. Let {\mathcal{U}}_p(t) be the pth-order Trotter-Suzuki formula in Eq. (2). Define {\tilde{\alpha }}_{\mathrm{c}\mathrm{o}\mathrm{m}\mathrm{m},p}:=\sum _{l_1,\dots ,l_{p+1}=1}^L\Vert [H_{l_1},[H_{l_2},\dots ,[H_{l_p},H_{l_{p+1}}]]]\Vert. Then, the additive Trotter error can be asymptotically bounded as

\Vert e^{-\mathrm{i}Ht}-{\mathcal{U}}_p(t)\Vert =\mathcal{O}({\tilde{\alpha }}_{\mathrm{c}\mathrm{o}\mathrm{m}\mathrm{m},p}{t}^{p+1}).

Theorem 3.1 relates the Trotter-Suzuki errors with the nested-commutator term of the target Hamiltonian. In fact, this nested-commutator representation offers advantages for estimating simulation errors in typical Hamiltonian models. As one of the most commonly studied models, the nearest-neighbor Hamiltonian serves as a good demonstration:

H=\sum _{(i,j)\in \mathrm{\Lambda }}{H}_{i,j},\text{with}\Vert H_{i,j}\Vert \le 1,

where \mathrm{\Lambda } denotes a D-dimensional lattice, and (i,j)\in \mathrm{\Lambda } means sites i and j belong to some nearest-neighbor edge. Suppose the system has n sites, and we can estimate the norm \Vert H\Vert =\mathcal{O}\left(n\right), leading to the simulation error \mathcal{O}\left((nt{)}^{p+1}\right) according to the ordinary analysis. In contrast, we get a better error estimation using the nested-commutator bound in Theorem 3.1. To view this, let us peel out the nested commutator layer by layer. For the innermost layer, every term H_{l_{p+1}} anti-commutes (or overlaps) with at most a constant number of choices for H_{l_p}. The resulting non-zero [H_{l_p},H_{l_{p+1}}]’s are still geometrically local, implying that there are only a constant number of choices for H_{l_{p-1}} that ensure the second level of the commutator to be non-zero. Therefore, we can continue this derivation and iterate over all layers. Given that the order p is a constant, the summation of nested commutators effectively sums over the first layer, thus implying {\tilde{\alpha }}_{\mathrm{c}\mathrm{o}\mathrm{m}\mathrm{m},p}=\mathcal{O}\left(n\right). The overall error is therefore \mathcal{O}\left(n{t}^{p+1}\right), which is a substantial improvement over the previous estimate.

3.3 Advanced analyses for Trotter-suzuki formula

Despite the theoretical advances for asymptotic error analyses introduced previously, empirical studies [26] have revealed that implementing these formulas often incurs substantial practical overheads, which can be prohibitive. This gap between theoretical efficiency and practical performance has motivated research to further improve the analysis of the simulation error. The idea stems from the observation that while the previous errors in the operator-norm distance represent worst-case estimates, we usually know the specific setting in simulation tasks. In this part, we introduce efforts to adopt this extra information to facilitate better and more accurate estimations of the simulation error.

As a natural complement to the worst-case analyses, we also want an average-case analysis of the simulation error. This setting was analyzed in [27,28], where the authors consider the average-case performance

(3)R={\mathbb{E}}_{\psi \sim \mathrm{H}\mathrm{a}\mathrm{a}\mathrm{r}}\Vert U^{\mathrm{\prime }}|\psi ⟩-e^{-\mathrm{i}Ht}|\psi ⟩{\Vert }_{l_2},

where the input states are randomly picked under the Haar measure. They proved the following average-case bound:

Theorem 3.2 (Average-case Trotter error, rephrased from [27]). Let H=\sum _{l=1}^L{H}_l be a Hermitian operator consisting of L summands, and let 0\le \Vert H\Vert t\le 1. Let {\mathcal{U}}_p(t) be the pth-order Trotter-Suzuki formula in Eq. (2). Define T_p:=\sum _{l_1,\dots ,l_{p+1}=1}^L\Vert [H_{l_1},[H_{l_2},\dots ,[H_{l_p},H_{l_{p+1}}]]]\Vert {}_F. Then, the average error in the {\ell }_2 norm for a 1-design input ensemble {\mu }_1 has the asymptotic upper bound

{\mathbb{E}}_{\psi \sim {\mu }_1}\Vert {\mathcal{U}}_p(t)|\psi ⟩-e^{-\mathrm{i}Ht}|\psi ⟩{\Vert }_{l_2}=\mathcal{O}\left(T_p{t}^{p+1}\right).

This average-case analysis switches from the nested commutators’ operator norms to normalized Frobenius norms, which always reduces the error. In the typical case of a nearest-neighbor Hamiltonian, we find that the new coefficient T_p scales sublinearly with the system size, namely T_p=\mathcal{O}\left(\sqrt{n}\right). The overall error under this setting is thus \mathcal{O}\left(\sqrt{n}{t}^{p+1}\right) versus the worst-case analysis \mathcal{O}\left(n{t}^{p+1}\right), implying a quadratic speedup.

The operator norm \Vert {\mathcal{U}}_p(t)-e^{-\mathrm{i}Ht}\Vert captures the worst-case deviation over all input states, which is the right notion for black-box simulation guarantees but can be overly pessimistic in settings where (i) the input is typical (e.g., drawn from a unitary design) or (ii) one ultimately cares about average performance or observable estimation. The quantity in Eq. (3) measures a state-dependent error \Vert ({\mathcal{U}}_p(t)-e^{-\mathrm{i}Ht})|\psi ⟩{\Vert }_{{\ell }_2} averaged over an ensemble. This choice also explains why the resulting bounds involve Frobenius-norm quantities rather than operator-norm quantities.

Beyond average-case analyses, there are also efforts focusing on practically relevant classes of input states. In [29], the authors consider states in the low-energy subspace of the Hamiltonian to obtain better simulation performance. They show that when the initial state has support only on energies up to some cutoff \mathrm{\Delta }, the error of standard Trotter-Suzuki product-formula simulation can be bounded using an effective low-energy norm {\mathrm{\Delta }}^{\prime } rather than the full operator norm \Vert H\Vert, and they prove that leakage out of the low-energy subspace during each Trotter step is exponentially suppressed in the energy gap to higher states.

Theorem 3.3 (Low-energy subspace error bound, rephrased from [29]). Let H=\sum _{\ell =1}^L{H}_{\ell } be a k-local Hamiltonian as above, H_{\ell }\ge 0, \mathrm{\Delta }\ge 0, 0\le \Vert H\Vert t\le 1, and {\mathcal{U}}_p(t) a p-th order product formula as in Eq. (2). Then,

\Vert (e^{-\mathrm{i}Ht}-{\mathcal{U}}_p(t)){\mathrm{\Pi }}_{\le \mathrm{\Delta }}\Vert =\mathcal{O}\left(({\mathrm{\Delta }}^{\prime }t{)}^{p+1}\right),

where

{\mathrm{\Delta }}^{\prime }=\mathrm{\Delta }+{\beta }_1\Vert H\Vert \log \left(\frac{{\beta }_2}{\Vert H\Vert t}\right)+{\beta }_3{\Vert H\Vert }^{2}nt,

and the {\beta }_i’s are positive constants, {\beta }_2\ge 1.

As a result, for local spin Hamiltonians with bounded degree, the number of Trotter steps needed to simulate time evolution scales better (often removing an explicit system-size factor) than the best known worst-case bounds, especially for long time evolutions or high-precision regimes. Overall, this establishes a general framework for Hamiltonian simulation “at low energy”, suggesting that physically relevant states—like ground states and low-lying excitations in many-body systems—can often be simulated on a quantum computer with asymptotically lower cost than generic states.

In [30], Zhao et al. discovered that entangled states can help accelerate quantum simulation, thereby connecting, for the first time, one of the most fundamental quantum tasks to one of the most remarkable quantum phenomena (see Theorem 3.4). The idea comes from a specific class known as k-uniform states, whose k-partite reduced density matrices are always maximally mixed. This type of entangled state masks its information from local measurements, which can be tailored to make local errors ineffective. Combining this observation with the nested-commutator errors in the Trotter-Suzuki formula, we can systematically eliminate the ineffective error terms.

Theorem 3.4 (Entanglement-based error bound [30]). For a given pure quantum state |\psi ⟩ and the pth-order Trotter-Suzuki formula {\mathcal{U}}_p, the error in a Trotter step of duration t has the upper bound

\begin{array}{rl}& \Vert (e^{-\mathrm{i}Ht}-{\mathcal{U}}_p)|\psi ⟩\Vert \\ & =\mathcal{O}\left(t^{p+1}\sqrt{\sum _{j,j^{\prime }}\Vert E_j^{†}{E}_{j^{\prime }}\Vert \sqrt{\log (d_{{\rho }_{j,j^{\prime }}})-S({\rho }_{j,j^{\prime }})}}+t^{p+1}\Vert E{\Vert }_F\right),\end{array}

where S(\cdot ) denotes the von-Neumann entropy, E=\sum _j{E}_j is the leading error term of {\mathcal{U}}_p, and {\rho }_{j,j^{\prime }}:={\mathrm{t}\mathrm{r}}_{[N]\setminus \mathrm{s}\mathrm{u}\mathrm{p}\mathrm{p}(E_j{E}_j^{\prime })}(|\psi ⟩⟨\psi |) is the reduced density matrix of |\psi ⟩⟨\psi | on the subsystem of \mathrm{s}\mathrm{u}\mathrm{p}\mathrm{p}(E_j{E}_j^{\prime }).

With this entropy-related error bound, the simulation error of a highly entangled state under the typical nearest-neighbor Hamiltonian model scales as \mathcal{O}\left(\sqrt{n}{t}^{p+1}\right). It is worth noting that classical simulation of quantum dynamics often requires low entanglement to be tractable, but this result suggests that entanglement accelerates quantum simulation. This novel finding can widen the gap between classical and quantum computational power, which is essential for exploring and demonstrating quantum advantage.

In addition to specific input states, certain families of observables provide significant advantages for advancing the analysis of product-formula simulations. A notable example is the work by Heyl et al. [31], who demonstrated error reductions associated with local observables in the long-time limit. They report a transition in the errors of time-averaged simulated states as the length of Trotter steps increases, indicating that the error becomes smaller than initially expected when the step length is below a certain threshold. Through perturbation theory, the authors elucidate the source of this advantage in the long-time regime and connect their findings to the well-known kicked-top model [32], offering valuable insights.

As an inverse direction, Childs et al. [3] derived short-time advantages from local observables by leveraging the local interactions of the target Hamiltonian. They effectively constrain the propagation of the light cone in Heisenberg’s picture during the simulation, which reduces the contributions of certain evolution terms in the final measurement, thereby minimizing simulation error. This concept of light-cone analysis is also foundational to understanding error reductions associated with local observables in analog simulations, as illustrated in [33].

Recent work [34,35] provides a comprehensive study of the benefits of observable information, addressing both short-time and arbitrary-time scenarios. For the short-time case, [34] analyzes light cones for both local observables and global observables composed of local summands, demonstrating the optimality of light-cone analyses using tailored product formulas. The error can be reduced so that it is independent of system size. In the context of arbitrary-time analysis, simulation can get further advantages from observables under average-case performance with random inputs or sufficient entanglement [34,35]. Typically, these results show that Pauli-summation observables significantly reduce simulation errors compared to more generic cases. Furthermore, [35] demonstrates that the Trotter error in observable growth is fundamentally bounded by cumulative operator scrambling, which provides a tighter upper bound on the Trotter error for general observable dynamics.

Theorem 3.5 (Accumulated scrambling-based bound, rephrased from [35]). Given ideal unitary U_0^r=e^{-\mathrm{i}Ht}, where t=r\delta t, and approximate Trotterized unitary {\mathcal{U}}_p^r with {\mathcal{U}}_p=U_0(I+\mathcal{M}), for the Hamiltonian H=\sum _{l=1}^L{H}_l and an input pure state |\psi ⟩, the additive Trotter error of an observable O can be bounded as

{ϵ}_O\le \sum _{k=1}^r\sqrt{⟨{\psi }_k||[O(k\delta t),\mathcal{M}]{|}^2|{\psi }_k⟩}\delta t^{p+1},

where |{\psi }_k⟩={\mathcal{U}}_p^{r-k}|\psi ⟩ denotes the kth step evolved state and O(k\delta t)=e^{\mathrm{i}Hk\delta t}O{e}^{-\mathrm{i}Hk\delta t}.

4 Variants of product formula

In this section, we review and introduce refinements of the standard Trotter–Suzuki product-formula approach to Hamiltonian simulation. Product-formula methods are among the most accessible and experimentally friendly techniques: they perform extremely well for geometrically local Hamiltonians and implement time evolution through sequences of simple, structured gates. Yet, in regimes with many non-local terms, naive deterministic Trotterization can require a large number of steps to control the error. This tension has driven lines of work that do not discard Trotterization, but instead systematically improve it to suppress coherent error and tighten worst-case bounds.

4.1 Advanced product formulas

In [17], Childs et al. show that by randomly permuting the summands in a Trotter-Suzuki decomposition, one can significantly tighten error bounds compared to a purely deterministic ordering. This enhancement is based on the so-called mixing lemma, which was first proposed in [15,16].

Lemma 4.1 (Mixing lemma [15]). Let V be a target unitary, with associated channel \mathcal{V}(\rho )=V\rho V^{†}. Let a,b>0 and \{U_1,U_2,\dots ,U_n\} be a set of unitaries such that

1. for all j\in \{1,\dots ,n\} we have \Vert U_j-V\Vert \le a;

2. there exist positive numbers \{p_j\} such that \sum _{j=1}^n{p}_j=1 and \Vert (\sum _j{p}_j{U}_j)-V\Vert \le b.

It follows that \mathcal{E}=\sum _j{p}_j{\mathcal{U}}_j satisfies

\Vert \mathcal{E}-\mathcal{V}{\Vert }_{\diamond }\le a^2+2b.

The mixing lemma states that an ideal unitary can be approximated with reduced error by randomly mixing a set of unitary channels rather than employing a single fixed ordering. Randomizing the ordering before each Trotter step mitigates systematic, ordering-dependent errors, leading to improved accuracy for a given circuit depth. This yields a quadratic reduction in gate complexity relative to standard product-formula constructions.

Underpinned by this lemma, Campbell [14] introduced the quantum stochastic drift (qDRIFT) protocol, which also incorporates randomness into Trotter-like product formulas but in a distinct manner. Instead of reordering the Hamiltonian terms, qDRIFT samples a single term at each Trotter step with probability proportional to its interaction strength, then “drifts” through the Hamiltonian in a stochastic fashion. This importance-sampling perspective means that each short evolution is drawn from a probability distribution dictated by the Hamiltonian coefficients.

Theorem 4.1 (Error bound of qDRIFT, rephrased from [14]). For a Hamiltonian H=\sum _{j=1}^L{h}_j{H}_j with \Vert H_j\Vert \le 1, h_j>0, and \lambda =\sum _{j=1}^L{h}_j, the qDRIFT protocol using r segments approximates the unitary channel \mathcal{U}(\rho )=e^{-\mathrm{i}Ht}\rho e^{\mathrm{i}Ht} with diamond-norm error bounded by

{\Vert {\mathcal{E}}_{\mathit{\text{qDRIFT}}}-\mathcal{U}\Vert }_{\diamond }\le \frac{2{\lambda }^2{t}^2}{r}.

Notably, qDRIFT achieves an L-independent complexity, making it a significant improvement over standard Suzuki-Trotter methods. The cost of this method is a large number of experimental repetitions to ensure that the approximate channel is close to the ideal one.

Randomized product-formula constructions do not define a single implemented unitary U^{\prime }. Instead, they define a distribution over unitaries and therefore implement a CPTP map

\mathcal{E}(\rho )={\mathbb{E}}_{\omega }\left[U_{\omega }\rho U_{\omega }^{†}\right].

In this setting, the operator norm \Vert U-U^{\prime }\Vert is not the appropriate metric, because there is no single U^{\prime } and the output is generally mixed; the standard worst-case error measure is the diamond norm \Vert \mathcal{E}-\mathcal{U}{\Vert }_{\diamond }, which quantifies distinguishability of channels.

Further, Nakaji et al. [36] introduced qSWIFT, a high-order randomized algorithm for Hamiltonian simulation. The number of gates needed for a specific precision in qSWIFT is independent of the Hamiltonian’s term count, and the systematic error decreases exponentially with the order parameter. qSWIFT is a higher-order randomized scheme; unlike many Trotter methods, it uses only a single ancilla qubit.

Lastly, Granet and Dreyer [37] introduced a Hamiltonian simulation algorithm that eliminates Trotter discretization error by reformulating small-angle gates in terms of randomly applied large-angle gates. Instead of requiring finely discretized time evolution steps, their approach constructs an unbiased estimator where small-angle rotations are probabilistically replaced with larger rotations, ensuring that the correct evolution is reconstructed in expectation. This approach achieves an ϵ-independent average gate count.

4.2 Multi-product formulas

One of the well-known variants of PF for quantum simulation is multi-product formulas (MPFs), which were initially introduced by Childs and Wiebe [7]. MPFs leverage the LCU technique to coherently superpose multiple low-order product formulas, enabling the approximation of higher-order product formulas. This approach offers a viable and practical framework for near-term quantum simulations. MPFs generalize the construction procedure of {\mathcal{U}}_p(t) by incorporating sums of product formulas. This approach simplifies the approximation process, as constructing a Taylor series through polynomial addition is more straightforward than relying solely on multiplication. This method requires only \mathcal{O}(p^2) matrix exponentials to achieve an approximation of U(t) with an error of \mathcal{O}(t^{2p+1}). Concretely, Childs and Wiebe [7] proved the following lemma for the MPFs:

Lemma 4.2 (Error bound for MPFs, rephrased from [7]). Let U(t)=e^{-\mathrm{i}Ht} be the exact time-evolution operator generated by a Hamiltonian H, and let {\mathcal{U}}_p(\lambda ) denote the pth-order Trotter–Suzuki product formula for approximating e^{-\mathrm{i}H\lambda }. Define

M_{k,p}(t)=\sum _{j=1}^{k+1}{C}_j{\mathcal{U}}_p(t/l_j{)}^{l_j},

where C_j and l_j are suitable real parameters. For sufficiently small t\ll 1, the approximation error satisfies

\Vert M_{k,p}(t)-U(t)\Vert \in \mathcal{O}\left(t^{2(k+p)+1}\right).

The MPFs M_{k,p}(t) have a smaller approximation error when approximating the target unitary U(t) than {\mathcal{U}}_p(\lambda ). The implementation of MPFs requires the LCU technique to coherently add k+1 terms of {\mathcal{U}}_p(t/l_j{)}^{l_j}. Since k is a constant factor, the number of ancillary qubits required to implement MPFs is \mathcal{O}(\log k), which may be practical for quantum simulations.

Nowadays, much attention has been paid to MPFs [38–40,97]. In particular, Carrera Vazquez et al. [39] and Rendon et al. [97] proposed a simpler implementation of MPFs for observable estimation without coherent control, i.e., without the LCU technique. Instead, one may implement each circuit {\mathcal{U}}_p(t/l_j{)}^{l_j} independently on a quantum device and then classically post-process the measured data. The idea of this variant of MPFs is to approximate the ideal channel using a linear combination of density matrices, i.e.,

\mu (t)=\sum _j{c}_j{\rho }_{l_j}(t)=\sum _j{c}_j{\mathcal{U}}_p(t/l_j{)}^{l_j}{\rho }_0{\mathcal{U}}_p^{†}(t/l_j{)}^{l_j},

where {\rho }_0 is the initial state of the quantum system. This suggests that the classical-mixture version of MPFs can be used to reduce the error in predicting the properties of many-body systems by evaluating \mathrm{T}\mathrm{r}(\mu (t)O).

More recently, based on [39,97], Zhuk, Robertson, and Bravyi demonstrated that MPFs could potentially achieve a quadratic reduction in Trotter error in 1-norm on arbitrary time intervals, without necessitating an increase in the circuit depth or qubit connectivity [40]. Furthermore, they proposed dynamic multi-product formulas, using time-dependent coefficients optimized to minimize a certain efficiently computable proxy for the Trotter error.

4.3 Compensation-based refinements

In this subsection, we review and introduce a complementary line of work that augments standard Trotter–Suzuki decompositions with an additional compensation stage. The common goal is to retain the simplicity and hardware accessibility of product-formula simulation while systematically cancelling its leading errors and thereby improving both asymptotic accuracy and gate complexity. The basic idea is to factor the ideal evolution over a timestep t as

e^{-\mathrm{i}Ht}\approx {\mathcal{U}}_p(t)C_p(t),

where {\mathcal{U}}_p(t) is the pth-order Trotter-Suzuki formula and C_p(t) is a tailored correction unitary. One derives C_p(t) by expanding both e^{-\mathrm{i}Ht} and {\mathcal{U}}_p(t) in Taylor (or Magnus-style) series and identifying the first nonvanishing error terms. After determining the explicit form of C_p(t), a naive Monte Carlo implementation using randomly sampled Pauli rotation gates yields the correction. Nevertheless, this can have a large normalization overhead (and hence large sampling cost).

To address this concern, Zeng et al. [41] introduce the Pauli-pairing technique. They deliberately group the higher-order anti-Hermitian Pauli operators with a lower-order identity operator to generate a new unitary, which can substantially reduce the sampling overhead. In effect, the simulator still looks like a product formula at hardware level, but its per-step error behaves like that of a much higher-order scheme.

Zeng et al. also propose Trotter–LCU compensation constructions that make the compensation part coherent. After each pth-order Trotter step {\mathcal{U}}_p(t), one appends a small ancilla-assisted block that implements a carefully organized linear combination of unitaries approximating the remainder C_p(t), using few ancilla qubits and a handful of controlled Pauli rotations. Because the lowest-order error terms of {\mathcal{U}}_p(t) are known analytically (they are fixed nested commutators of the local Hamiltonian terms), the correction can be sculpted so that all contributions up through order p cancel and even the next orders are paired into anti-Hermitian rotations with suppressed weight. As a result, the composite step achieves better time and accuracy scalings than an uncompensated 2Kth-order Suzuki formula: the effective time-step exponent improves from 1+1/K to 1+1/(2K+1), and the dependence on the target precision \epsilon is nearly polylogarithmic rather than polynomial. Crucially, this preserves the geometric locality and nested-commutator cost of Trotterization (favorable for lattice Hamiltonians), while importing the high-precision behavior usually associated with truncated-Taylor algorithms which we will introduce shortly.

4.4 Algorithmic error mitigation

In this subsection, we review error-mitigation-based Hamiltonian simulation methods [42–44] that achieve the optimal logarithmic dependence on 1/ϵ. Observing that Trotter error admits a smooth expansion in the inverse step size s=1/r for time-evolved observables, one can iteratively perform Richardson extrapolation to improve the quality of the estimate. By evaluating the same observable at multiple step sizes and forming calibrated linear combinations, one can cancel the leading orders at each iteration and effectively extrapolate to the s\to 0 limit. Intuitively, product formulas incur deterministic commutator errors whose sign and magnitude vary systematically with step size; combining runs at different step sizes cancels these structured biases, leading to an optimal logarithmic dependence of the Trotter error.

Endo et al. [42] introduce algorithmic error mitigation, extending hardware error mitigation techniques to cancel Trotter error. They characterize the fundamental trade-off between algorithmic error and accumulated physical noise, identifying an optimal Trotter step N_{\mathrm{o}\mathrm{p}\mathrm{t}}=\sqrt{\alpha /\beta } where \alpha ,\beta quantify algorithmic and physical error scales. By measuring observables at m steps and applying a Richardson-type linear combination, the method cancels the first m terms of the error expansion, yielding a residual error O({ϵ}^{(m+1)}). However, the variance increases by an exponentially large factor, resulting in a sampling cost that can be significant for high-order cancellation.

Watson and Watkins [43] establish rigorous complexity guarantees for Richardson and Chebyshev–grid extrapolation applied to pth–order Trotter–Suzuki formulas. They show that for time t and precision ϵ, one can estimate observables with circuit depth scaling as

O\left(t^{1+1/p}\mathrm{p}\mathrm{o}\mathrm{l}\mathrm{y}\mathrm{l}\mathrm{o}\mathrm{g}(1/ϵ)\right),

which yields an exponential improvement in ϵ–dependence relative to fixed–step Trotterization. Their analysis quantifies the conditioning of interpolation procedures and compares incoherent (classical sampling) versus coherent (amplitude–estimation) variants.

Watson [44] further extends extrapolation techniques to the randomly compiled setting, proposing qFLO for expectation-value estimation under random Hamiltonian sampling. The protocol applies Richardson extrapolation to qDRIFT, whose depth depends only on the simulation time and Hamiltonian norm, not on the number of Hamiltonian terms. This approach inherits the hardware friendliness of randomized compilations but does not exploit higher-order commutator structure; asymptotically, higher-order Trotter formulas can still provide better scaling in t.

Theorem 4.2 (qFLO resource scaling [44]). Suppose we wish to obtain the expectation value of an observable A after evolving an initial state {\rho }_0 for time t under a Hamiltonian H=\sum _{i=1}^L{h}_j{H}_j with \Vert H_j\Vert =1. By varying the time–step size of the qDRIFT mapping and taking measurements at m different values of the time–step, with high probability it is possible to construct an order–m Richardson estimator \hat{A} satisfying

|\hat{A}-\mathrm{t}\mathrm{r}\left[{\rho }_0{e}^{\mathrm{i}Ht}A{e}^{-\mathrm{i}Ht}\right]|\le ϵ\Vert A\Vert ,

using m=\tilde{O}(\log (1/ϵ)) sampling points and circuit depths of no more than

O\left((\lambda t{)}^2\log \left(\frac{1}{ϵ}\right)\right),

where \lambda =\sum _j{h}_j. Each sample point used to construct the estimator requires O(1/{ϵ}^2) runs of the qDRIFT mapping. As a result the total gate count scales as

O\left(\frac{1}{{ϵ}^2}(\lambda t{)}^2{\log }^2\left(\frac{1}{ϵ}\right)\right).

5 Truncated Taylor series method

In addition to PF, another promising class of quantum simulation algorithms is post-Trotter methods. Given an exponential function of an operator H, post-Trotter protocols can simulate e^{-\mathrm{i}Ht}=\sum _k{g}_k{f}_k(H), where f_k(H) is the kth function of H and these functions can be power or trigonometric functions. For example, Taylor expansion gives e^{-\mathrm{i}Ht}=\sum _{k=0}^{\mathrm{\infty }}\frac{(-\mathrm{i}t{)}^k}{k!}{H}^k, and the Jacobi-Anger expansion gives e^{-\mathrm{i}Ht}=\sum _{k=-\mathrm{\infty }}^{\mathrm{\infty }}{\mathrm{i}}^k{J}_k(t)e^{\mathrm{i}k{\cos }^{-1}H}, where J_k(t) is a Bessel function of the first kind. These two expansions correspond to the TS method [5] and QSP [8,45] for quantum simulation, respectively.

In this section, we review the TS method [5]. Post-Trotter approaches are characterized by the fact that the expansion of the objective function—such as the exponential function—decomposes into a linear combination of functions of H. Consequently, the exponential function e^{-\mathrm{i}Ht} can be implemented by multiple accesses to H. In general, Hamiltonian simulation using post-Trotter protocols requires a block-encoding of the matrix H in a quantum circuit or quantum state. Typically, the LCU technique [7] is the standard way to implement a block-encoding of H [8]. LCU is probabilistic and therefore requires oblivious amplitude amplification [46].

5.1 Linear combination of unitaries

Here we briefly review the LCU method [5,6]. Consider the task of implementing an arbitrary n-qubit unitary operator V on a quantum computer. It is known that V can be expressed as a linear combination of unitary operators: V=\sum _{j=1}^J{\alpha }_j{U}_j, where each U_j is a unitary, and {\alpha }_j\in {\mathbb{R}}^{+} are coefficients. The core idea of the LCU technique involves introducing an ancillary register to encode the coefficients {\alpha }_j. Specifically, assume the existence of a unitary G acting on the ancillary space such that

G|\mathbf{0}⟩=\frac{1}{\sqrt{\alpha }}\sum _{j=1}^J\sqrt{{\alpha }_j}|\mathbf{j}⟩,

where |\mathbf{0}⟩ is the initial ancillary state, |\mathbf{j}⟩ forms an orthonormal basis of the ancillary system, and \alpha =\sum _{j=1}^J{\alpha }_j is a normalization constant. Next, define the select operation \mathrm{select}(U):=\sum _{j=1}^J|\mathbf{j}⟩⟨\mathbf{j}|\otimes U_j, which acts conditionally on the ancillary register. Using these components, we construct the operator

W:=(G^{†}\otimes \mathbb{I})\mathrm{select}(U)(G\otimes \mathbb{I}),

acting on the combined ancillary and system Hilbert space. Applying W to the initial state |\mathbf{0}⟩\otimes |\psi ⟩, we obtain

W|\mathbf{0}⟩\otimes |\psi ⟩=\frac{1}{\alpha }|\mathbf{0}⟩\otimes V|\psi ⟩+\frac{\sqrt{{\alpha }^2-1}}{\alpha }|\varphi ⟩,

where |\varphi ⟩ is a state orthogonal to the subspace spanned by |\mathbf{0}⟩\otimes {\mathcal{H}}_{\text{sys}}. Projection onto the ancillary state |\mathbf{0}⟩ thus probabilistically yields the desired transformation V on the target system, with success probability 1/{\alpha }^2. In scenarios where the success probability 1/{\alpha }^2 is insufficient, amplitude amplification techniques can be employed to boost the probability. When \alpha =2, the oblivious amplitude amplification protocol [83] enables deterministic implementation via the iteration

(-WR{W}^{†}RW)|\mathbf{0}⟩\otimes |\psi ⟩=|\mathbf{0}⟩\otimes V|\psi ⟩,

where R:=(\mathbb{I}-2|\mathbf{0}⟩⟨\mathbf{0}|)\otimes \mathbb{I} acts as a reflection operator on the combined ancillary and system Hilbert space. For cases where \alpha >2, the unitary V can be decomposed into a product of unitaries V^{(j)}, each corresponding to a sublinear combination with {\alpha }^{(j)}=2. Employing this decomposition, the overall operation can be implemented deterministically using on the order of \mathcal{O}(\alpha ) queries.

5.2 Hamiltonian simulation with truncated Taylor series

Considering the truncated Taylor series algorithm [5], one can expand the operator e^{-\mathrm{i}Ht} to the Kth order as follows: e^{-\mathrm{i}Ht}\approx \sum _{k=0}^K\frac{(-\mathrm{i}t\sum _l{\alpha }_l{H}_l{)}^k}{k!}+O(t^K), where H=\sum _{l=1}^L{\alpha }_l{H}_l and each H_l is a Pauli matrix. The first summation term is realized via the LCU method. The approximation error arises from truncating the Taylor series expansion.

Given a Hamiltonian H=\sum _{l=1}^L{\alpha }_l{H}_l, where all H_l are n-fold tensor products of Pauli matrices. According to TS, one has

e^{-\mathrm{i}Ht}\approx \sum _{k=0}^K\sum _{l_1,...l_k=1}^L\frac{(-\mathrm{i}t{)}^k}{k!}{\alpha }_{l_1}\dots {\alpha }_{l_k}{H}_{l_1}\dots H_{l_k}.

Without loss of generality, we can assume each {\alpha }_l>0. One can rewrite the expansion as V=\sum _j{\beta }_j{U}_j, where {\beta }_j={\beta }_{(k,l_1,\dots ,l_k)}:= \frac{t^k}{k!}{\alpha }_{l_1}\cdots {\alpha }_{l_k} and U_j=U_{(k,l_1,\dots ,l_k)}:=(-i{)}^k{H}_{l_1}\cdots H_{l_k}. That is, one can make use of the LCU technique to implement quantum simulation using the above expansion.

To ensure that the truncated Taylor series approximates the evolution operator with the desired accuracy ϵ, the order K for the truncated Taylor series should be [5]

K=\mathcal{O}\left(\frac{\log (\alpha t/ϵ)}{\log \log (\alpha t/ϵ)}\right),

where \alpha =\sum _l\left|{\alpha }_l\right|. In order to effectively achieve oblivious amplitude amplification in LCU, one needs to make the segment normalization satisfy {\alpha }_{\mathrm{T}\mathrm{S}}\approx 2. In general, however, {\alpha }_{\mathrm{T}\mathrm{S}} need not equal 2 a priori. Similar to the product formula method, here we should divide the simulation time t into r segments such that the simulation time in each segment \tau =t/r=\frac{\ln 2}{\sum _l{\alpha }_l} satisfies \sum _{k=0}^K\frac{(t\sum _l{\alpha }_l/r{)}^k}{k!}\approx e^{t\sum _l{\alpha }_l/r}=2.

For the remaining time , where \tau =\frac{\ln 2}{\sum _l{\alpha }_l} and t= (r-1)\tau +{\tau }_{re}, the whole evolution can be implemented as U(t)=V({\tau }_{re})V(\tau {)}^{r-1}. The algorithmic error mainly comes from the finite truncated Taylor series \Vert U(\tau )-e^{-\mathrm{i}H\tau }\Vert \le 2\frac{(\ln 2{)}^{K+1}}{(K+1)!}. It is easy to verify that

V(\tau )=\frac{3}{{\alpha }_{\mathrm{T}\mathrm{S}}}U(\tau )-\frac{4}{{\alpha }_{\mathrm{T}\mathrm{S}}^3}U(\tau )U(\tau )U^{†}(\tau ).

Now suppose the truncation error is \delta =2\frac{(\ln 2{)}^{K+1}}{(K+1)!}. Then the error for each segment can be bounded as

\Vert V(\tau )-e^{-\mathrm{i}H\tau }\Vert \le \frac{{\delta }^2+3\delta +4}{2}\delta .

Thus, the total error accumulates with r segments as \Vert V(t)-e^{-\mathrm{i}Ht}\Vert \le \frac{{\delta }^2+3\delta +4}{2}r\delta.

5.3 Unary encoding setting and gate complexity

For the analysis of the gate complexity of the Taylor Series (TS) method, the unary-encoding setting should be considered [5]. The idea behind unary encoding for the TS method is intuitive. It involves preparing an auxiliary system and applying control operations that allow the linear combination of different H^k with coefficients \frac{(\mathrm{i}t{)}^k}{k!}.

Specifically, the overall ancillary system (with K+1 registers) prepares B_{\text{unary}}|\mathbf{0}⟩=G_{\mathrm{T}\mathrm{S}}\otimes \underset{i=1}{\overset{K}{⨂}}{G}_{(i)}|\mathbf{0}⟩, where G_{\mathrm{T}\mathrm{S}}=\frac{1}{\sqrt{{\alpha }_{\mathrm{T}\mathrm{S}}}}\sum ^{K}k= 0\sqrt{\frac{(\alpha t/r{)}^k}{k!}}|1^k{0}^{K-k}⟩ and G_{(i)}=\sum _{l=1}^L\frac{1}{\sqrt{\alpha }}\sqrt{{\alpha }_l}{|l⟩}_i. Here, {\alpha }_{\mathrm{T}\mathrm{S}}=\sum _{k=0}^K\frac{(\alpha t/r{)}^k}{k!} and \alpha =\sum _{l=1}^L\Vert H_l\Vert.

Since the first register makes use of unary encoding, it consists of K=\mathcal{O}\left(\frac{\log (\alpha t/ϵ)}{\log \log (\alpha t/ϵ)}\right) qubits. The remaining K registers are used to encode the coefficient information of the Hamiltonian H=\sum _{l=1}^L{\alpha }_l{H}_l, with each register containing \log (L) qubits. Therefore, the total number of qubits required for the ancillary state is K\log (L)+K, i.e., \mathcal{O}\left(\frac{\log (L)\log (\alpha t/ϵ)}{\log \log (\alpha t/ϵ)}\right).

The elementary gate complexity to implement the unitary B_{\text{unary}} is \mathcal{O}(KL). This is because G_{\text{TS}} and G_i can be implemented using \mathcal{O}(K) and \mathcal{O}(L) single-qubit and CNOT gates, respectively [47]. Since B_{\text{unary}} is the tensor product of K+1 unitaries, the overall gate complexity for B_{\text{unary}} is \mathcal{O}\left(L\frac{\log (\alpha t/ϵ)}{\log \log (\alpha t/ϵ)}\right).

Now, let’s analyze the circuit construction and gate complexity of the \mathrm{s}\mathrm{e}\mathrm{l}\mathrm{e}\mathrm{c}\mathrm{t}(U_{\text{unary}}). The \mathrm{s}\mathrm{e}\mathrm{l}\mathrm{e}\mathrm{c}\mathrm{t}(U_{\text{unary}}) is given as

\mathrm{s}\mathrm{e}\mathrm{l}\mathrm{e}\mathrm{c}\mathrm{t}(U_{\text{unary}})=\underset{i=1}{\overset{K}{⨂}}(|0⟩{⟨0|}_i\otimes \mathbb{I}-\mathrm{i}|1⟩{⟨1|}_i\otimes \mathrm{s}\mathrm{e}\mathrm{l}\mathrm{e}\mathrm{c}\mathrm{t}(U_i)),

where \mathrm{s}\mathrm{e}\mathrm{l}\mathrm{e}\mathrm{c}\mathrm{t}(U_i)=\sum _{j=1}^L|\mathbf{j}⟩{⟨\mathbf{j}|}^i\otimes H_j^{(i)}. Clearly, this unitary is the tensor product of K controlled-select(H) operations, and it maps |k⟩|l⟩|\psi ⟩=|k⟩|l⟩(-i{H}_l{)}^k|\psi ⟩. The gate complexity for a single controlled-select(H) operation is \mathcal{O}(L(n+\log L)) [5], resulting in the overall gate complexity for all controlled-select(H) operations \mathcal{O}(L(n+\log L)K).

Now one can define

W_{\text{unary}}:=B_{\text{unary}}^{†}\otimes {\mathbb{I}}_s(\mathrm{s}\mathrm{e}\mathrm{l}\mathrm{e}\mathrm{c}\mathrm{t}(U_{\text{unary}}))B_{\text{unary}}\otimes {\mathbb{I}}_s,

such that

W_{\text{unary}}|\mathbf{0}⟩{|\psi ⟩}_s=\frac{1}{{\alpha }_{\mathrm{T}\mathrm{S}}}|\mathbf{0}⟩\otimes \left(\sum _{k=0}^K\frac{(-\mathrm{i}{)}^k(t/r{)}^k}{k!}{H}^k|\psi ⟩\right)+|\varphi ⟩,

where H=\sum _l^L{H}_l. For oblivious amplitude amplification [46], one chooses r such that

{\alpha }_{\mathrm{T}\mathrm{S}}=\sum _{k=0}^K\frac{(\alpha t/r{)}^k}{k!}\approx e^{\alpha t/r}\approx 2,

which is achieved by taking t/r=\ln 2/\alpha. Specifically, one can approximate |\mathbf{0}⟩\otimes e^{-\mathrm{i}Ht/r}|\psi ⟩ by

-W_{\text{unary}}R{W}_{\text{unary}}^{†}R{W}_{\text{unary}}|\mathbf{0}⟩\otimes |\psi ⟩=|\mathbf{0}⟩\otimes \sum _{k=0}^K\frac{(-\mathrm{i}{)}^k(t/r{)}^k}{k!}{H}^k|\psi ⟩.

By repeating the above process, the quantum dynamics U\approx (e^{-\mathrm{i}Ht/r}{)}^r=e^{-\mathrm{i}Ht} can be implemented. Combined with the gate counts for ancillary states, the total gate complexity for r segments (r=\mathcal{O}(\alpha t)) will be

\mathcal{O}\left(\alpha tL(n+\log L)\frac{\log (\alpha t)/ϵ}{\log \log (\alpha t/ϵ)}\right).

Overall, the gate complexity of truncated TS with LCU is summarized as follows:

Theorem 5.1 (Gate complexity of truncated Taylor series method, rephrased from [5]). For a Hamiltonian H=\sum _{l=1}^L{\alpha }_l{H}_l with {\alpha }_l>0 and \alpha =\sum _{l=1}^L{\alpha }_l, the truncated Taylor series method approximates e^{-\mathrm{i}Ht} with error ϵ using a truncation order of

K=\mathcal{O}\left(\frac{\log (\alpha t/ϵ)}{\log \log (\alpha t/ϵ)}\right),

and divides the evolution into r=\mathcal{O}(\alpha t) segments. The gate complexity is

\mathcal{O}\left(\alpha tL(n+\log L)\frac{\log (\alpha t/ϵ)}{\log \log (\alpha t/ϵ)}\right),

where each segment time satisfies \tau =t/r=\frac{\ln 2}{\alpha } to enable efficient oblivious amplitude amplification.

5.4 Advanced truncated Taylor series

Intuitively, Hamiltonians with more commutative terms are also easier to simulate on a quantum computer, and anti-commutative relations generally cause more errors, such as in the product formula method. In [48], Zhao and Yuan explored how anticommutation relations can be leveraged to bolster error bounds in truncated Taylor series quantum algorithms. Their key observation is that Hamiltonians with mutually anti-commuting terms can sometimes be simulated as efficiently as commutative ones, defying the usual intuition that commutation alone leads to simpler dynamics.

Theorem 5.2 (Exact simulation for anti-commuting Hamiltonians [48]). Consider a Hamiltonian H=\sum _{l=1}^L{\alpha }_l{H}_l acting on n qubits, where each H_l is a unitary operator and all terms are pairwise anti-commuting: \{H_{l_1},H_{l_2}\}=0 for l_1\ne l_2. The evolution operator U(t)=e^{-\mathrm{i}Ht} can be implemented exactly without algorithmic error using the LCU method with a gate complexity of \mathcal{O}(L^{\frac{3}{2}}(n+\log L)).

Their method reduces the required gate cost for Hamiltonians with strong anti-commutation properties, making it particularly effective for quantum chemistry applications. This property is leveraged to decrease algorithmic errors and gate complexity in the truncated Taylor series quantum algorithm for general problems.

6 Quantum signal processing and qubitization

Quantum Signal Processing (QSP) [8,45,49] and Quantum Singular Value Transformation (QSVT) [9] have become a central primitive for optimal-query-complexity Hamiltonian simulation algorithms. In particular, qubitization enables simulation with query complexity \mathcal{O}\left(t\alpha +\frac{\log (1/ϵ)}{\log (e+\log (1/ϵ)/(\alpha |t|))}\right) [8–9], which is both asymptotically optimal in time and achieves the known additive optimality in precision. The essential mechanism is the ability to implement bounded polynomial transformations of eigenvalues and singular values using a sequence of controlled phase rotations and reflections.

We begin by revisiting the original construction of Low, Yoder, and Chuang, where discrete single-qubit rotations realize polynomial responses of a two-level system. This construction arises from a discrete, analytically-synthesized analogue of optimal quantum control, producing a target polynomial transformation. We then show how this single-qubit formalism extends to arbitrary operators via block-encoding, leading to the qubitization framework.

Throughout, we highlight the algebraic constraints that characterize feasible QSP polynomials (parity, boundedness, and SU(2) structure) and their connections to Chebyshev expansions and Jacobi–Anger series. Finally, we review practical advances that enhance the implementability of QSP-based Hamiltonian simulation, including reliable phase-vector computation, generalized SU(2) signal operators, and improved \mathrm{S}\mathrm{E}\mathrm{L}\mathrm{E}\mathrm{C}\mathrm{T}/\mathrm{P}\mathrm{R}\mathrm{E}\mathrm{P}\mathrm{A}\mathrm{R}\mathrm{E} oracles. These developments lower classical compilation cost, reduce circuit depth, and improve oracle efficiency, thereby strengthening the practicality of QSP and qubitization for high-precision Hamiltonian simulation and other matrix-function algorithms.

6.1 Optimal quantum control for response function

Before we introduce a general framework of QSP (and its extensions like QSVT and qubitization), we review a result of Low, Yoder, and Chuang [49] for performing a transformation of a 2\times 2 matrix.

Suppose there is a single-qubit rotation about the x-axis

W(x):=\left(\begin{array}{cc}x& \mathrm{i}\sqrt{1-x^2}\\ \mathrm{i}\sqrt{1-x^2}& x\end{array}\right)=e^{\mathrm{i}{\sigma }_x\text{arccos(x)}}.

The goal of QSP is to generate a matrix whose first entry is a polynomial in x. To do this, one needs to introduce another rotation e^{{\varphi }_i{\sigma }_z} and intersperse it with W(x), giving

W_{\mathrm{\Phi }}(x):=e^{\mathrm{i}{\varphi }_0{\sigma }_{\mathrm{z}}}W(x)e^{\mathrm{i}{\varphi }_1{\sigma }_{\mathrm{z}}}W(x)\cdots W(x)e^{\mathrm{i}{\varphi }_{\mathrm{k}}{\sigma }_{\mathrm{z}}},

where \mathrm{\Phi }:=({\varphi }_0,{\varphi }_1,{\varphi }_2,...,{\varphi }_k). The following theorem shows that the desired polynomial response in x can be realized in this way:

Theorem 6.1 (Available functions [9]). There exists \mathrm{\Phi }\in {\mathbb{R}}^{k+1} such that

W_{\mathrm{\Phi }}(x)=\left(\begin{array}{cc}P(x)& \mathrm{i}Q(x)\sqrt{1-x^2}\\ \mathrm{i}Q(x)\sqrt{1-x^2}& P^{\ast }(x)\end{array}\right),

if and only if P,Q\in \mathbb{C}[x] satisfy (i) \mathit{\text{deg}}(P)\le k and \mathit{\text{deg}}(Q)\le k-1, (ii) P has parity k mod 2 and Q has parity k-1 mod 2, and (iii) \mathrm{\forall }x\in [-1,1],|P(x){|}^2+(1-x^2)|Q(x){|}^2=1.

Note that W(x) can be decomposed into a single-qubit reflection operator R(x) and single-qubit phase gates so that

W_{\mathrm{\Phi }}(x)={\mathrm{i}}^k{e}^{\mathrm{i}{\varphi }_0{\sigma }_z}\prod _{j=1}^k{e}^{-\mathrm{i}{\sigma }_z\pi /4}R(x)e^{\mathrm{i}{\sigma }_z({\varphi }_j+\pi /4)},

where

R(x):=\left(\begin{array}{cc}x& \sqrt{1-x^2}\\ \sqrt{1-x^2}& -x\end{array}\right).

We will see that this decomposition is an essential observation for the generalization of QSP to higher dimensions.

Low and Chuang initially proposed Quantum Signal Processing in [45]. The main idea was to implement quantum algorithms using quantum control, but with analytic rather than empirical synthesis. Starting from the theoretical framework in [49], which showed that a tuple of polynomials (A,B,C,D) can be generated by a chain of discrete single-qubit rotations {\hat{R}}_{\varphi }(\theta )=e^{-\mathrm{i}(\theta /2)({\hat{\sigma }}_x\cos \varphi +{\hat{\sigma }}_y\sin \varphi )}, the full chain gives

\begin{array}{rl}\hat{V}(\theta )& ={\hat{R}}_{{\varphi }_N}(\theta ){\hat{R}}_{{\varphi }_{N-1}}(\theta )\cdot \cdot \cdot {\hat{R}}_{{\varphi }_1}(\theta ),\overrightarrow{\varphi }\in {\mathbb{R}}^N,\\ & =A(\theta )\hat{\mathbb{I}}+\mathrm{i}B(\theta ){\hat{\sigma }}_z+\mathrm{i}xC(\theta ){\hat{\sigma }}_x+\mathrm{i}D(\theta ){\hat{\sigma }}_y,\end{array}

where the four arbitrary functions should satisfy certain constraints and depend on \overrightarrow{\varphi }\in {\mathbb{R}}^N. [49] established the constraints for \hat{V}(\theta ) to exist and suggested that it can be used to perform Hamiltonian simulation.

6.2 Quantum signal processing

Another useful lemma for QSP is as follows:

Lemma 6.1 (Lemma 14 in Ref. [8]). For any even integer Q>0, a choice of functions \mathcal{A}(\theta ) and \mathcal{C}(\theta ) is achievable by the framework of QSP if and only if the following are true:

1. \mathcal{A}(\theta )=\sum _{k=0}^K{a}_k\cos (k\theta ) is a real cosine Fourier series of degree at most K, where a_k are coefficients;

2. \mathcal{C}(\theta )=\sum _{k=1}^K{c}_k\sin (k\theta ) is a real sine Fourier series of degree at most K, where c_k are coefficients;

3. \mathcal{A}(0)=1+{ϵ}_1, where |{ϵ}_1|\le 1;

4. \mathrm{\forall }\theta \in \mathbb{R},{\mathcal{A}}^2(\theta )+{\mathcal{C}}^2(\theta )\le 1+{ϵ}_2, where {ϵ}_2\in [0,1].

Further, [45] shows that the Jacobi-Anger expansion of the exponential function satisfies these criteria and yields the stated Hamiltonian-simulation cost in the sparse-input model. Specifically,

\begin{array}{rl}{e}^{-\mathrm{i}\lambda t}& =J_0(t)+2\sum _{\mathrm{e}\mathrm{v}\mathrm{e}\mathrm{n}k>0}^{\mathrm{\infty }}(-1{)}^{\frac{k}{2}}{J}_k(t)T_k(\lambda )+2\mathrm{i}\sum _{\mathrm{o}\mathrm{d}\mathrm{d}k>0}^{\mathrm{\infty }}(-1{)}^{\frac{k-1}{2}}{J}_k(t)T_k(\lambda )\\ & =A(\lambda )+\mathrm{i}C(\lambda ),\end{array}

where J_k is the Bessel function of the first kind and T_k are Chebyshev’s polynomials.

Based on QSP, [45] proposed an optimal Hamiltonian simulation that achieves optimal complexity in all parameters. This approach has several drawbacks. First, the access model is the sparse oracle, which means that there would be an extremely high overhead if H cannot be expressed as a sparse matrix. Moreover, we do not have access to the black-box oracle, as we cannot guarantee that the non-zero elements are efficiently row-computable. Second, its black-box oracles can be challenging to realize. Avoiding the \mathcal{O}(2^n) blowup by exploiting sparsity requires that positions of non-zero elements are efficiently row-computable, which is not always the case. Lastly, the iterator requires a quantum walk [50], which doubles the number of qubits.

6.3 Qubitization

Following this work, Low and Chuang overcame these problems and proposed qubitization in [8]. In this work, they construct the iterator from a signal oracle and show that the iterator acts as a quantum walk in an invariant subspace. Therefore, the polynomial transformation of eigenvalues within such a subspace is not affected by the dynamics in other subspaces.

The major difference between QSP and qubitization is that QSP is purely a sum of SU(2) submatrices, whereas qubitization is a framework wherein a unitary matrix can be transformed into a standard form.

One of the key steps in qubitization is to implement a block-encoding of a signal matrix. A unitary U is a block-encoding of a matrix V if

(4)U=\left(\begin{array}{cc}V& \cdot \\ \cdot & \cdot \end{array}\right)=|\mathbf{0}⟩⟨\mathbf{0}|\otimes V+\cdots ,

where |0⟩ denotes the first computational basis state of the ancillary qubit. It is easy to verify that V=(⟨\mathbf{0}|\otimes \mathbb{I})U(|\mathbf{0}⟩\otimes \mathbb{I}). Obviously, LCU is one method for block-encoding matrices. Suppose V=\sum _i{\beta }_i{U}_i, where U_i is a unitary operator. According to Eq. (4), the signal operator V may be given as

⟨0{|}_{a}W|0{⟩}_a=⟨0{|}_a(G^{†}\otimes \mathbb{I})\text{select}(U)(G\otimes \mathbb{I})|0{⟩}_a=\frac{V}{\alpha },

where \alpha =\sum _i|{\alpha }_i|. Here we call W a signal oracle that encodes V using one query each to select(U), G, and G^{†}. Generally, V can be a linear combination of Hermitian operators, i.e., V=\sum _i{H}_i.

Finally, this construction achieves

Theorem 6.2 (Optimal Hamiltonian simulation by Qubitization [8]). Let (⟨G{|}_a\otimes {\hat{I}}_s)\hat{U}(|G{⟩}_a\otimes {\hat{I}}_s)=\hat{H}\in {\mathbb{C}}^{N\times N} be Hermitian for some unitary \hat{U}\in {\mathbb{C}}^{Nd\times Nd} and some state-preparation unitary \hat{G}|0{⟩}_a=|G{⟩}_a\in {\mathbb{C}}^d. Then e^{-\mathrm{i}\hat{H}t} can be simulated for time t, error ϵ in spectral norm, and failure probability \mathcal{O}(ϵ), using at most \lceil {\log }_2(N)\rceil +\lceil {\log }_2(d)\rceil +2 qubits in total, \mathrm{\Theta }(Q) queries to controlled-\hat{G}, controlled-\hat{U}, and their inverses, and \mathcal{O}(Q\log (d)) additional two-qubit quantum gates where

\begin{array}{rl}Q& =min\left\{q\in {\mathbb{Z}}^{+}:ϵ\ge \frac{4|t{|}^q}{q!2^q}=\mathcal{O}\left({\left(\frac{e\tau }{2q}\right)}^q\right)\right\}\\ & =\mathcal{O}\left(t+\log (1/ϵ)\right).\end{array}

The error tolerance ϵ in Theorem 6.2 and other block-encoding-based Hamiltonian simulations is measured in the partial operator norm. However, the overall procedure constitutes a quantum channel if one keeps the failure flag as part of the output, and it should therefore be assessed using the diamond norm. In qubitization and related block-encoding approaches, intermediate guarantees are often stated for the encoded block, i.e., scenarios where the measurement indicates failure are discarded. Consequently, the partial operator norm \Vert \cdot {\Vert }_{\mathrm{p}\mathrm{a}\mathrm{r}} is sufficient because the algorithm’s action on the system, conditioned on the ancilla indicating success, is precisely the intended transformed block. Once one conditions on success, the implemented operation on the system can be compared with e^{-iHt} under the standard operator norm, which corresponds exactly to the partial operator norm distance.

[51] exploits the structure in detail to achieve better performance in Hamiltonian simulation. Specifically, it achieves a polynomial speedup in some parameters for simulating Hamiltonian dynamics in the low-energy subspace.

Theorem 6.3 (Uniform spectral amplification of low-energy subspaces [51]). Given Hermitian standard-form (\hat{H},\alpha ,\hat{U},d) with eigenstates \hat{H}/\alpha |\lambda ⟩=\lambda |\lambda ⟩, let \mathrm{\Delta }\in (0,1) be a positive constant, and \hat{\mathrm{\Pi }}=\sum _{\lambda \in [-1,-1+\mathrm{\Delta }]}|\lambda ⟩⟨\lambda | be a projector onto the low-energy subspace of \hat{H}. Then there exists a standard-form ({\hat{H}}_{\mathrm{a}\mathrm{m}\mathrm{p}},\mathrm{\Delta }\alpha ,\hat{V},4d) such that

\Vert \hat{\mathrm{\Pi }}\left(\frac{{\hat{H}}_{\mathrm{a}\mathrm{m}\mathrm{p}}}{\mathrm{\Delta }\alpha }-\frac{\hat{H}+\alpha \hat{I}(1-\mathrm{\Delta })}{\mathrm{\Delta }\alpha }\right)\hat{\mathrm{\Pi }}\Vert \le ϵ,

and \hat{V} requires

\mathcal{O}\left({\mathrm{\Delta }}^{-1/2}{\log }^{3/2}\left(\frac{1}{\mathrm{\Delta }ϵ}\right)\right)

queries to controlled-\hat{U}.

This speedup stems from the construction of a polynomial that achieves spectral gap amplification by a factor {\mathrm{\Delta }}^{-1} using only \mathcal{O}\left({\delta }^{1/2}\right) queries. The result is stored in a standard form of block-encoding. Using a normal procedure to execute Hamiltonian simulation reduces the query complexity from \mathcal{O}\left(t\alpha +\frac{\log (1/ϵ)}{\log (e+\log (1/ϵ)/(\alpha |t|))}\right) [8,9] to \mathcal{O}\left(t\alpha \sqrt{\mathrm{\Delta }}{\log }^{3/2}\left(\frac{t\alpha }{ϵ}\right)+{\mathrm{\Delta }}^{-1/2}{\log }^{5/2}\left(\frac{t\alpha }{ϵ}\right)\right). Moreover, by regarding the block-encoding as state overlap, i.e.,

\frac{{\hat{H}}_{jk}}{\alpha }=⟨u_j|\frac{\hat{H}}{\alpha }|u_k⟩=\left(⟨0{|}_a⟨u_j{|}_s{\hat{U}}_{\mathrm{r}\mathrm{o}\mathrm{w}}^{†}\right)\left({\hat{U}}_{\mathrm{c}\mathrm{o}\mathrm{l}}|0{⟩}_a|u_k{⟩}_s\right)=⟨{\chi }_{0,j}|{\psi }_{0,k}{⟩}_{as},

Hamiltonian simulation can be accelerated in instances where the state-preparation procedure in amplitude amplification is costly. Instead of preparing |{\psi }_j⟩ and |{\chi }_j⟩, one may reduce the difficulty of state preparation by preparing |{\tilde{\psi }}_j⟩=\sqrt{{\lambda }_{\beta }{\beta }_j}|{\psi }_j⟩|0⟩+|{\psi }_{\mathrm{b}\mathrm{a}\mathrm{d}}⟩|1⟩ and |{\tilde{\chi }}_j⟩=\sqrt{{\lambda }_{\gamma }{\gamma }_j}|{\chi }_j⟩|0⟩+|{\chi }_{\mathrm{b}\mathrm{a}\mathrm{d}}j⟩|2⟩. Consequently, the normalization factor is reduced from the spectral norm to \alpha \sqrt{{\lambda }_{\gamma }{\lambda }_{\beta }}, which yields a speedup.

Haah et al. (2023) [52] found that for the evolution of lattice Hamiltonians, one can take advantage of the Lieb-Robinson bound [53] to achieve optimal gate complexity in Hamiltonian simulation. Specifically, if the terms in the Hamiltonian have supports that barely intersect, the entire Hamiltonian evolution can be split into evolving each term (or group of terms) and then compensating for the intersection terms by other local Hamiltonian evolutions. This method achieves the least gate complexity for local lattice Hamiltonian simulation to date.

Theorem 6.4 (Hamiltonian simulation in gate complexity [52]). Let H(t)=\sum _{X\subseteq \mathrm{\Lambda }}{h}_X(t) be a time-dependent Hamiltonian on a lattice \mathrm{\Lambda } of n qubits, embedded in the Euclidean metric space {\mathbb{R}}^D. Assume that every unit ball contains \mathcal{O}(1) qubits and h_X=0 if \mathrm{d}\mathrm{i}\mathrm{a}\mathrm{m}(X)>1. Also assume that every local term h_X(t) is efficiently computable (e.g., analytic), piecewise slowly varying on the time domain [0,T], and satisfies \Vert h_X(t)\Vert \le 1 for any X and t. Then, there exists a quantum algorithm that can approximate the time evolution of H for time T to accuracy ϵ using

\mathcal{O}\left(Tn\mathrm{p}\mathrm{o}\mathrm{l}\mathrm{y}\mathrm{l}\mathrm{o}\mathrm{g}(Tn/ϵ)\right),

2-qubit local gates, and has depth

\mathcal{O}\left(T\mathrm{p}\mathrm{o}\mathrm{l}\mathrm{y}\mathrm{l}\mathrm{o}\mathrm{g}(Tn/ϵ)\right).

[54] used generalized QSP to double the efficiency of Hamiltonian simulation. This improvement stems from the walk operator U being flipped to its inverse by moving the reflection before the block-encoding unitary V. Then, selecting between U and U^{†} results in R(2\theta ) instead of R(\theta ). This trick raises the issue of complex-valued functions for odd orders, which the original QSP cannot implement. However, generalized QSP, which performs general SU(2) rotations, accommodates this because the functions that GQSP implements are complex-valued. [55] used QSVT to simulate a random Hamiltonian to benchmark the ability of a quantum computer to perform Hamiltonian simulation.

A major problem in QSP is the classical computation of the angle vector \overrightarrow{\varphi }\in {\mathbb{R}}^N. Childs’s work [26] states in Appendix H.3 that the angle vector is uncomputable for degrees larger than 31. To circumvent this problem, one can segment the Hamiltonian simulation so that the order of the polynomials required for each small time step can be computed efficiently. However, QSP only exhibits superiority over other Hamiltonian simulation algorithms when the whole evolution is implemented in one unitary.

Many subsequent works have focused on resolving the classical computational hardness regarding this angle. [56] used a product decomposition to establish a classical algorithm with a random-access memory model to compute the polynomial in time \mathcal{O}(N^3\text{polylog}(N/ϵ)). The angles were found by identifying a Laurent polynomial from the input functions. In this step, coefficients with magnitudes that are too small are discarded because they cause instability in the angle-finding algorithm. The angles can then be computed efficiently using a Fast Fourier Transform. [57] improved this algorithm and obtained a similar cubic dependence on N through numerical observation.

Dong et al. [58] proposed an optimization-based method and achieved classical complexity of \mathcal{O}(N^2) through numerical observation. This optimization method heavily relies on the initial guess of the phase vector, as the landscape of the optimization space is complex and contains a large number of local minima. A heuristic initial guess was provided in the paper. Although it is numerically robust, a theoretical certification of its effectiveness is still lacking.

Later, generalized QSP [59] further reduced the complexity of finding phase angles by allowing for arbitrary SU(2) signal operators. Therefore, the phase angles can be obtained by solving a simple optimization problem, and the classical complexity is \mathcal{O}(N\log N), which is almost linear with the truncation order N.

Optimization of the oracle construction in the QSP framework is also important. Zhang et al. [60] improved the \mathrm{S}\mathrm{E}\mathrm{L}\mathrm{E}\mathrm{C}\mathrm{T} and \mathrm{P}\mathrm{R}\mathrm{E}\mathrm{P}\mathrm{A}\mathrm{R}\mathrm{E} oracles for qubitization by utilizing excess ancilla qubits and preparing the rotation angle in parallel. By using product unitary memory (PUM), the \mathrm{S}\mathrm{E}\mathrm{L}\mathrm{E}\mathrm{C}\mathrm{T} and \mathrm{P}\mathrm{R}\mathrm{E}\mathrm{P}\mathrm{A}\mathrm{R}\mathrm{E} oracles can be implemented with circuit depths of \mathcal{O}(\log (nL)) and \mathcal{O}(\log (L)), respectively. This achieves an exponential speedup in depth. Another parallel quantum algorithm for Hamiltonian simulation [61] exponentially improved the simulation precision and achieved a doubly (poly-)logarithmic circuit depth of \mathcal{O}(\text{polylog}\log (1/ϵ)).

Babbush et al. [62] discussed how to encode an electronic Hamiltonian into a quantum circuit with fewer gates than before, especially with a linear dependence on the T-gate count. The reduction is achieved by finding better ways to perform the \mathrm{S}\mathrm{E}\mathrm{L}\mathrm{E}\mathrm{C}\mathrm{T} and \mathrm{P}\mathrm{R}\mathrm{E}\mathrm{P}\mathrm{A}\mathrm{R}\mathrm{E} oracles, as some gates cancel each other. We, therefore, can implement fewer gates than faithfully construct all multi-control gates. Further, by introducing unary iteration and leveraging the fact that the electronic Hamiltonian has only four types of Pauli strings and \mathcal{O}(N) unique values of coefficients, significant improvements were made in simplifying these heavily used oracles. Ultimately, the circuit size and depth are linear with N for both Clifford and T gates.

Wan et al. [63] found that for arbitrary fermionic Hamiltonians, if each term in the Hamiltonian involves at most k=\mathcal{O}(1) spin-orbitals under the Jordan-Wigner transformation, the circuit from Babbush et al. (2018) [62] could be parallelized. The depth of the \mathrm{S}\mathrm{E}\mathrm{L}\mathrm{E}\mathrm{C}\mathrm{T} oracle can be exponentially reduced to \mathcal{O}(\text{poly}(\log (N))) using Clifford+T gates without using ancilla qubits, where N is the total number of spin-orbitals.

6.4 Randomization accelerated Hamiltonian simulation

Over the last few years, randomization has proven to be a potent tool for reducing costs in Hamiltonian simulation, revealing new ways to handle the complexity of quantum evolution. A particularly illustrative example comes from the juxtaposition of two methods proposed in 2019 that both rely on randomization but implement it quite differently [14,17]. We reviewed these in the section on variants of PF. Here we introduce randomized frameworks for post-Trotter methods.

Meister et al. [64] introduced an adaptive truncation scheme for LCU-based simulation. Rather than eliminating terms purely by polynomial order, their algorithm prioritizes high-magnitude terms in a truncated Taylor series, using a randomized iterative selection to maximize accuracy per gate. This method offers a constant-factor improvement.

Later, Wang and Zhao [18] brought randomization into the realm of truncated expansions in another way. Their Randomized Truncated Series (RTS) framework probabilistically mixes two truncated expansions of the time evolution operator, thereby achieving a continuously adjustable truncation order. RTS obtains a quadratic reduction in truncation error over conventional single-cutoff approaches. This approach is flexible and extends across multiple truncated-polynomial-based algorithms, including simulation algorithms.

In a closely related vein, Martyn and Rall [65] turned their attention specifically to QSP-based Hamiltonian simulation, proposing stochastic QSP. Here, instead of mixing two expansions, one randomizes over an entire ensemble of polynomial functions. Although narrower in scope—being tailored to the QSP framework—this approach can effectively halve the query complexity of QSP-based simulations. By incorporating probabilistic mixtures of polynomials, Martyn and Rall again illustrate the powerful role of randomness in slashing the cost of truncated-polynomial-based algorithms.

6.5 Experimental realization

Reference [66] addresses one of the central challenges in implementing Quantum Signal Processing (QSP) on noisy intermediate-scale quantum (NISQ) devices, namely the block-encoding of the target Hamiltonian. To alleviate the heavy circuit requirements typically associated with block-encoding, the authors propose two complementary strategies. First, they introduce a variational quantum eigensolver (VQE) ansatz to approximate the block-encoding for small system sizes, thus avoiding the need for exact (and often deep) construction of the block-encoded operator. Second, they exploit multiplexor circuit compilation to reduce the overhead of multi-controlled gates, substantially lowering the circuit depth. These methods were experimentally validated through Hamiltonian simulation of the Ising model with n=3 and n=4 qubits on the Quantinuum H1-1 trapped-ion quantum computer, where circuits involving up to 452 two-qubit gates were successfully executed. The authors also discuss how gate errors, decoherence, and other practical noise sources impact the overall fidelity, highlighting the potential of these techniques to render QSP more tractable on present-day quantum hardware.

7 Hamiltonian simulation via quantum singular value transformation

This section presents Hamiltonian simulation in the framework of quantum singular value transformation [9] (QSVT). Conceptually, QSVT subsumes qubitization and quantum signal processing (QSP) by formulating them as polynomial transformations of operators given through (projected) block-encodings. In addition to clarifying the mechanism behind optimal simulation algorithms, the QSVT formalism applies uniformly to more general operators than the original QSP setting, including non-Hermitian or rectangular contractions.

QSP implements a bounded polynomial transformation p(\lambda ) of the eigenvalues \lambda \in [-1,1] of a Hermitian operator, whereas QSVT applies the same idea to the singular values of a general operator A encoded as a projected unitary A=\tilde{\mathrm{\Pi }}U\mathrm{\Pi }; when A is Hermitian, in which case the singular values coincide with |\lambda |, QSVT reduces to the usual QSP eigenvalue transformation.

QSVT is implemented by alternating conditional phases with applications of U and U^{†}, as set up in Section 2. For a phase vector \mathrm{\Phi }=({\varphi }_0,{\varphi }_1,\dots ,{\varphi }_m), define an alternating sequence of the form

U_{\mathrm{\Phi }}:=e^{i{\varphi }_0(2\tilde{\mathrm{\Pi }}-I)}U{e}^{i{\varphi }_1(2\mathrm{\Pi }-I)}{U}^{†}{e}^{i{\varphi }_2(2\tilde{\mathrm{\Pi }}-I)}U\cdots ,

where the pattern alternates between \tilde{\mathrm{\Pi }}- and \mathrm{\Pi }-controlled phases and alternates U and U^{†}. Such a sequence can be implemented using O(m) uses of U and U^{†} and O(m) uses of the subspace-checking primitives for \mathrm{\Pi } and \tilde{\mathrm{\Pi }}, plus O(m) single-qubit phase gates.

7.1 Hamiltonian simulation from QSVT

The central statement behind QSVT is that the above alternating sequence can implement a polynomial transformation of the singular values of the encoded operator A=\tilde{\mathrm{\Pi }}U\mathrm{\Pi }. In the special case where A is Hermitian, singular values coincide with absolute eigenvalues, and the transformation may be phrased as an eigenvalue polynomial.

Informally, let P be a real polynomial of degree d satisfying the QSVT admissibility conditions: (i) |P(x)|\le 1 for all x\in [-1,1]; (ii) the parity of P matches the parity required by the alternating construction (even or odd, depending on whether the sequence starts and ends on \mathrm{\Pi } or \tilde{\mathrm{\Pi }}). Then there exists a phase vector \mathrm{\Phi } of length O(d) such that the corresponding U_{\mathrm{\Phi }} implements a new projected encoding whose effective operator is P(A), up to a controllable approximation error:

\tilde{\mathrm{\Pi }}{U}_{\mathrm{\Phi }}\mathrm{\Pi }\approx P(A).

Equivalently, in the standard block-encoding setting, U_{\mathrm{\Phi }} yields a block-encoding of P(H/\alpha ) provided one starts from a block-encoding of H/\alpha.

Corollary 7.1 (Complexity of block-Hamiltonian simulation [9]). For any ϵ\in (0,1/2) and t\in \mathbb{R}, it is necessary and sufficient (up to constant factors) to use the block-encoding U a total number of times

\mathrm{\Theta }\left(\alpha |t|+\frac{\log (1/ϵ)}{\log (e+\frac{\log (1/ϵ)}{\alpha |t|})}\right),

to implement such a block-encoding of e^{\mathrm{i}tH}.

The QSVT construction reduces Hamiltonian simulation to implementing a bounded polynomial approximation of

x\mapsto e^{i(\alpha t)x}\text{on}x\in [-1,1],

inside the invariant \mathrm{S}\mathrm{U}(2) subspaces induced by the qubitization framework. A convenient way to design such approximants is to expand \cos (\alpha tx) and \sin (\alpha tx) in a Chebyshev series and truncate at some degree. The truncation error is

ϵ\le \mathcal{O}\left({\frac{\alpha |t|}{r}}^r\right),

where r is the polynomial degree.

It is therefore natural to introduce a degree parameter r(t,ϵ) as the smallest integer r for which the truncation error is at most ϵ. In the heuristic model r\log (r/t)\approx \log (1/ϵ), one obtains the Lambert-W expression

r(t,ϵ)\approx \frac{\log (1/ϵ)}{W(\log (1/ϵ)/t)}.

However, this expression fails for some ranges of t. [9] developed a revised bound that is valid for all t:

r(\tau ,ϵ)=\mathrm{\Theta }\left(\tau +\frac{\log (1/ϵ)}{\log (e+\frac{\log (1/ϵ)}{\tau })}\right).

This cleanly interpolates between two regimes. In the long-time regime (t\gg \log (1/ϵ)), the degree is linear: r(t,ϵ)=\mathrm{\Theta }(t). In the short-time regime (t\ll \log (1/ϵ)), the degree is sublinear in \log (1/ϵ), i.e., r(t,ϵ)=\mathrm{\Theta }\left(\frac{\log (1/ϵ)}{\log (\log (1/ϵ)/t)}\right).

For Hamiltonian simulation in the block-encoding model, the construction shows that an ϵ-precise block-encoding of e^{\mathrm{i}tH} can be implemented using 3r\left(\frac{e\alpha |t|}{2},\frac{ϵ}{6}\right) uses of U and U^{†}, see Theorem 58 in [9].

8 Theoretical advantages and limitation of Hamiltonian simulation

While the previous sections have demonstrated the efficiency of various quantum algorithms for Hamiltonian simulation, two fundamental questions are whether classical computers can solve this problem efficiently and whether there are limitations on the quantum efficiency of Hamiltonian simulation. In this section, we will first review that Hamiltonian simulation is classically hard unless quantum computation has no super-polynomial advantage over classical computation [1,67]. And then we will review that quantum algorithms for Hamiltonian simulation cannot be fast-forwarded, i.e., any generic algorithm cannot run in sublinear time with respect to the actual evolution time [22].

8.1 Advantages: BQP-completeness

Along with the first proposal of quantum simulation [1], Feynman proposed the concept of quantum computers as well as the circuit-to-Hamiltonian construction [67], which established the equivalence between the quantum circuit model and the Hamiltonian dynamics. This result can be stated in the modern complexity theory language as:

HamiltoniandynamicssimulationisBQP-complete.

Informally, \mathsf{\text{BQP}} is the class of problems that can be solved with bounded error by quantum algorithms in polynomial time, meaning that they are tractable for quantum computers [68–71]. A \mathsf{\text{BQP}}-complete problem is among the hardest problems in \mathsf{\text{BQP}} in the sense that every problem in \mathsf{\text{BQP}} can be reduced to it efficiently. Assuming quantum computers are more powerful than classical computers, a \mathsf{\text{BQP}}-complete problem is classically intractable and could therefore demonstrate quantum advantage. Here, we briefly review Feynman’s circuit-to-Hamiltonian construction, which is a fundamental tool for encoding a quantum circuit into a local Hamiltonian.

Lemma 8.1 (Circuit-to-Hamiltonian construction, rephrased from [67,72–73]). Given a quantum circuit U=U_N\dots U_1 composed of N gates, there exists a Hamiltonian H_U\in {\mathcal{H}}_{\mathbf{c}}\otimes {\mathcal{H}}_{\mathbf{o}} acting on a clock register and an output register, namely,

H_U=\sum _{j=1}^N\sqrt{j(N-j+1)}\left(|j⟩{⟨j-1|}_{\mathbf{c}}\otimes U_j+\mathit{\text{H. C.}}\right),

such that the Hamiltonian dynamics governed by H_U at t=\pi /2 reproduce the original quantum-circuit output as

e^{-\mathrm{i}{H}_{U}t}|0^{\otimes \log (N)}⟩{}_{\mathbf{c}}{|{\psi }_0⟩}_{\mathbf{o}}={|N⟩}_{\mathbf{c}}\otimes U{|{\psi }_0⟩}_{\mathbf{o}},

on the output register.

For concrete examples, we refer the reader to pedagogical references [74–75] or specific references [72–73] on perfect state transfer. In summary, any polynomial-size quantum circuit can be encoded in a local Hamiltonian evolution with a polynomial evolution time. This result has profound implications for quantum computation and lays the foundation for quantum advantage arising from Hamiltonian simulation. In addition, there have been recent works [76–80] showing the rigorous advantage of Hamiltonian simulation from the perspective of sampling complexity.

8.2 Limitations: No fast-forwarding theorem

While we have efficient quantum algorithms for Hamiltonian simulation, there are fundamental limits (lower bounds) on the complexity of these quantum algorithms. The No fast-forwarding theorem [22] is such a lower bound on the gate complexity of quantum dynamics simulation in terms of evolution time t. It states that the gate complexity of any quantum algorithm to simulate the evolution of a local Hamiltonian for time t is at least linear with t.

Theorem 8.1 (No fast-forwarding theorem [22]). For all positive integers N\le 2^n, there exists a (row-computable) 2-sparse n-qubit Hamiltonian H_{\mathbf{x}} such that simulating the evolution of H_{\mathbf{x}} for scaled time N\pi /2 within precision (trace distance) 1/4 requires at least N/4 queries to H_{\mathbf{x}}.

The proof starts from the lower bound for the query complexity of Hamiltonian simulation. To prove it, Ref. [22] made a reduction from the PARITY problem to Hamiltonian simulation. The PARITY problem PARITY{}_N(\mathbf{x}) is to evaluate the parity of an input bit-string \mathbf{x} of length N by querying the bits as few as possible. It has been proved by Refs. [81–82] that any quantum algorithm requires \mathrm{\Omega }(N) queries to evaluate the parity, even allowing bounded error. One can leverage the circuit-to-Hamiltonian construction in Lemma 8.1 again to construct a sparse Hamiltonian H_{\mathbf{x}} for any input bit-string \mathbf{x} such that

e^{-\mathrm{i}{H}_{\mathbf{x}}\pi /2}|0^{\otimes n}⟩{}_{\mathbf{c}}{|0⟩}_{\mathbf{o}}\mapsto {|N⟩}_{\mathbf{c}}\otimes |{\text{PARITY}}_N(\mathbf{x})⟩{}_{\mathbf{o}},

where the output register encodes the parity of the input bit-string. In this way, if there exists a generic algorithm for Hamiltonian simulation with o(t) query complexity, one can solve PARITY{}_N(\mathbf{x}) by o(N) queries to the bit-string \mathbf{x}\in {0,1}^N. This leads to a contradiction with the known lower bound of parity. Together with the upper bound on the complexity by the high-order product formula and the lower bound by \mathrm{N}\mathrm{o}\mathrm{f}\mathrm{a}\mathrm{s}\mathrm{t}-\mathrm{f}\mathrm{o}\mathrm{r}\mathrm{w}\mathrm{a}\mathrm{r}\mathrm{d}\mathrm{i}\mathrm{n}\mathrm{g}, [22] showed that the product formula algorithms are nearly optimal in the parameter time t.

With ideas from the proof of the \mathrm{N}\mathrm{o}\mathrm{f}\mathrm{a}\mathrm{s}\mathrm{t}-\mathrm{f}\mathrm{o}\mathrm{r}\mathrm{w}\mathrm{a}\mathrm{r}\mathrm{d}\mathrm{i}\mathrm{n}\mathrm{g} theorem, Refs. [46,83] proved that any sparse Hamiltonian simulation method must use \mathrm{\Omega }\left(t+\frac{\log (1/ϵ)}{\log \log (1/ϵ)}\right) discrete queries to obtain error at most ϵ, so the dependence of the query complexity of post-Trotter algorithms (LCU, also QSP) on ϵ is tight up to constant factors, which is exponentially better than the Trotter algorithms.

In addition to the results in [22,46,83] based on reductions from PARITY, Haah et al. [52] proved a similar lower bound from the perspective of gate synthesis. In addition, Atia and Aharonov [84] showed the impossibility of a generic fast-forwarding procedure for realizable Hamiltonians from the perspective of complexity classes: if the 2-sparse row-computable Hamiltonians can be exponentially fast-forwarded, the complexity class \mathsf{\text{PSPACE}} equals \mathsf{\text{BQP}} which is highly unlikely.

Nevertheless, these lower bounds do not rule out the possibilities of Hamiltonian simulation with large but “low-depth” circuits by running things in parallel. Chia et al. [85] closed this gap by showing that sparse Hamiltonians and (geometrically) local Hamiltonians cannot be fast-forwarded in parallel.

On the other hand, the No fast-forwarding theorem does not directly cover special cases, such as simulating low-energy states. For this case, [86] proved that even though the Trotter simulation is restricted to a low-energy subspace, a similar No fast-forwarding theorem holds. As a complement, Zlokapa and Somma [87] studied post-Trotter algorithms, including LCU and QSP within low-energy subspace. They also proved the No fast-forwarding theorem of low-energy states with the corresponding access models.

While the No fast-forwarding theorem sets the fundamental limit on simulation time of general Hamiltonians, Atia and Aharonov [84] showed that commuting local Hamiltonians and quadratic fermionic Hamiltonians can be exponentially fast-forwarded by quantum algorithms. Later, Gu, Somma, and Şahinoğlu [88] coined the definition of “fast-forwarding” that accounts for any asymptotic complexity improvement of the general case, going beyond the exponential fast-forwarding studied previously.

In a similar manner, Feng et al. gave a lower bound on the communication complexity of distributed quantum simulation by reducing from the INNERPRODUCT problem to distributed quantum simulation [89].

Theorem 8.2 (No fast-forwarding theorem for distributed quantum simulation [89]). For any positive integer \tau \le n, there exists a sparse Hamiltonian H acting on 2n+\log (n)+1 qubits with \Vert H\Vert =\tau and a constant evolution time t=\pi /2 such that the (bounded-error) quantum communication complexity ({\mathsf{\text{Qcc}}}_2) of any generic \mathrm{\Gamma }-partite protocol for distributed quantum simulation scales at least linearly in scaled evolution time, that is, {\mathsf{\text{Qcc}}}_2=\mathrm{\Omega }(\mathrm{\Gamma }\Vert H\Vert t).

The proof idea is as follows. The distributed version of PARITY, namely INNERPRODUCT, is a Boolean function denoted I{P}_n(\mathbf{x},\mathbf{y}):{0,1}^n\times {0,1}^n\to 0,1 that evaluates the inner product of two bit-strings. Refs. [90–91] proved that quantum communication complexity {\mathsf{\text{Qcc}}}_2 of evaluating I{P}_n is {\mathsf{\text{Qcc}}}_2(I{P}_n)=\mathrm{\Omega }(n). If two parties could simulate quantum dynamics distributively with o(\Vert H\Vert t) quantum communication, then the Boolean function I{P}_n can be evaluated with o(n) communication, which would contradict the known lower bound.

For a \mathrm{\Gamma }-partite network, [92] proved that the lower bound of \mathrm{\Gamma }-partite INNERPRODUCT can be extended to \mathrm{\Omega }(\mathrm{\Gamma }n). So, the proof can be directly extended to the \mathrm{\Gamma }-partite case by replacing the Toffoli gates with \mathrm{\Gamma }-partite controlled-NOT gates, and the lower bound becomes \mathrm{\Omega }(\mathrm{\Gamma }\Vert H\Vert t). Since the quantum communication complexity of our post-Trotter protocol is linear with t and \mathrm{\Gamma }, this theorem suggests that our distributed quantum simulation algorithms achieve optimal dependence on evolution time t and the number of partitions \mathrm{\Gamma } provided \Vert H\Vert t\le n.

9 Discussion and conclusions

This review has provided a comprehensive historical and technical overview of the development of fault-tolerant digital Hamiltonian simulation. We discussed the evolution from early product-formula approaches to modern post-Trotter techniques such as the truncated Taylor-series and quantum signal-processing methods [2–3,5,8,93]. Together, these advances define the current landscape of digital quantum simulation and continue to shape the search for resource-optimal algorithms under realistic fault-tolerant constraints.

A central goal in this pursuit is to derive tighter upper bounds on simulation error, thereby reducing overall gate complexity (resource cost) and improving scalability. As reviewed above, recent progress has demonstrated multiple promising strategies in this direction. For example, Watson’s extrapolation-based approach achieves an exponentially reduced-depth product formula [43], representing a major conceptual advance in Hamiltonian-simulation design. Beyond circuit depth and gate cost, quantum communication complexity is emerging as an equally critical factor for large-scale or distributed implementations of quantum computing. New methods for reducing communication overhead are enabling more scalable distributed quantum simulation architectures [89].

Although fully fault-tolerant implementations of these algorithms remain beyond current hardware capabilities, there is increasing progress toward near-term realizations. Approaches such as variational quantum simulation [94], perturbative quantum simulation [95], shadow-based Hamiltonian simulation [96], and error-mitigated digital simulation [42,97] exemplify how hybrid and error-aware paradigms can extend simulation capabilities on noisy intermediate-scale quantum devices.

Looking ahead, quantum simulation is poised to become a cornerstone for scientific discovery in chemistry, materials science, and condensed-matter physics [98–99]. As algorithmic theory continues to mature and quantum hardware steadily improves, the synergy between fault-tolerant design principles and application-oriented methodologies will be essential for translating formal quantum-algorithmic advances into practical, real-world simulations. In this sense, the ongoing progress in Hamiltonian simulation not only refines our understanding of quantum computational complexity but also marks a decisive step toward realizing the broader vision of quantum-enabled scientific computing.

References

[1] Feynman R P . Simulating physics with computers. International Journal of Theoretical Physics, 1982, 21( 6−7): 467–488

[2] Lloyd S . Universal quantum simulators. Science, 1996, 273( 5278): 1073–1078

[3] Childs A M, Su Y, Tran M C, Wiebe N, Zhu S . Theory of trotter error with commutator scaling. Physical Review X, 2021, 11( 1): 011020

[4] Suzuki M . Quantum Monte Carlo methods — recent developments. Physica A: Statistical Mechanics and its Applications, 1993, 194( 1−4): 432–449

[5] Berry D W, Childs A M, Cleve R, Kothari R, Somma R D . Simulating Hamiltonian dynamics with a truncated Taylor series. Physical Review Letters, 2015, 114( 9): 090502

[6] Long G L . Duality quantum computing and duality quantum information processing. International Journal of Theoretical Physics, 2011, 50( 4): 1305–1318

[7] Childs A M, Wiebe N . Hamiltonian simulation using linear combinations of unitary operations. Quantum Information and Computation, 2012, 12( 11−12): 901–924

[8] Low G H, Chuang I L . Hamiltonian simulation by qubitization. Quantum, 2019, 3: 163

[9] Gilyén A, Su Y, Low G H, Wiebe N. Quantum singular value transformation and beyond: exponential improvements for quantum matrix arithmetics. In: Proceedings of the 51st Annual ACM SIGACT Symposium on Theory of Computing. 2019, 193−204

[10] Martyn J M, Rossi Z M, Tan A K, Chuang I L. A grand unification of quantum algorithms. 2021, arXiv preprint arXiv: 2105.02859

[11] Berry D W, Wan K, Baczewski A D, Eklund E C, Tikku A, Babbush R. Quantum simulation of chemistry via quantum fast multipole method. 2025, arXiv preprint arXiv: 2510.07380v1

[12] Berry D W, Rubin N C, Elnabawy A O, Ahlers G, DePrince III A E, Lee J, Gogolin C, Babbush R . Quantum simulation of realistic materials in first quantization using non-local pseudopotentials. npj Quantum Information, 2024, 10( 1): 130

[13] Ballarin M, Cataldi G, Magnifico G, Jaschke D, Liberto M D, Siloi I, Montangero S, Silvi P . Digital quantum simulation of lattice fermion theories with local encoding. Quantum, 2024, 8: 1460

[14] Campbell E . Random compiler for fast Hamiltonian simulation. Physical Review Letters, 2019, 123( 7): 070503

[15] Campbell E . Shorter gate sequences for quantum computing by mixing unitaries. Physical Review A, 2017, 95( 4): 042306

[16] Hastings M B . Turning gate synthesis errors into incoherent errors. Quantum Information and Computation, 2017, 17( 5−6): 488–494

[17] Childs A M, Ostrander A, Su Y . Faster quantum simulation by randomization. Quantum, 2019, 3: 182

[18] Wang Y, Zhao Q. Faster quantum algorithms with “fractional”-truncated series. 2024, arXiv preprint arXiv: 2402.05595v2

[19] Childs A M, Berry D W . Black-box Hamiltonian simulation and unitary implementation. Quantum Information and Computation, 2012, 12( 1−2): 29–62

[20] Aharonov D, Ta-Shma A. Adiabatic quantum state generation and statistical zero knowledge. In: Proceedings of the 35th Annual ACM Symposium on Theory of Computing. 2003, 20−29

[21] Childs A M, Goldstone J . Spatial search by quantum walk. Physical Review A, 2004, 70( 2): 022314

[22] Berry D W, Ahokas G, Cleve R, Sanders B C . Efficient quantum algorithms for simulating sparse Hamiltonians. Communications in Mathematical Physics, 2007, 270( 2): 359–371

[23] Suzuki M . General theory of fractal path integrals with applications to many-body theories and statistical physics. Journal of Mathematical Physics, 1991, 32( 2): 400–407

[24] Somma R D . A Trotter-Suzuki approximation for Lie groups with applications to Hamiltonian simulation. Journal of Mathematical Physics, 2016, 57( 6): 062202

[25] Childs A M, Su Y . Nearly optimal lattice simulation by product formulas. Physical Review Letters, 2019, 123( 5): 050503

[26] Childs A M, Maslov D, Nam Y, Ross N J, Su Y . Toward the first quantum simulation with quantum speedup. Proceedings of the National Academy of Sciences of the United States of America, 2018, 115( 38): 9456–9461

[27] Zhao Q, Zhou Y, Shaw A F, Li T, Childs A M . Hamiltonian simulation with random inputs. Physical Review Letters, 2022, 129( 27): 270502

[28] Chen C F, Brandão F G S L . Average-case speedup for product formulas. Communications in Mathematical Physics, 2024, 405( 2): 32

[29] Şahinoğlu B, Somma R D . Hamiltonian simulation in the low-energy subspace. npj Quantum Information, 2021, 7( 1): 119

[30] Zhao Q, Zhou Y, Childs A M . Entanglement accelerates quantum simulation. Nature Physics, 2025, 21( 8): 1338–1345

[31] Heyl M, Hauke P, Zoller P . Quantum localization bounds trotter errors in digital quantum simulation. Science Advances, 2019, 5( 4): eaau8342

[32] Sieberer L M, Olsacher T, Elben A, Heyl M, Hauke P, Haake F, Zoller P . Digital quantum simulation, Trotter errors, and quantum chaos of the kicked top. npj Quantum Information, 2019, 5( 1): 78

[33] Trivedi R, Franco Rubio A, Cirac J I . Quantum advantage and stability to errors in analogue quantum simulators. Nature Communications, 2024, 15( 1): 6507

[34] Yu W, Xu J, Zhao Q . Observable-driven speed-ups in quantum simulations. Communications Physics, 2025, 8( 1): 340

[35] Feng T, Cao Y, Zhao Q. Trotterization, operator scrambling, and entanglement. 2025, arXiv preprint arXiv: 2506.23345

[36] Nakaji K, Bagherimehrab M, Aspuru-Guzik A . High-order randomized compiler for Hamiltonian simulation. PRX Quantum, 2024, 5( 2): 020330

[37] Granet E, Dreyer H . Hamiltonian dynamics on digital quantum computers without discretization error. npj Quantum Information, 2024, 10( 1): 82

[38] Low G H, Kliuchnikov V, Wiebe N. Well-conditioned multiproduct Hamiltonian simulation. 2019, arXiv preprint arXiv: 1907.11679

[39] Carrera Vazquez A, Egger D J, Ochsner D, Woerner S . Well-conditioned multi-product formulas for hardware-friendly Hamiltonian simulation. Quantum, 2023, 7: 1067

[40] Zhuk S, Robertson N F, Bravyi S . Trotter error bounds and dynamic multi-product formulas for Hamiltonian simulation. Physical Review Research, 2024, 6( 3): 033309

[41] Zeng P, Sun J, Jiang L, Zhao Q . Simple and high-precision Hamiltonian simulation by compensating trotter error with linear combination of unitary operations. PRX Quantum, 2025, 6( 1): 010359

[42] Endo S, Zhao Q, Li Y, Benjamin S, Yuan X . Mitigating algorithmic errors in a Hamiltonian simulation. Physical Review A, 2019, 99( 1): 012334

[43] Watson J D, Watkins J . Exponentially reduced circuit depths using trotter error mitigation. PRX Quantum, 2025, 6( 3): 030325

[44] Watson J D. Randomly compiled quantum simulation with exponentially reduced circuit depths. 2025, arXiv preprint arXiv: 2411.04240

[45] Low G H, Chuang I L . Optimal Hamiltonian simulation by quantum signal processing. Physical Review Letters, 2017, 118( 1): 010501

[46] Berry D W, Childs A M, Kothari R. Hamiltonian simulation with nearly optimal dependence on all parameters. In: Proceedings of the 56th Annual Symposium on Foundations of Computer Science. 2015, 792−809

[47] Shende V V, Bullock S S, Markov I L . Synthesis of quantum-logic circuits. IEEE Transactions on Computer-Aided Design of Integrated Circuits and Systems, 2006, 25( 6): 1000–1010

[48] Zhao Q, Yuan X . Exploiting anticommutation in Hamiltonian simulation. Quantum, 2021, 5: 534

[49] Low G H, Yoder T J, Chuang I L . Methodology of resonant equiangular composite quantum gates. Physical Review X, 2016, 6( 4): 041067

[50] Childs A M, Cleve R, Deotto E, Farhi E, Gutmann S, Spielman D A. Exponential algorithmic speedup by a quantum walk. In: Proceedings of the 35th Annual ACM Symposium on Theory of Computing. 2003, 59−68

[51] Low G H, Chuang I L. Hamiltonian simulation by uniform spectral amplification. 2017, arXiv preprint arXiv: 1707.05391

[52] Haah J, Hastings M B, Kothari R, Low G H . Quantum algorithm for simulating real time evolution of lattice Hamiltonians. SIAM Journal on Computing, 2023, 52( 6): FOCS18-250–FOCS18-284

[53] Lieb E H, Robinson D W . The finite group velocity of quantum spin systems. Communications in Mathematical Physics, 1972, 28( 3): 251–257

[54] Berry D W, Motlagh D, Pantaleoni G, Wiebe N . Doubling the efficiency of Hamiltonian simulation via generalized quantum signal processing. Physical Review A, 2024, 110( 1): 012612

[55] Dong Y, Whaley K B, Lin L . A quantum Hamiltonian simulation benchmark. npj Quantum Information, 2022, 8( 1): 131

[56] Haah J . Product decomposition of periodic functions in quantum signal processing. Quantum, 2019, 3: 190

[57] Chao R, Ding D, Gilyen A, Huang C, Szegedy M. Finding angles for quantum signal processing with machine precision. 2020, arXiv preprint arXiv: 2003.02831v2

[58] Dong Y, Meng X, Whaley K B, Lin L . Efficient phase-factor evaluation in quantum signal processing. Physical Review A, 2021, 103( 4): 042419

[59] Motlagh D, Wiebe N . Generalized quantum signal processing. PRX Quantum, 2024, 5: 020368

[60] Zhang X M, Li T, Yuan X . Quantum state preparation with optimal circuit depth: implementations and applications. Physical Review Letters, 2022, 129( 23): 230504

[61] Zhang Z, Wang Q, Ying M . Parallel quantum algorithm for Hamiltonian simulation. Quantum, 2024, 8: 1228

[62] Babbush R, Gidney C, Berry D W, Wiebe N, McClean J, Paler A, Fowler A, Neven H . Encoding electronic spectra in quantum circuits with linear T complexity. Physical Review X, 2018, 8( 4): 041015

[63] Wan K . Exponentially faster implementations of Select(H) for fermionic Hamiltonians. Quantum, 2021, 5: 380

[64] Meister R, Benjamin S C, Campbell E T . Tailoring term truncations for electronic structure calculations using a linear combination of unitaries. Quantum, 2022, 6: 637

[65] Martyn J M, Rall P . Halving the cost of quantum algorithms with randomization. npj Quantum Information, 2025, 11( 1): 47

[66] Kikuchi Y, Mc Keever C, Coopmans L, Lubasch M, Benedetti M . Realization of quantum signal processing on a noisy quantum computer. npj Quantum Information, 2023, 9( 1): 93

[67] Feynman R P . Quantum mechanical computers. Optics News, 1985, 11( 2): 11–20

[68] Bernstein E, Vazirani U . Quantum complexity theory. SIAM Journal on Computing, 1997, 26( 5): 1411–1473

[69] Kitaev A Y, Shen A H, Vyalyi M N. Classical and Quantum Computation. Providence: American Mathematical Society, 2002

[70] Wocjan P, Zhang S. Several natural BQP-complete problems. 2006, arXiv preprint arXiv: quant-ph/0606179

[71] Arora S, Barak B. Computational Complexity: A Modern Approach. New York: Cambridge University Press, 2009

[72] Christandl M, Datta N, Ekert A, Landahl A J . Perfect state transfer in quantum spin networks. Physical Review Letters, 2004, 92( 18): 187902

[73] Kay A . Perfect, efficient, state transfer and its application as a constructive tool. International Journal of Quantum Information, 2010, 8( 4): 641–676

[74] Manenti R, Motta M. Quantum Information Science. Oxford: Oxford University Press, 2023

[75] Childs A M. Lecture Notes on Quantum Algorithms. University of Maryland, 2025. Available at: https://www.cs.umd.edu/~amchilds/qa/qa.pdf

[76] Hangleiter D, Kliesch M, Schwarz M, Eisert J . Direct certification of a class of quantum simulations. Quantum Science and Technology, 2017, 2( 1): 015004

[77] Bermejo-Vega J, Hangleiter D, Schwarz M, Raussendorf R, Eisert J . Architectures for quantum simulation showing a quantum speedup. Physical Review X, 2018, 8( 2): 021010

[78] Haferkamp J, Hangleiter D, Bouland A, Fefferman B, Eisert J, Bermejo-Vega J . Closing gaps of a quantum advantage with short-time Hamiltonian dynamics. Physical Review Letters, 2020, 125( 25): 250501

[79] Liu Z, Devulapalli D, Hangleiter D, Liu Y K, Kollár A J, Gorshkov A V, Childs A M . Efficiently verifiable quantum advantage on near-term analog quantum simulators. PRX Quantum, 2025, 6( 1): 010341

[80] Quek Y. Quantum advantage from random geometrically-two-local Hamiltonian dynamics. 2025, arXiv preprint arXiv: 2510.06321

[81] Farhi E, Goldstone J, Gutmann S, Sipser M . Limit on the speed of quantum computation in determining parity. Physical Review Letters, 1998, 81( 24): 5442

[82] Beals R, Buhrman H, Cleve R, Mosca M, de Wolf R . Quantum lower bounds by polynomials. Journal of the ACM, 2001, 48( 4): 778–797

[83] Berry D W, Childs A M, Cleve R, Kothari R, Somma R D. Exponential improvement in precision for simulating sparse Hamiltonians. In: Proceedings of the 46th Annual ACM Symposium on Theory of Computing. 2014, 283−292

[84] Atia Y, Aharonov D . Fast-forwarding of Hamiltonians and exponentially precise measurements. Nature Communications, 2017, 8( 1): 1572

[85] Chia N H, Chung K M, Hsieh Y C, Lin H H, Lin Y T, Shen Y C. On the impossibility of general parallel fast-forwarding of Hamiltonian simulation. In: Proceedings of the 38th Computational Complexity Conference. 2023, 33

[86] Gong W, Zhou S, Li T . Complexity of digital quantum simulation in the low-energy subspace: applications and a lower bound. Quantum, 2024, 8: 1409

[87] Zlokapa A, Somma R D . Hamiltonian simulation for low-energy states with optimal time dependence. Quantum, 2024, 8: 1449

[88] Gu S, Somma R D, Şahinoğlu B . Fast-forwarding quantum evolution. Quantum, 2021, 5: 577

[89] Feng T, Xu J, Yu W, Ye Z, Yao P, Zhao Q. Distributed quantum simulation. 2024, arXiv preprint arXiv: 2411.02881

[90] Kremer I. Quantum communication. Hebrew University, Dissertation, 1995

[91] de Wolf R. Quantum computing and communication complexity. University of Amsterdam, Dissertation, 2001

[92] Le Gall F, Suruga D. Bounds on oblivious multiparty quantum communication complexity. In: Proceedings of the 15th Latin American Symposium on LATIN 2022: Theoretical Informatics. 2022, 641−657

[93] Low G H. Hamiltonian simulation with nearly optimal dependence on spectral norm. In: Proceedings of the 51st Annual ACM SIGACT Symposium on Theory of Computing. 2019, 491−502

[94] Yuan X, Endo S, Zhao Q, Li Y, Benjamin S C . Theory of variational quantum simulation. Quantum, 2019, 3: 191

[95] Sun J, Endo S, Lin H, Hayden P, Vedral V, Yuan X . Perturbative quantum simulation. Physical Review Letters, 2022, 129( 12): 120505

[96] Somma R D, King R, Kothari R, O’Brien T E, Babbush R . Shadow Hamiltonian simulation. Nature Communications, 2025, 16( 1): 2690

[97] Rendon G, Watkins J, Wiebe N . Improved accuracy for trotter simulations using Chebyshev interpolation. Quantum, 2024, 8: 1266

[98] Zhang Y, Zhang X, Sun J, Lin H, Huang Y, Lv D, Yuan X. Fault-tolerant quantum algorithms for quantum molecular systems: a survey. 2025, arXiv preprint arXiv: 2502.02139

[99] McArdle S, Endo S, Aspuru-Guzik A, Benjamin S C, Yuan X . Quantum computational chemistry. Reviews of Modern Physics, 2020, 92( 1): 015003

Rights & permissions

The Author(s) 2026. This article is published with open access at link.springer.com and journal.hep.com.cn