Simple and high-precision Hamiltonian simulation by compensating Trotter error with linear combination of unitary operations

The major idea of the proposed Trotter-LCU algorithm is illustrated in Fig. 1. In a normal $K$th-order Trotter circuit, we decompose the time evolution $U(t)=e^{-iHt}$ to $\nu$ segments, each with a small evolution time $x=t/\nu$. For consistency, we denote the $0$th-order Trotter formula as $S_{0}(x)=I$. After we perform the $K$th-order Trotter formula $S_{K}(x)(K=0,1,2k,k\in\mathbb{N}_{+})$, there will be a remaining Trotter error $V_{K}(x):=U(x)S_{K}(x)^{\dagger}$, which affects the simulation accuracy. To address this problem, we introduce a random LCU formula to compensate the Trotter error using one ancilla and simple gates, which achieves a high-precision Hamiltonian simulation with low cost.

Consider the following LCU formula of an operator $V$,

such that the spectral norm distance $\|V-\tilde{V}\|\leq\varepsilon$. Here, $\mu>0$ is the $1$-norm (i.e., $l_{1}$-norm) of the coefficient vector, $\Pr(i)$ is a probability distribution over different unitaries $V_{i}$, and $\{V_{i}\}_{i=0}^{\Gamma-1}$ is a set of unitaries. There are usually two ways to implement the LCU formula: the coherent implementation Berry et al. (2014, 2015) and the random-sampling implementation Yang et al. (2021); Wan et al. (2022); Faehrmann et al. (2022). Our major focus is on the random-sampling implementation, where we can estimate the properties of the target state $\rho=U(t)\rho_{0}U(t)^{\dagger}$ with only one ancillary qubit. In Appendix H, we discuss the potential use of the coherent implementation of our algorithm. In the random-sampling implementation, we can use Eq. 1 to estimate an arbitrary observable value $O$ on $\rho$,

As is shown in Fig. 2(b), since the estimation of $\mathrm{Tr}(OV_{i}\rho_{0}V_{j}^{\dagger})$ can be implemented using Hadamard-test-type circuits Kitaev (1995), we only need to sample $V_{i}$ and $V_{j}$ based on the LCU formula in Eq. 1 to estimate $\mathrm{Tr}(O\rho)$ with $\varepsilon$ accuracy using $\mathcal{\mathcal{missing}}{O}(\mu^{4}/\varepsilon^{2})$ sampling resource, which owns an extra $\mu^{4}$ overhead compared to the normal Hamiltonian simulation algorithms Faehrmann et al. (2022). To make the algorithm efficient, we need to set $\mu$ to be a constant. We also provide a variant in Fig. 2(c) where the ancillary qubit is measured and reset in each segment, which is equivalent to Fig. 2(b) for the observable estimation. In this case, the expectation value of $\mu^{2}X^{(tot)}_{A}\otimes O_{S}$ provides an unbiased estimation of $\mathrm{Tr}(O\rho)$ where $X^{(tot)}:=\bigotimes_{k=1}^{\nu}X_{k}$ is the multiplication of all the ancillary measurement values. This variant reduces the need to store the ancilla qubit, simplifying the implementation on a fault-tolerant quantum computer.

The construction of an appropriate LCU formula for the $K$th-order Trotter remainder, $V_{K}(x)$, is crucial for developing an efficient Hamiltonian simulation algorithm. In the following subsections, we briefly introduce two approaches for constructing the LCU formula for $V_{K}(x)$. The resulting composite Trotter-LCU algorithms are referred to as paired Taylor-series compensation (PTSC) and nested-commutator compensation (NCC), respectively. Detailed analysis and performance proofs for these two algorithms can be found in Sec. IV and Sec. V.

Without loss of generality, we focus on the case of an $n$-qubit Hamiltonian $H$, which can be written as

where $\{P_{l}\}_{l}$ are different $n$-qubit Pauli operators. We set all the coefficients $\{\alpha_{l}\}_{l}$ to be positive and absorb the signs into Pauli operators $\{P_{l}\}_{l}$. $\lambda:=\sum_{l}\alpha_{l}$ is the $l_{1}$-norm of the Hamiltonian coefficient vector. We consider the Hamiltonians where $\lambda$ increases polynomially with respect to $n$.

We first consider to construct the LCU formulas for $V_{K}(x)$ from Taylor-series expansion Berry et al. (2015). In the $0$th-order case, when no Trotter formula is introduced, the Trotter remainder $V_{0}(x)$ is the short-time evolution $U(x)$ itself, which can be expanded as

Here, $\mathrm{Poi}(s;a):=e^{-a}a^{s}/s!$ is the Poisson distribution. $F_{0,s}(x)$ denotes the $s$-order expansion term. Eq. 4 illustrated in Fig. 3(a) is a LCU formula of $V_{0}(x)$ with $1$-norm $\mu_{x}=e^{\lambda x}$. The $1$-norm of the overall evolution $U(t)=U(x)^{\nu}$ is then $\mu=e^{\lambda t}$, which, unfortunately, grows exponentially with respect to $t$ regardless of how much we increase the segment number $\nu$. This implies that the direct random-sampling implementation of the LCU formula $V_{0}(x)$ in Eq. 4 is not feasible. In Ref. Berry et al. (2015), the authors discuss the coherent implementation of $V_{0}(x)$ instead and find that one can achieve good time and accuracy dependence in that scenario.

When focusing on the random-sampling implementation, we need to suppress the $1$-norm $\mu_{x}$ of each segment. To this end, we first consider the usage of the Trotter formula. For example, if we apply first-order Trotter formula $S_{1}(x)=\prod_{l=1}^{L}e^{-ixH_{l}}$ in each segment, the first-order remainder $V_{1}(x):=U(x)S_{1}(x)^{\dagger}$ can be expanded as

where $\{r_{1},...,r_{L}\}$ denotes $L$ expansion variables related to the Trotter formula and $s:=r+\sum_{l=1}^{L}r_{l}$. $F_{1,s}$ denotes the $s$-order expansion term of $V_{1}(x)$. The $1$-norm of $V_{1}(x)$ in Eq. 5 is $e^{2\lambda x}$, which seems to be even larger than the $0$th-order case. However, since $V_{1}(x)$ denotes the (multiplicative) Trotter error, we have $F_{1,1}(x)=0$. Using this condition, we can rewrite $V_{1}(x)$ as

as illustrated in Fig. 3. From the Taylor-series expansion, we can bound the $1$-norm of the new LCU formula in Eq. 6 by $\mu_{x}=e^{2\lambda x}-(2\lambda x)\leq e^{(2\lambda x)^{2}}$. In this way, we reduce $\mu_{x}$ from $1+\mathcal{O}(\lambda x)$ to $1+\mathcal{O}((\lambda x)^{2})$. The $1$-norm of the overall time evolution becomes $\mu=\mu_{x}^{\nu}\leq\exp((2\lambda t)^{2}/\nu)$. As a result, by increasing the segment number $\nu$ — or equivalently reducing the unit evolution time $x$ — we can decrease the 1-norm $\mu$, leading to a lower sampling cost. If we choose the segment number as $\nu=\mathcal{O}((\lambda t)^{2})$, $\mu$ will remain constant.

From the above discussion, it is clear that, to reduce the 1-norm $\mu$ of the overall LCU formula for $U=e^{-iHt}$, our main objective is to suppress the leading-order term of the $1$-norm remainder $\mu_{x}-1$ for each segment, which determines the number of segments $\nu$ and hence circuit depth when $\mu$ is set to be a constant.

We can further reduce the $1$-norm $\mu_{x}$ of $V_{1}(x)$ by taking advantage of the structure of Trotter errors. For an anti-Hermian Pauli operator $\pm iP$ where $P\in\{I,X,Y,Z\}^{\otimes n}$, we have the following Euler’s formula,

Here, $\theta=\tan^{-1}(y)$ and we suppose $0<y<1$. The 1-norm of the left-hand side of Eq. 7 is $1+y$, while the 1-norm of the right-hand side is $\sqrt{1+y^{2}}<1+\frac{y^{2}}{2}=1+\mathcal{O}(y^{2})$. As a result, the exponent of $\mu_{x}-1$ is effectively doubled. In subsection IV.2, we prove that the expansion terms $F_{1,2}$ and $F_{1,3}$ in the LCU formula Eq. 5 are anti-Hermitian. As a result, we can further suppress $\mu_{x}$ by pairing the terms in $F_{1,s}(x)$ with $F_{1,0}=I$ using Euler’s formula in Eq. 7. When $y=x^{s}$, the paired formula $R_{1,s}(x)$ as a summation of Pauli rotation unitaries owns the $1$-norm of $\mu_{x}=1+\mathcal{O}(x^{2s})$, whose $x$ dependence is doubled, as illustrated inFig. 3(c).

To generalize the discussion, for the $K$th-order Trotter error $V_{K}(x)$ ($K=2k$, $k\in\mathbb{N}_{+}$), we have $F_{K,1}=...F_{K,K}=0$ Suzuki (1990); Childs et al. (2021). Moreover, we can show that $F_{K,K+1}$, $F_{K,K+2}$,…, $F_{K,2K+1}$ are anti-Hermitian. We call the term with $s=K+1$ to $2K+1$ the leading-order terms. The algorithm utilizing $K$th-order Trotter formula and the paired idea in the LCU construction is called the $K$th-order PTSC algorithm. We provide the detailed algorithm description and gate-complexity analysis in Sec. IV.

The PTSC algorithm is generic for the $L$-sparse Hamiltonian $H=\sum_{l=1}^{L}H_{l}=\lambda\sum_{l=1}^{L}p_{l}P_{l}$ with $\lambda=\sum_{l}\|H_{l}\|$ and $\{P_{l}\}$ are Pauli matrices. It can be implemented using a simple and universal classical random-sampling procedure: first, we sample the order $s$ from the Taylor-series expansion, and then we sample a Pauli string based on the Hamiltonian coefficients. With the random-sampling implementation, we prove that by appending few gates after the $K$th-order Trotter formula with only one ancillary qubit, one can improve the time scaling from $1+1/K$ to $1+1/(2K+1)$ and exponentially improves the accuracy scaling of the $K$th-order Trotter formula compared to the $K$th-order Trotter. We have the following theorem.

Theorem 1 (Informal, see Theorem 1 in Sec. IV) In a $K$th-order paired Taylor-series compensation algorithm ($K=0,1$ or $2k$, $k\in\mathbb{N}_{+}$), the gate complexity in a single round is $\mathcal{O}\left((\lambda t)^{1+\frac{1}{2K+1}}(\kappa_{K}L+\frac{\log(1/\varepsilon)}{\log\log(1/\varepsilon)})\right)$, where $\lambda=\sum_{l}\|H_{l}\|$, $\kappa_{K}=K$ when $K=0$ or $1$, $\kappa_{K}=2\times 5^{K/2-1}$ otherwise.

From Theorem 1, we observe that setting $K=0$, i.e., not using Trotter formulas, still yields a valid PTSC algorithm by pairing $F_{0,1}$ with $I$. In this case, the gate complexity is independent of the sparsity $L$ but quadratically dependent on $t$, similar to the algorithm in Ref. Wan et al. (2022). Conversely, when using a $K$th-order Trotter formula, the PTSC algorithm becomes $L$ dependent with an almost linear dependence on $t$. In both cases, the PTSC algorithms achieve high simulation accuracy $\varepsilon$. We expect the $0$th-order algorithm to be particularly useful for quantum chemistry Hamiltonians with large $L$, while higher-order algorithms are better suited for generic $L$-sparse Hamiltonians with long simulation times $t$.

The PTSC algorithms above are generic and applicable to any Hamiltonian. When we consider the detailed structure of Hamiltonians, we could make the compensation algorithms more efficient by taking advantage of the commutation relationship of the terms in the Hamiltonians, which was formerly also studied in the Trotter algorithms Childs et al. (2021).

We will take the first-order Trotter remainder $V_{1}(x)$ as an illustrative example. Following the Taylor-series expansion in Eq. 5, the second-order term $F_{1,2}(x)$ in $V_{1}(x)$ can be written as

Since $F_{1,2}(x)$ is anti-Hermitian, all the Hermitian expansion terms in Eq. 8 will cancel out. We can then simplify $F_{1,2}(x)$ as,

which is a summation of $L(L-1)$ commutators. Since the commutators of Hermitian operators are always anti-Hermitian, this implies that the nested-commutator expansion of $F_{1,2}(x)$ in Eq. 9 is compact enough since there is no Hermitian expansion terms in it.

For a common physical Hamiltonian with locality constraints, we can take advantage of the commutator-form expression like Eq. 9. For example, for an $n$-qubit lattice Hamiltonian with the form,

where the summand $H_{j,j+1}$ acts on the $j$th and $(j+1)$th vertices. We can split the Hamiltonian to two components $H=A+B$ where $A:=\sum_{j:\mathrm{even}}H_{j,j+1},\quad B:=\sum_{j:\mathrm{odd}}H_{j,j+1}$ , so that the summands $H_{j,j+1}$ commute with each other in each component. We denote the norm of each Hamiltonian summand as

Now, suppose we estimate the $1$-norm of $F_{1,2}(x)$ of the lattice Hamiltonian based on Eq. 9, we can see that there are only $n$ nonzero terms: for any given Hamiltonian component $H_{j,j+1}$, only $H_{j-i,j}$ and $H_{j+1,j+2}$ do not commute with it. Then, the norm of $F_{1,2}(x)$ is bounded by

Comparing with the original bound in Eq. 62, $\|F_{1,2}(x)\|_{1}\leq\eta_{2}=\mathcal{O}((\lambda x)^{2})=\mathcal{O}(n^{2}x^{2})$, we improve the system-size-related factor $n$. The improved system-size $n$ dependence of the 1-norm $\|F_{1,2}(x)\|_{1}$ suggests a corresponding improvement in the system-size dependence of the gate complexity for the nested commutator algorithm.

Now, we are going to generalize the idea above. In Sec. V, we show how to expand the second- and third-order terms of $V_{1}(x)$ as a summation of nested commutators,

where $\mathcal{A}_{\mathrm{com}}^{(s)}(A_{L},...,A_{1};B)$ is defined to be

We also use the adjoint notation $\mathrm{ad}_{A_{L}}...\mathrm{ad}_{A_{1}}B:=[A_{L},...[A_{1},B]...]$. It is easy to check that the form of $F_{1,s}(x)$ ($s=2$ or $3$) in Eq. 13 is anti-Hermitian, which is consistent with the discussion in the paired Taylor-series algorithm in the previous section. We can also generalize the method to the case of $K$th-order Trotter remainder, that is, to express the expansion terms $F_{K,s}$ of the $K$th-order Trotter remainder based on the nested commutators.

Based on the nested-commutator expansion, we propose the $K$th-order nested-commutator compensation (NCC) algorithm. As is illustrated in Fig. 4, in the construction of NCC, we first utilize the nested-commutator forms of Trotter error terms from $K+1$ to $2K+1$ order, i.e., the leading-order terms; then we apply the order-pairing techniques similar to PTSC in Fig. 3(c) to further suppress the $1$-norm of $\mu_{x}$. A key difference between NCC and PTSC algorithms is that in PTSC algorithms, we compensate the Trotter error $V_{K}(x)$ up to arbitrary order; while in NCC algorithms, we only compensate $V_{K}(x)$ for leading-order terms, which shrinks the error from $\mathcal{\mathcal{missing}}{O}(x^{K+1})$ to $\mathcal{O}(x^{2K+2})$ in one slice with the sampling cost $\mu=1+\mathcal{O}(\frac{\|C_{K}\|_{1}}{(K+1)!}x^{2K+2})$. The gate complexity estimation is then converted to the calculation of the $1$-norm of the commutator $\|C_{K}\|_{1}$. For instance, if we consider the $n$-qubit lattice Hamiltonian models in Eq. 10, then we can prove that $\|C_{K}\|_{1}=\mathcal{O}(n)$. We can then provide the following performance guarantee for the NCC algorithms.

Theorem 2 (Informal, see Theorem 2 in Sec. V) In a $K$th-order nested-commutator compensation (NCC) algorithm ($K=1$ or $2k$) with $n$-qubit lattice Hamiltonians, the gate complexity in a single round is $\mathcal{O}(n^{1+\frac{2}{2K+1}}t^{1+\frac{1}{2K+1}}\varepsilon^{-\frac{1}{2K+1}})$.

Compared to the performance of $K$th-order Trotter algorithm $\mathcal{O}((nt)^{1+1/K}\varepsilon^{-1/K})$ Childs and Su (2019), we achieve $t$- and $\varepsilon$-dependence better than $2K$th-order Trotter using only $K$th-order Trotter formula with simple compensation gates of Pauli-rotation operators. To generalize the result in Theorem 2, we also study the performance of Nested Commutator (NCC) algorithms when applied to a general Hamiltonian in subsection V.3.

A simple implementation of the Trotter-LCU algorithm in Fig. 2(b) or Fig. 2(c) requires not only an easy-to-implement quantum circuit but also efficient classical random sampling of Pauli operators from the Trotter remainder $V_{K}(x)$. We now briefly discuss how to realize an efficient classical random sampling in PTSC and NCC algorithms, that is, with a space resource of $\mathcal{O}(\kappa_{K}K)$ and time resource of $\mathcal{O}(K(\log{L}+\log\kappa_{K}))$ where $L$ is the sparsity of the Hamiltonian.

A key idea for achieving efficient sampling is to use a multistage hierarchical sampling algorithm. Rather than fully expanding the Trotter remainder into a direct summation of Pauli operators, we structure the LCU formula into multiple layers. This allowss us to decompose the overall Pauli operator sampling process into a series of simpler, more manageable sampling tasks.

In the PTSC algorithm, the Trotter remainder in Eq. 5 is derived by expanding the time-evolution $e^{-ixH}$ and each Hamiltonian summand term $e^{ixH_{l}}$ by Taylor series independently. As a result, the sampling can be done by first sampling the overall expansion order $s$, then sample the individual expansion order $r$ of Hamiltonian $H$ or the expansion order $\vec{r}:=[r_{1},r_{2},...,r_{L}]$ of the summands $\{H_{l}\}$. The sampling of $r$ and $\vec{r}$, following the analysis in Sec. IV, can be done based on a multinomial distribution $\text{Mul}(r,\vec{r};\{\frac{1}{2},\frac{\vec{p}}{2}\};s)$ where $\vec{p}:=[p_{1},p_{2},...,p_{L}]$ denotes the normalized coefficient factor of $H$ defined in Eq. 3. For the sampled Hamiltonian $H$, we further sample the summands $H_{l}$ inside based on $\vec{p}$. We summarize the sampling algorithm in Fig. 5.

In the NCC algorithm, the Trotter remainder is expanded based on the adjoint operators. For example, for the lattice Hamiltonian $H=A+B$ in Eq. 10, we can write the second- and third-order Trotter remainder of $V_{1}(x)$ as,

We can first sample the expansion order $s=2$ or $3$. If $s=3$, we then sample the specific commutator, i.e., $\mathrm{ad}_{B}\mathrm{ad}_{A}B$ or $\mathrm{ad}_{A}^{2}B$. For the given commutator, for example, $\mathrm{ad}_{B}\mathrm{ad}_{A}B=[B,[A,B]]$, we first randomly sample a summand $H_{j,j+1}$ for the rightmost $B$ as the starting point of the adjoint operator. The action of the subsequent $\mathrm{ad}_{A}$ and $\mathrm{ad}_{B}$ will enlarge the support of $H_{j,j+1}$, but within a “light-cone” region shown in Fig. 6. We then sample the Hamiltonian summand $H_{j_{1},j_{1}+1}$ and $H_{j_{2},j_{2}+1}$ for the adjoint operators $\mathrm{ad}_{B}$ and $\mathrm{ad}_{A}$, but within the light-cone region. This will ensure our sampling to be efficient and with nested-commutator scaling.

A problem of sampling the Hamiltonian summands $H_{j,j+1}$ in the commutator is that, since the Hamiltonian $H$ may be nonhomogeneous, the $1$-norm of $H_{j,j+1}$ with different $j$ may be different. This will complicate the sampling algorithm, since we need to first calculate the $1$-norm of all the elementary commutators with the form of

for all $j,j_{1}$ and $j_{2}$. When the Trotter order $K$ and the expansion order $s=K+1,...,2K+1$ gets larger, the number of elementary adjoint operators will increase exponentially. Consequently, we cannot estimate the $1$-norm of all the elementary commutator.

To solve this problem but still keep the advantage of the NCC algorithm, we introduce the following Hamiltonian “padding” technique to ensure all the elementary nested commutator own the same $1$-norm. Consider a Hamiltonian summand $H_{j,j+1}$, which can be expanded to some Pauli operators,

where $\alpha_{j}^{(\omega)}$ is a positive number and $P_{j,j+1}^{(\omega)}$ is a normalized Pauli operator whose support is on qubit $j$ and $j+1$. Recall that $\Lambda_{1}:=\max_{j}\|H_{j,j+1}\|_{1}$ and $\|H_{j,j+1}\|_{1}:=\sum_{\omega}\alpha_{j}^{(\omega)}$. When $\|H_{j,j+1}\|_{1}$ is smaller than $\Lambda_{1}$, we add extra trivial terms $\pm I$ in the Pauli decomposition of $H_{j,j+1}$,

where $\delta\Lambda_{1}:=\Lambda_{1}-\|H_{j,j+1}\|_{1}$. Eq. 18 holds naturally, but now with a manually predetermined $1$-norm value $\Lambda_{1}$. Similarly, we can pad all the elementary commutators with the form

so that their $1$-norms are all $2\Lambda_{1}^{2}$. In this way, we ignore the commutator relationship between $H_{j,j+1}$ and $H_{j_{1},j_{1}+1}$ as long as they are in the light-cone region.

After the padding procedure described above, all elementary nested commutators of the same order will have the same $1$-norm. This property allows us to uniformly sample these commutators: the starting summand $H_{j,j+1}$ is sampled uniformly from those in $B$, and the subsequent $H_{j_{1},j_{1}+1}$ and $H_{j_{2},j_{2}+1}$ are sampled uniformly within the light-cone region. However, when the starting summand $H_{j,j+1}$ is near the boundary, applying a few adjoint operators may cause it to touch the boundary. This reduces the number of possible elementary nested commutators compared to those starting from the center, resulting in different $1$-norms for $\text{ad}_{A}H_{j,j+1}$ depending on $j$. This complicates the sampling of the starting summand $H_{j,j+1}$. To resolve this issue and ensure uniform sampling of $H_{j,j+1}$, we introduce virtual qubits at the boundary, as illustrated in Fig. 6, and pad the virtual qubits with $0$-summed $\pm I$ terms. Since we only perform $\pm I$ operation on the virtual qubits, we do not need to introduce it in the real experiments.

We remark that, our padding method preserve the locality structure. As a result, the performance guarantee in Theorem 2 still holds. We summarize the sampling algorithm in Fig. 7. If we consider the Heisenberg Hamiltonian

as an example, where $\vec{\sigma}_{i}:=(X_{i},Y_{i},Z_{i})$ is the vector of Pauli operators on the $i$th qubit, we can define the summand $H_{j,j+1}$ to be

In this case, $\Lambda_{1}=\|H_{j,j+1}\|=4$. The probability distribution $\bar{p}_{j}^{(\omega)}$ in Fig. 7 is to uniformly sample the $XX$, $YY$, $ZZ$ and $ZI$ term. As a demonstration, we explicitly present the algorithm to sample the Pauli-rotation operator in the first-order NCC algortihm for the Heisenberg Hamiltonian in Eq. 20 in Algorithm 1.

As a final remark, the sampling procedure in both PTSC and NCC algorithms are independent of the implementation of the quantum circuit and the measurement outcome. Thanks to this property, we can perform the classical sampling during the quantum circuit implementation or even generate the sampled Pauli matrices before the implementation of the quantum circuits.

In Table 1, we compare the implementation complexity and the gate complexity in a single round of experiment of the $0$th-order PTSC, $K$th-order PTSC, and $K$th-order NCC algorithms to previous Hamiltonian simulation algorithms. For a fair comparison, we set the 1-norm of all LCU formulas $\mu$ to be constant. In this case, the sample complexity of PTSC and NCC algorithms incurs a $\mu^{4}$ overhead compared to standard sampling from the $K$th-order Trotter or post-Trotter algorithms.

We show that by inserting a few randomly sampled Pauli-rotation gates after each Trotter segment, as illustrated in Fig. 1, both PTSC and NCC achieve improved accuracy and time dependence. The gate counts for PTSC exhibit logarithmic dependence on accuracy, $\log(1/\varepsilon)$, while the NCC gate counts show improved system-size dependence.

To demonstrate how the Trotter-LCU algorithms can help reduce gate costs in practical scenarios, we estimate the single-shot gate count of random-sampling Trotter-LCU algorithms and compare it with the state-of-the-art Trotter algorithm, i.e., fourth-order Trotter Suzuki (1991); Childs et al. (2019, 2021). To ensure a fair comparison, we set $1$-norm of the LCU formula to be $\mu=2$. Based on Proposition 1, this implies that the random-sampling Trotter-LCU algorithm will require an additional factor of $16$ in the sample number to estimate the properties of the target state $U\rho U^{\dagger}$ where $U:=e^{iHt}$ to a given precision compared to the normal Trotter or coherent implementation of LCU algorithms.

We compile their quantum circuits to $\mathcal{CNOT}$ gates, single-qubit Clifford gates, and single-qubit $Z$-axis rotation gates $R_{z}(\theta)=e^{i\theta Z}$. Here, we mainly compare the number of $R_{z}(\theta)$ gates since they are the most resource-consuming part on a fault-tolerant quantum computer Litinski (2019). The $\mathcal{CNOT}$ gate-number comparison results can be found in Appendix G, which is similar to the $R_{z}(\theta)$ gate comparison. We also compare the gate counts of the Trotter-LCU algorithms to the coherent-implementation of LCU Berry et al. (2014, 2015) and QSP Low and Chuang (2017, 2019) in Appendix G.

In the first comparison, we consider the simulation of generic Hamiltonians without the usage of commutator information, for which we choose the $2$-local Hamiltonian, $H=\sum_{i,j}X_{i}X_{j}+\sum_{i}Z_{i}$ where $X_{i}$ and $Z_{i}$ are the Pauli matrices on the $i$th qubit. Fig. 8(a,b) show the gate counts for the fourth-order Trotter formula and the PTSC Trotter-LCU algorithms with different orders with an increasing time $t$ and increasing system size $n$, respectively. The gate counting method for fourth-order Trotter with random permutation is based on the analytical bounds in Ref. Childs et al. (2019).

From Fig. 8(a,b), we can clearly see the advantage of composition of Trotter and LCU formulas: if we do not use LCU and merely apply Trotter formulas, the gate resource of fourth-order Trotter suffers from a large overhead that is 2 orders of magnitudes larger than the PTSC algorithms. Moreover, if we increase $\varepsilon$ from $10^{-3}$ to $10^{-5}$, we can see a clear increase of the gate resources for the Trotter algorithm. For the PTSC algorithms, however, the gate number is almost not affected by $\varepsilon$ since they enjoy a logarithmic $\varepsilon$-dependence.

On the other hand, if we do not use Trotter formula, the $0$th-order PTSC algorithm owns a quadratically worse $t$-dependence ($\mathcal{O}(t^{2})$) than the fourth-order Trotter ($\mathcal{O}(t^{1.25})$), second-order PTSC ($\mathcal{O}(t^{1.2})$) and fourth-order PTSC ($\mathcal{O}(t^{1.11})$) algorithms in Fig. 8(a). For a short-time evolution $t=n$, the system-size dependence of second- or fourth-order PTSC in Fig. 8(b) outperforms $0$th-order PTSC algorithm. In the case when long-time Hamiltonian simulation is required, for example, $t$ should be set to $10^{3}$ for the phase estimation Campbell (2019), the advantage of $2$nd or fourth-order PTSC to $0$th-order PTSC will be more obvious. The composition of Trotter and LCU formulas enables $2$nd or fourth-order PTSC to enjoy good $t$ and $\varepsilon$ dependence and small gate-resource overhead simultaneously.

Next, we compare the gate count when simulating the lattice models, where the commutator analysis will help remarkably reduce the gate count. We consider the Heisenberg Hamiltonian $H=\sum_{i}\vec{\sigma}_{i}\vec{\sigma}_{i+1}+\sum_{i}Z_{i}$ using the nested-commutator bounds. In Fig. 8(c), we choose $n=t=12$ and $50$ and show the gate number with respect to the accuracy requirement $\varepsilon$. The fourth-order Trotter error analysis is based on the nested-commutator bound (Proposition M.1 in Ref. Childs et al. (2021)), which is currently the tightest Trotter error analysis. The performance of our second-order NCC algorithm is based on the analytical bound in Sec. H in the Appendix. We mainly present results for the second-order NCC algorithm due to its simplicity and leave precise higher-order NCC gate count analysis for future study.

From Fig. 8(c), we can see that, while enjoying near-optimal system-size scaling similar to the fourth-order Trotter algorithm which is currently the best one for lattice Hamiltonians Childs et al. (2021); Childs and Su (2019), the second-order NCC algorithm shows better accuracy dependence than fourth-order Trotter algorithm. Particularly, using the same gate number as the fourth-order Trotter, we are able to achieve a 3 to 4 orders of magnitudes higher accuracy $\varepsilon$.

The most natural way to approximate $U(x)$ is to apply the Lie-Trotter-Suzuki formulas Suzuki (1990, 1991). Hereafter, we refer to them as Trotter formulas. The first-order Trotter formula is

Here, $\vec{H}:=[H_{1},H_{2},...,H_{L}]$. In Eq. 23, we simplify the notation of the sequential products from $l=1$ to $L$ with the same form, using the arrow to denote the ascending direction of the dummy index $l$. The Hermitian conjugation of $S_{1}(x)$ can similarly be written as $S_{1}(x)^{\dagger}=\overset{\rightarrow}{\prod}e^{ix\vec{H}}$.

The second-order Trotter formula is