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

Pei Zeng<EMAIL_ADDRESS>Pritzker School of Molecular Engineering, The
University of Chicago, Illinois 60637, USA  Jinzhao Sun
<EMAIL_ADDRESS>QOLS, Blackett Laboratory, Imperial College London,
London SW7 2AZ, United Kingdom  Quantum Advantage Research, Beijing 100080,
China  Liang Jiang<EMAIL_ADDRESS>Pritzker School of Molecular
Engineering, The University of Chicago, Illinois 60637, USA  Qi Zhao
<EMAIL_ADDRESS>QICI Quantum Information and Computation Initiative,
Department of Computer Science, The University of Hong Kong, Pokfulam Road,
Hong Kong

###### Abstract

Trotter and linear-combination-of-unitary (LCU) are two popular Hamiltonian
simulation methods. We propose Hamiltonian simulation algorithms using LCU to
compensate Trotter error, which enjoy both of their advantages. By adding few
gates after the Kth-order Trotter, we realize a better time scaling than 2Kth-
order Trotter. Our first algorithm exponentially improves the accuracy scaling
of the Kth-order Trotter formula. In the second algorithm, we consider the
detailed structure of Hamiltonians and construct LCU for Trotter errors with
commutator scaling. Consequently, for lattice Hamiltonians, the algorithm
enjoys almost linear system-size dependence and quadratically improves the
accuracy of the Kth-order Trotter.

Introduction. Hamiltonian simulation, i.e., to simulate the real-time
evolution $U(t)=e^{-iHt}$ of a physical Hamiltonian $H=\sum_{l}H_{l}$, is
considered to be a natural and powerful application of quantum computing
Feynman et al. (1982). It can also be used as an important subroutine in many
other quantum algorithms like ground-state preparation Abrams and Lloyd
(1999); Aspuru-Guzik et al. (2005), optimization problems Farhi et al. (2014);
Zhou et al. (2020), quantum linear solvers Harrow et al. (2009). To pursue
real-world applications of Hamiltonian simulation with near-term quantum
devices, we need to design feasible algorithms with small space complexity
(i.e., qubit number) and time complexity (i.e., circuit depth and gate
number).

One of the most natural Hamiltonian simulation methods is based on Trotter
formulas Lloyd (1996); Suzuki (1990, 1991); Berry et al. (2007); Campbell
(2019); Childs et al. (2019); Childs and Su (2019); Childs et al. (2021),
which approximate the real-time evolution of $H=\sum_{l=1}^{L}H_{l}$ by the
product of the simple evolution of its summands $e^{-iH_{l}t}$. Besides its
prominent advantage of simple realization without ancillas, Trotter methods
are recently rigorously shown to enjoy commutator scaling Childs and Su
(2019); Childs et al. (2021), i.e., the Trotter error is only related to the
nested commutators of the Hamiltonian summands $\\{H_{l}\\}$. This is very
helpful for the Hamiltonians with strong locality constraints. Hereafter, we
refer to the error analysis utilizing the detailed structure of Hamiltonians
as commutator scaling analysis. For example, when we consider $n$-qubit
lattice Hamiltonians, the gate cost of high-order Trotter methods is almost
linear to the system size $n$, which is nearly optimal Childs and Su (2019).
The major drawback of the Trotter methods is its polynomial gate cost to the
inversed accuracy $1/\varepsilon$, $\rm Poly(1/\varepsilon)$. In many
applications where high-precision simulation is demanded to obtain practical
advantages over the existing classical algorithms Reiher et al. (2017), the
gate cost of Trotter methods is large.

In recent years, we have seen the developments of “post-Trotter” algorithms
with exponentially improved accuracy dependence Berry et al. (2014, 2015); Low
and Chuang (2019, 2017). Due to the smart choice of the expansion formulas
(i.e., Taylor series Berry et al. (2015) or Jacobi-Anger expansion Low and
Chuang (2017)), these post-Trotter methods are able to catch the dominant
terms in the time evolution operator $U(t)$ with polynomially increasing gate
resources, leading to a logarithmic gate-number dependence on the accuracy
requirement $1/\varepsilon$. Unlike Trotter methods, these advanced algorithms
are not able to utilize the specific structure of Hamiltonians due to the lack
of commutator-based error form. Consequently, for instance, for $n$-qubit
lattice Hamiltonians, their gate complexities are $\mathcal{O}(n^{2})$, which
is worse than those in Trotter algorithms $\mathcal{O}(n^{1+o(1)})$.
Furthermore, these post-Trotter algorithms require the implementation of
linear-combination-of-unitary (LCU) formulas Childs and Wiebe (2012) or block
encoding of Hamiltonians Low and Chuang (2019) which often costs many
ancillary qubits and multi-controlled Toffoli gates. This is still
experimentally challenging in a near-term or early fault-tolerant quantum
computer Lin and Tong (2021). To reduce the hardware requirement of compiling
the LCU formulas, recent studies focus on a random-sampling implementation of
a LCU formula Childs and Wiebe (2012); Yang et al. (2021); Wan et al. (2022);
Faehrmann et al. (2021), where the elementary unitaries are sampled to realize
the LCU formula statistically. Similar ideas have also been studied in the
ground-state algorithms Lin and Tong (2020); Zeng et al. (2021); Zhang et al.
(2022).

Here, we develop the theory of composite algorithms that combine the inherent
advantages of Trotter and LCU methods — easy implementation, high precision,
and commutator scaling — by performing the Trotter method and then
compensating the Trotter error with the LCU methods. We provide versatile
realizations of this composite algorithm based on either random-sampling Yang
et al. (2021); Wan et al. (2022); Faehrmann et al. (2021) or coherent
implementation Berry et al. (2014, 2015) of LCU formulas. Based on the methods
to construct the LCU formulas, we propose two types of composite algorithms
with unique advantages in different scenarios. Our first algorithm is generic
for the $L$-sparse Hamiltonian $H=\sum_{l=1}^{L}H_{l}$ with
$\lambda=\sum_{l}\|H_{l}\|$. 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 achieve a gate cost of $\mathcal{O}((\lambda
t)^{1+1/(2K+1)}(L+\log(1/\varepsilon)))$, which improves the time scaling from
$1+1/K$ to $1+1/(2K+1)$ and exponentially improves the accuracy scaling of the
Kth-order Trotter formula compared to the $K$th-order Trotter. In the second
algorithm, we consider the locality structure of Hamiltonians and take
advantage of the nest-commutator forms of Trotter errors Childs and Su (2019);
Childs et al. (2021) in the LCU formula construction. As a result, our
algorithm enjoys almost linear system-size dependence, and achieves a better
time and accuracy dependence than 2Kth-order Trotter only with Kth-order
Trotter’s gate complexity.

Figure 1: (a) In the $K$th-order Trotter-LCU algorithm, we first implement
$K$th-order Trotter formula, then compensate the remainder using LCU methods.
(b) Random-sampling implementation of the LCU formula. We sample the
elementary unitaries $V_{i}$ and $V_{j}$ independently based on the LCU
formula of $V_{K}(x)$ in Eq. 1 and implement the controlled $V_{i}$ and
$V_{j}$ gate. Then the Trotter remainder $V_{K}(x)$ can be compensated in a
Hadamard-test type circuit.

Trotter-LCU algorithms. The algorithm is based on a series connection of
Trotter and LCU algorithms, shown in Fig. 1. We call the composite algorithm
with the $K$th-order Trotter and LCU formula the $K$th-order Trotter-LCU
method. For consistency, we call the algorithm with no Trotter formula the
$0$th-order Trotter-LCU method. In our algorithms, we decompose the time
evolution $U(t)=e^{-iHt}$ to $\nu$ segments, each with a small evolution time
$x=t/\nu$. In each segment, we first perform the $K$th-order Trotter formula
$S_{K}(x)$ ($K=0,1,2k$, $k\in\mathbb{N}_{+}$). For consistency, we denote the
$0$th-order Trotter formula as $S_{0}(x)=I$. Then, we compensate the
multiplicative Trotter remainder $V_{K}(x):=U(x)S_{K}(x)^{\dagger}$ by
implementing the LCU formula of $V_{K}(x)$. A $(\mu,\varepsilon)$-LCU formula
of an operator $V$ is defined by a set of unitaries $\\{V_{i}\\}$ and related
probability distribution $\vec{p}(V)=\\{p_{i}\\}_{i}$ with

$\tilde{V}=\mu\sum_{i=0}^{\Gamma-1}p_{i}V_{i},$ (1)

such that the spectral norm distance $\|V-\tilde{V}\|\leq\varepsilon$. We
denote $\mu:=\|\tilde{V}\|_{1}$ as the $1$-norm of this LCU formula.

We consider two ways to implement the LCU formula. In the random-sampling
implementation Yang et al. (2021); Wan et al. (2022); Faehrmann et al. (2021),
we sample the elementary unitaries $\\{V_{i}\\}$ based on the corresponding
probability $\\{p_{i}\\}$ and embed the realization of the LCU formula into a
Hadamard test Kitaev (1995). In the coherent implementation Berry et al.
(2014, 2015), we introduce $m=\lceil\log_{2}\Gamma\rceil$ ancillary qubits to
realize the superposition of different $\\{V_{i}\\}$. Our major focus is on
the random-sampling implementation shown in Fig. 1(b), where we can estimate
the properties of the target state $\rho=U(t)\rho_{0}U(t)^{\dagger}$ with only
one ancillary qubit. Given a $(\mu,\varepsilon)$-LCU formula in Eq. 1, we can
estimate an arbitrary observable value $O$ on $\rho$ based on the formula,

$\mathrm{Tr}(O\rho)\approx\mu^{2}\sum_{i,j}p_{i}p_{j}\mathrm{Tr}(OV_{i}\rho_{0}V_{j}^{\dagger}).$
(2)

While 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{O}(\mu^{4}/\varepsilon^{2})$ sampling resource, which owns an extra
$\mu^{4}$ overhead compared to the normal Hamiltonian simulation algorithms
Faehrmann et al. (2021). To make the algorithm efficient, we need to set $\mu$
to be a constant.

To construct the LCU formulas, we can directly apply the Taylor-series
expansion Berry et al. (2015) on $V_{K}(x)$. We start our discussion from
$0$th-order case, where no Trotter formula is performed. The $0$th-order
remainder $V_{0}(x):=U(x)$ can be expanded as

$\displaystyle V_{0}(x)$
$\displaystyle=e^{-ixH}=\sum_{s=0}^{\infty}F_{0,s}(x)=\sum_{s=0}^{\infty}\frac{(-ix)^{s}}{s!}H^{s},$
(3) $\displaystyle=e^{\lambda x}\sum_{s=0}^{\infty}\mathrm{Poi}(s;\lambda
x)\sum_{l_{1},...,l_{s}}p_{l_{1}}...p_{l_{s}}P_{l_{1}}...P_{l_{s}}.$

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. 3 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$.

To reduce the $1$-norm $\mu_{x}$ in each segment, we consider introducing the
Trotter formula. For example, if we apply $1$st-order Trotter formula
$S_{1}(x)=\prod_{l=1}^{L}e^{-ixH_{l}}$ in each segment, the $1$st-order
remainder $V_{1}(x):=U(x)S_{1}(x)^{\dagger}$ can be expanded as

$\displaystyle V_{1}(x)$
$\displaystyle=e^{-ixH}\prod_{l=L}^{1}e^{ixH_{l}}=\sum_{s=0}^{\infty}F_{1,s}(x)$
(4)
$\displaystyle=\sum_{r;r_{1},r_{2},...,r_{L}}\frac{(-ix)^{r}}{r!}H^{r}\prod_{l=L}^{1}\frac{(ix)^{r_{l}}}{r_{l}!}H^{r_{l}},$

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. 4 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 $F_{1,1}$ as

$V_{1}(x)=I+\sum_{s=2}^{\infty}F_{1,s}(x).$ (5)

From the Taylor-series expansion, we can bound the $1$-norm of the new LCU
formula in Eq. 5 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)$. If we choose the
segment number $\nu=\mathcal{O}((\lambda t)^{2})$, $\mu$ will 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. In Appendix G, we show how to expand the
$2$nd- and $3$rd-order terms of $V_{1}(x)$ as a summation of nested
commutators,

$\displaystyle F_{1,s}(x)$
$\displaystyle=\frac{x^{s}}{s!}C_{s-1}=\frac{(-ix)^{s}}{s!}\Big{(}\mathcal{A}_{com}^{(s-1)}(H_{L},...,H_{1};H)$
(6)
$\displaystyle\,-\sum_{l=1}^{L}\mathcal{A}_{com}^{(s-1)}(H_{L},...,H_{l+1};H_{l})\Big{)},\quad
s=2,3,$

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

$\displaystyle\sum_{m_{1}+...+m_{L}=s}\binom{s}{m_{1},...,m_{L}}\mathrm{ad}_{A_{L}}^{m_{L}}...\mathrm{ad}_{A_{1}}^{m_{1}}B.$
(7)

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 $F_{1,s}(x)$ ($s=2$ or $3$) in Eq. 6 is anti-
Hermitian which can be written as,

$F_{1,s}(x)=\frac{i\|C_{s-1}\|_{1}}{s!}x^{s}\sum_{\gamma}\Pr(\gamma)P_{\gamma},$
(8)

where $\\{P_{\gamma}\\}$ are Pauli operators with real coefficient $\pm 1$. 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
$I+iyP=\sqrt{1+y^{2}}e^{i\theta P}$ with $\theta=\tan^{-1}(y)$. When
$y=x^{s}$, the paired formula as a summation of Pauli rotation unitaries owns
the $1$-norm of $\mu_{x}=1+\mathcal{O}(x^{2s})$, whose $x$-dependence is
doubled.

Similarly, 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 expand $F_{K,K+1}$, $F_{K,K+2}$,…, $F_{K,2K+1}$
by the summation of nested commutators which 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 paired Taylor-series compensation (PTSc) algorithm. We
provide the detailed algorithm description and gate-complexity analysis in
Appendix F.

Theorem 1 (Informal, see Theorem 1 in Appendix F) 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\cdot 5^{K/2-1}$ otherwise.

From Theorem 1 we can see that, if we set $K=0$, i.e., do not use the Trotter
formula, the PTSc algorithm is also valid by pairing $F_{0,1}$ with $I$. In
this case, the gate complexity is independent of the sparsity $L$ while
quadratically dependent on $t$, similar to the algorithm in Ref. Wan et al.
(2022). On the other hand, if we use $K$th-order Trotter formula, the PTSc
algorithm is $L$-dependent and has an almost linear $t$ dependence. In both
cases, the PTSc algorithms enjoy a high simulation accuracy $\varepsilon$. We
expect the $0$th-order algorithm to be especially useful for the quantum
chemistry Hamiltonian with large $L$, while higher-order algorithms to be
useful for a generic $L$-sparse Hamiltonian with a long simulation time $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 nested-
commutator form like $\mathcal{A}_{com}^{(s-1)}(H_{L},...,H_{1};H)$ in Eq. 6.
For some physical Hamiltonians with locality constraints, the number of
nonzero nested commutators in Eq. 6 can be significantly reduced.

We propose the $K$th-order nested-commutator compensation (NCc) algorithm that
utilizes the detailed nested-commutator forms of Trotter error terms from
$K+1$ to $2K+1$ order. 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)$ from
$(K+1)$th order to $(2K+1)$th order, which shrinks the error from
$\mathcal{O}(x^{K+1})$ to $\mathcal{O}(x^{2K+2})$ in one slice with the
sampling cost $\mu=1+O(\frac{\|C_{K}\|_{1}}{K+1!}x^{2K+2})$. The gate
complexity estimation is then converted to the calculation of $\|C_{K}\|_{1}$.
For instance, if we consider the $n$-qubit lattice Hamiltonian models
$H=\sum_{j=0}^{n-1}H_{j,j+1}$ where the summand $H_{j,j+1}$ acts on the $j$th
and $(j+1)$th qubits, then we can prove that $\|C_{K}\|_{1}=\mathcal{O}(n)$.
This helps to reduce the system-size dependence of NCc algorithms. More
details are provided in Appendix G.

Theorem 2 (Informal, see Theorem 2 in Appendix G) 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. In Theorem 2, we
only present the gate complexity result for the lattice Hamiltonian due to its
simplicity. In Appendix G, we also study the performance of NCc algorithms
under a generic Hamiltonian by expanding the Trotter error to operators
$C_{s}$ with different order $s$. Combining with the former calculation of
nested commutators for different Hamiltonian models Childs et al. (2021), we
can easily extend the analysis to many typical Hamiltonian models.

In Table 1, we compare the properties of $0$th-order PTSc, $K$th-order PTSc,
and $K$th-order NCc algorithms with the random-sampling implementation of LCU
and compare the results to former algorithms. The PTSc gate counts are able to
achieve logarithmic accuracy dependence $\log(1/\varepsilon)$ while the NCc
gate counts are able to achieve a better system-size dependence. Both types of
Trotter-LCU algorithms can be realized with random-sampling implementation
using only one ancillary qubit.

Algorithm | Hardness | Accuracy | Size scaling
---|---|---|---
$K$th-order Trotter Suzuki (1990) | Easy | $\mathcal{O}(\varepsilon^{-1/K})$ | $\mathcal{O}(n^{1+o(1)})$
Post-Trotter Berry et al. (2015); Low and Chuang (2019) | Hard | $\mathcal{O}(\tilde{\log}(1/\varepsilon))$ | $\mathcal{O}(n^{2})$
$0$th-order PTSc | Easy | $\mathcal{O}(\tilde{\log}(1/\varepsilon))$ | $\mathcal{O}(n^{2})$
$K$th-order PTSc | Easy | $\mathcal{O}(\tilde{\log}(1/\varepsilon))$ | $\mathcal{O}(n^{2+o(1)})$
$K$th-order NCc | Easy | $\mathcal{O}(\varepsilon^{-1/(2K+1)})$ | $\mathcal{O}(n^{1+o(1)})$
Table 1: Compare Trotter-LCU methods with former algorithms. Here,
$\tilde{\log}(x):=\log(x)/\log\log(x)$. The (implementation) hardness refers
to whether one needs to implement multi-controlled gates with plenty of
ancillary qubits. The algorithms with commutator scaling can use advanced
analysis to reduce the gate complexity of $n$-qubit lattice Hamiltonians from
$\mathcal{O}(n^{2})$ to $\mathcal{O}(n)$. More details of the comparisons and
gate complexities can be found in Appendix C.

Numerical results. We estimate the gate cost of Trotter-LCU algorithms and
compare it with the state-of-the-art Trotter algorithm, i.e., $4$th-order
Trotter Suzuki (1991); Childs et al. (2019, 2021). We compile their quantum
circuits to CNOT gates, single-qubit Clifford gates, and single-qubit $Z$-axis
rotation gates $R_{z}(\theta)=e^{i\theta Z}$. Since the non-Clifford $Z$-axis
rotation gates $R_{z}(\theta)$ are the most resource-consuming part on a
fault-tolerant quantum computer, here we mainly compare the $R_{z}(\theta)$
gate number in the main text. The CNOT gate number comparison results are
similar, shown in Appendix L. 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}$. Fig. 2(a,b) show the gate counts for
the $4$th-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 $4$th-order Trotter with random
permutation is based on the analytical bounds in Ref. Childs et al. (2019).

Figure 2:  Non-Clifford $R_{z}(\theta)$ gate number estimation for simulating
real-time dynamics with Trotter and Trotter-LCU algorithms. (a) and (b) show
the $R_{z}(\theta)$ gate counts when simulating the Hamiltonian
$H=\sum_{i,j}X_{i}X_{j}+\sum_{i}Z_{i}$. (a) shows the $R_{z}(\theta)$ gate
count with an increasing time and fixed system-size $n=20$. (b) shows the
$R_{z}(\theta)$ gate count with an increasing system size with the time $t=n$.
We list the performance of $0$th, $2$nd, and $4$th-order PTSc algorithms with
$\mu=2$. The $4$th-order Trotter analytical bound is from Ref. Childs et al.
(2019). (c) The $R_{z}(\theta)$ gate count for the nearest-interaction
Heisenberg model using the nested commutator bound with $n=12$ and $n=50$. The
$4$th-order Trotter commutator bound is from Ref. Childs et al. (2021).

From Fig. 2(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 $4$th-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 $4$th-order Trotter ($O(t^{1.25})$), $2$nd-order PTSc
($\mathcal{O}(t^{1.2})$) and $4$th-order PTSc ($\mathcal{O}(t^{1.11})$)
algorithms in Fig. 2(a). For a short-time evolution $t=n$, the system-size
dependence of $2$nd- or $4$th-order PTSc in Fig. 2(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 $4$th-order PTSc to $0$th-order PTSc will be
more obvious. The composition of Trotter and LCU formulas enables $2$nd or
$4$th-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}$, where
$\vec{\sigma}_{i}:=(X_{i},Y_{i},Z_{i})$ is the vector of Pauli operators on
the $i$th qubit, using the nested commutator bounds. In Fig. 2(c), we choose
$n=t=12$ and $50$ and show the gate number with respect to the accuracy
requirement $\varepsilon$. The $4$th-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
$2$nd-order NCc algorithm is based on the analytical bound in Appendix H. We
mainly present results for the $2$nd-order NCc algorithm due to its simplicity
and leave precise higher-order NCc gate count analysis for future study.

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

Summary and outlook. We develop a unified framework for the Trotter-LCU
algorithms with a series connection of Trotter and LCU algorithms. In fact,
both Trotter and LCU algorithms can be regarded as special cases in the
framework. This provides a possible composition method between different
Hamiltonian simulation algorithms Hagan and Wiebe (2022). Recently, Cho et al.
Cho et al. (2022) consider similar idea to compensate the Trotter error using
randomized unitary operators. There are still many possible variants to be
explored. For example, we may replace the deterministic Trotter algorithms
with the ones with random permutation Childs et al. (2019). We may also
consider the usage of quantum signal processing to compensate the Trotter
errors. Using anti-commutative cancellation Zhao and Yuan (2021), we can
further reduce the compensation terms. Furthermore, we may consider using LCU
to compensate both Trotter and physical errors to pursue a robust Hamiltonian
simulation.

###### Acknowledgements.

We thank Xiaoming Zhang, Xiao Yuan, Min-Hsiu Hsieh, Yuan Su, Kaiwen Gui, Ming
Yuan, Senrui Chen and Ying Li for helpful discussion and suggestions. P. Z.
and L. J. acknowledge support from the ARO MURI (W911NF-21-1-0325), AFOSR MURI
(FA9550-19-1-0399, FA9550-21-1-0209), AFRL (FA8649-21-P-0781), DoE Q-NEXT, NSF
(OMA-1936118, ERC-1941583, OMA-2137642), NTT Research, and the Packard
Foundation (2020-71479). Q. Z. acknowledges the support from HKU Seed Fund for
New Staff.

## References

  * Feynman et al. (1982) R. P. Feynman et al., Int. j. Theor. phys 21 (1982), URL http://dx.doi.org/10.1007/BF02650179.
  * Abrams and Lloyd (1999) D. S. Abrams and S. Lloyd, Phys. Rev. Lett. 83, 5162 (1999), URL https://link.aps.org/doi/10.1103/PhysRevLett.83.5162.
  * Aspuru-Guzik et al. (2005) A. Aspuru-Guzik, A. D. Dutoi, P. J. Love, and M. Head-Gordon, Science 309, 1704 (2005), eprint https://www.science.org/doi/pdf/10.1126/science.1113479, URL https://www.science.org/doi/abs/10.1126/science.1113479.
  * Farhi et al. (2014) E. Farhi, J. Goldstone, and S. Gutmann, _A quantum approximate optimization algorithm_ (2014), URL https://arxiv.org/abs/1411.4028.
  * Zhou et al. (2020) L. Zhou, S.-T. Wang, S. Choi, H. Pichler, and M. D. Lukin, Phys. Rev. X 10, 021067 (2020), URL https://link.aps.org/doi/10.1103/PhysRevX.10.021067.
  * Harrow et al. (2009) A. W. Harrow, A. Hassidim, and S. Lloyd, Phys. Rev. Lett. 103, 150502 (2009), URL https://link.aps.org/doi/10.1103/PhysRevLett.103.150502.
  * Lloyd (1996) S. Lloyd, Science 273, 1073 (1996), eprint https://www.science.org/doi/pdf/10.1126/science.273.5278.1073, URL https://www.science.org/doi/abs/10.1126/science.273.5278.1073.
  * Suzuki (1990) M. Suzuki, Physics Letters A 146, 319 (1990), ISSN 0375-9601, URL https://www.sciencedirect.com/science/article/pii/037596019090962N.
  * Suzuki (1991) M. Suzuki, Journal of Mathematical Physics 32, 400 (1991), URL https://aip.scitation.org/doi/10.1063/1.529425.
  * Berry et al. (2007) D. W. Berry, G. Ahokas, R. Cleve, and B. C. Sanders, Communications in Mathematical Physics 270, 359 (2007), URL https://doi.org/10.1007/s00220-006-0150-x.
  * Campbell (2019) E. Campbell, Phys. Rev. Lett. 123, 070503 (2019), URL https://link.aps.org/doi/10.1103/PhysRevLett.123.070503.
  * Childs et al. (2019) A. M. Childs, A. Ostrander, and Y. Su, Quantum 3, 182 (2019), ISSN 2521-327X, URL https://doi.org/10.22331/q-2019-09-02-182.
  * Childs and Su (2019) A. M. Childs and Y. Su, Phys. Rev. Lett. 123, 050503 (2019), URL https://link.aps.org/doi/10.1103/PhysRevLett.123.050503.
  * Childs et al. (2021) A. M. Childs, Y. Su, M. C. Tran, N. Wiebe, and S. Zhu, Phys. Rev. X 11, 011020 (2021), URL https://link.aps.org/doi/10.1103/PhysRevX.11.011020.
  * Reiher et al. (2017) M. Reiher, N. Wiebe, K. M. Svore, D. Wecker, and M. Troyer, Proceedings of the National Academy of Sciences 114, 7555 (2017), eprint https://www.pnas.org/doi/pdf/10.1073/pnas.1619152114, URL https://www.pnas.org/doi/abs/10.1073/pnas.1619152114.
  * Berry et al. (2014) D. W. Berry, A. M. Childs, R. Cleve, R. Kothari, and R. D. Somma, in _Proceedings of the forty-sixth annual ACM symposium on Theory of computing_ (ACM, 2014), URL https://doi.org/10.1145%2F2591796.2591854.
  * Berry et al. (2015) D. W. Berry, A. M. Childs, R. Cleve, R. Kothari, and R. D. Somma, Phys. Rev. Lett. 114, 090502 (2015), URL https://link.aps.org/doi/10.1103/PhysRevLett.114.090502.
  * Low and Chuang (2019) G. H. Low and I. L. Chuang, Quantum 3, 163 (2019), ISSN 2521-327X, URL https://doi.org/10.22331/q-2019-07-12-163.
  * Low and Chuang (2017) G. H. Low and I. L. Chuang, Phys. Rev. Lett. 118, 010501 (2017), URL https://link.aps.org/doi/10.1103/PhysRevLett.118.010501.
  * Childs and Wiebe (2012) A. M. Childs and N. Wiebe, Quantum Information and Computation 12 (2012), ISSN 1533-7146, URL http://dx.doi.org/10.26421/QIC12.11-12.
  * Lin and Tong (2021) L. Lin and Y. Tong, _Heisenberg-limited ground state energy estimation for early fault-tolerant quantum computers_ (2021), eprint 2102.11340, URL https://arxiv.org/abs/2102.11340.
  * Yang et al. (2021) Y. Yang, B.-N. Lu, and Y. Li, PRX Quantum 2, 040361 (2021), URL https://link.aps.org/doi/10.1103/PRXQuantum.2.040361.
  * Wan et al. (2022) K. Wan, M. Berta, and E. T. Campbell, Phys. Rev. Lett. 129, 030503 (2022), URL https://link.aps.org/doi/10.1103/PhysRevLett.129.030503.
  * Faehrmann et al. (2021) P. K. Faehrmann, M. Steudtner, R. Kueng, M. Kieferova, and J. Eisert, _Randomizing multi-product formulas for improved hamiltonian simulation_ (2021), URL https://arxiv.org/abs/2101.07808.
  * Lin and Tong (2020) L. Lin and Y. Tong, Quantum 4, 372 (2020), ISSN 2521-327X, URL https://doi.org/10.22331/q-2020-12-14-372.
  * Zeng et al. (2021) P. Zeng, J. Sun, and X. Yuan, _Universal quantum algorithmic cooling on a quantum computer_ (2021), URL https://arxiv.org/abs/2109.15304.
  * Zhang et al. (2022) R. Zhang, G. Wang, and P. Johnson, Quantum 6, 761 (2022), ISSN 2521-327X, URL https://doi.org/10.22331/q-2022-07-11-761.
  * Kitaev (1995) A. Y. Kitaev, _Quantum measurements and the abelian stabilizer problem_ (1995), URL https://arxiv.org/abs/quant-ph/9511026.
  * Hagan and Wiebe (2022) M. Hagan and N. Wiebe, _Composite quantum simulations_ (2022), URL https://arxiv.org/abs/2206.06409.
  * Cho et al. (2022) C. H. Cho, D. W. Berry, and M.-H. Hsieh, _Doubling the order of approximation via the randomized product formula_ (2022), URL https://arxiv.org/abs/2210.11281.
  * Zhao and Yuan (2021) Q. Zhao and X. Yuan, Quantum 5, 534 (2021), ISSN 2521-327X, URL https://doi.org/10.22331/q-2021-08-31-534.
  * Childs et al. (2018) A. M. Childs, D. Maslov, Y. Nam, N. J. Ross, and Y. Su, Proceedings of the National Academy of Sciences 115, 9456 (2018), URL https://www.pnas.org/content/115/38/9456.
  * Low et al. (2019) G. H. Low, V. Kliuchnikov, and N. Wiebe, _Well-conditioned multiproduct hamiltonian simulation_ (2019), URL https://arxiv.org/abs/1907.11679.
  * Zeng et al. (2022) P. Zeng, Q. Zhao, J. Sun, S. Chen, and L. Jiang, _Indirect shadow estimation for hamiltonian simulation_ (2022), in preparation.
  * Aaronson and Gottesman (2004) S. Aaronson and D. Gottesman, Phys. Rev. A 70, 052328 (2004), URL https://link.aps.org/doi/10.1103/PhysRevA.70.052328.
  * Schuch et al. (2008) N. Schuch, M. M. Wolf, K. G. H. Vollbrecht, and J. I. Cirac, New Journal of Physics 10, 033032 (2008), URL https://doi.org/10.1088/1367-2630/10/3/033032.
  * Knapp (2005) A. W. Knapp, _Basic real analysis_ (Springer Science & Business Media, 2005), URL https://doi.org/http://dx.doi.org/10.3792/euclid/9781429799997.
  * Canonne (2016) C. Canonne, _A short note on poisson tail bounds_ (2016), URL http://www.cs.columbia.edu/~ccanonne/files/misc/2017-poissonconcentration.pdf.
  * Hoorfar and Hassani (2008) A. Hoorfar and M. Hassani, J. Inequal. Pure and Appl. Math 9, 5 (2008), URL https://www.emis.de/journals/JIPAM/article983.html?sid=983.
  * Babbush et al. (2018) R. Babbush, C. Gidney, D. W. Berry, N. Wiebe, J. McClean, A. Paler, A. Fowler, and H. Neven, Phys. Rev. X 8, 041015 (2018), URL https://link.aps.org/doi/10.1103/PhysRevX.8.041015.
  * Sun et al. (2022) J. Sun, X. Yuan, M. Kim, V. Vlatko, and P. Zeng, _High-precision eigenenergy estimation and eigenstate property estimation: a framework and resource estimation_ (2022), in preparation.
  * Bocharov et al. (2015) A. Bocharov, M. Roetteler, and K. M. Svore, Phys. Rev. Lett. 114, 080502 (2015), URL https://link.aps.org/doi/10.1103/PhysRevLett.114.080502.

###### Contents

  1. A Outline
  2. B Related works
  3. C Summary of results
  4. D Preliminaries
    1. D.1 Trotter formulas
    2. D.2 Linear combination of unitaries
    3. D.3 Construct LCU formulas for Trotter remainders
  5. E Trotter-LCU formula with direct Taylor-series compensation
    1. E.1 Zeroth-order case
    2. E.2 First-order case
    3. E.3 Second- and higher-order case
  6. F Trotter-LCU formula with paired Taylor-series compensation
    1. F.1 Zeroth-order case
    2. F.2 First-order case
    3. F.3 Second- and higher-order case
  7. G Trotter-LCU formula with nested-commutator compensation
    1. G.1 Derivation of the nested-commutator compensation formula
    2. G.2 First-order case (lattice Hamiltonians)
    3. G.3 Second- and higher-order case (lattice Hamiltonians)
    4. G.4 General Hamiltonians
  8. H Tighter nested-commutator analysis for second-order Trotter-LCU algorithm
    1. H.1 Bound the 1-norm of recurrence function
    2. H.2 Estimating the gate counts
  9. I Explicit nested-commutator compensation for higher-order Trotter remainders
    1. I.1 Derivation of LCU formula with nested-commutator form
    2. I.2 Norm bounds for LCU formula
  10. J Explicit nested-commutator compensation for generic Hamiltonians
    1. J.1 Derivation of LCU formula in nested-commutator form
    2. J.2 Norm bounds for LCU formula
  11. K Some useful formulas in the proof
  12. L Additional numerical results
    1. L.1 Generic L-sparse Hamiltonians
    2. L.2 Lattice Hamiltonians

##  Appendix A Outline

We give a synopsis of the appendix here. In Appendix B, we review the related
historical works of the Trotter, LCU and composte algorithms. In Appendix C,
we summarize the gate complexity of different Trotter-LCU algorithms and
compare it with former algorithms.

In Appendix D, we introduce the background knowledges related to the Trotter
and LCU formulas. We review their definition and their gate complexity. With
these information, in Appendix D.3 we introduce three methods to construct LCU
formulas for the Trotter remainders: the direct Taylor-series compensation
(DTSc), paired Taylor-series compensation (PTSc), and nested-commutator
compensation (NCc).

In Appendix E, F, and G, we introduce the DTSc, PTSc, and NCc Trotter-LCU
methods in detail. We derive the explicit LCU formulas, propose detailed
random-sampling procedures and analyze the performance of the LCU formulas. We
will first analyze the LCU formula construction for the zeroth-order and
first-order Trotter remainders and then generalize it to general $K$th-order
case. Proofs of Theorem 1 and 2 are provided in Appendix F.3 and Appendix G.3,
respectively. As the complementary of Appendix G, Appendix H, I, and J
provides more details realted to the NCc methods.

Finally, in Appendix K we provide some useful formulas in the proof. In
Appendix L, we provide more numerical results, including the gate cost
comparison with the quantum signal processing (QSP) algorithm.

##  Appendix B Related works

Here, we summarize the related historical works of the Trotter, LCU
algorithms, and their composite ones.

Trotter formulas. There has been a long history for the study of Trotter-type
formulas. We refer the readers to Sec. IB in Ref. Childs et al. (2021) for the
history of Trotter formulas. The Trotter-type formulas are first analyzed in
the study of Lie group in 19th century. In 1990, Suzuki extends the Trotter
formulas to the high-order case and finds their application in the path-
integral formulation of quantum mechanics Suzuki (1990, 1991). In 1996, Lloyd
considers the usage of first-order Trotter formula in Hamiltonian simulation
Lloyd (1996). In 2007, Berry et. al. provide a detailed error analysis Berry
et al. (2007) for the Trotter formulas. After that, in 2018, Childs et. al.
tighten the Trotter error analysis using the using the commutator information
in the Hamiltonian Childs et al. (2018) and derived the nested commutator
forms in the following analysis Childs and Su (2019); Childs et al. (2021). To
improve the Hamiltonian sparsity dependence, Campbell proposes the qDRIFT
algorithm which randomly sample the summands to realize the product formulas
Campbell (2019).

Linear-combination-of-unitary. In 2012, Childs and Wiebe clarify the idea of
linear-combination-of-unitary in Hamiltonian simulation Childs and Wiebe
(2012). They design LCU formulas based on the summation of different Trotter
formulas and consider a non-deterministic realization. The authors adopt
multi-product formulas $M_{R,2k}(t)$ with $R$ times of different $2k$th-order
Trotter formulas and carefully design the linear combination coefficients to
cancel out the low-order time-dependent terms $t^{q}$ with $q<2R+1$ in the
additive errors. In 2019, Low et al. improves the multi-product formula
construction by reducing the 1-norm of the LCU of $K$th-order formula from
$e^{\Omega}(K)$ to $\mathcal{O}(\log K)$ Low et al. (2019). They also discuss
the commutator scaling of the multi-product formula. Besides the multi-product
formula construction, in 2015, Berry et. al. introduce new LCU formula based
on Taylor-series expansion Berry et al. (2015). By introducing the oblivious
amplitude amplification Berry et al. (2014), they make the LCU-based
simulation algorithm nearly-deterministic and prove the gate complexity to be
logarithmic to inverse accuracy $1/\varepsilon$. On the other hand, there have
been discussion about realizing LCU formula based on random sampling. In 2021,
Wan et al. consider the realization of Taylor-series expansion of $U(t)$ based
on random sampling Wan et al. (2022). Combined with the random-sampling
version of phase estimation algorithm by Lin and Tong Lin and Tong (2021),
they propose a randomized phase-estimation algorithm whose gate number in each
sample is independent of $L$. At the same time, Faehrmann et al. propose new
multi-product LCU formulas Faehrmann et al. (2021) following Ref. Childs and
Wiebe (2012) and show that the simulation error can be exponentially
suppressed by increasing the number of Trotter components $R$ in the LCU
formula. They consider a realization of multi-product LCU formula based on
random sampling.

Composite Hamiltonian simulation. There have been discussion about combining
different Hamiltonian simulation algorithms. In 2021, Yang et al. propose
Hamiltonian simulation algorithms by randomly sampling the terms in Taylor-
series expansion of the remainder in the Trotter formulas Yang et al. (2021).
They mainly work on the asymptotic performance without considering the effect
of finite gate number. We generalize the discussion and systematically study
different Trotter-LCU algorithms and provide rigorous gate complexity
analysis. In 2022, Hegam and Wiebe propose a composite Hamiltonian simulation
method by combining deterministic Trotter formulas with qDRIFT methods Hagan
and Wiebe (2022). They also discuss the composite algorithms between different
types of Hamiltonian simulation algorithm, such as Trotter and LCU. We remark
that, the concept of “composition” there is different from our discussion —
They consider to perform different types of Hamiltonian simulation algorithms
for different disjoint physical subspaces.

##  Appendix C Summary of results

In this work, we discuss various types of Trotter-LCU methods. We consider
three methods to construct the LCU formulas: direct Taylor-series compensation
(DTSc) (Sec. E), paired Taylor-series compensation (PTSc) (Sec. F), and
nested-commutator compensation (NCc) (Sec. G). We consider two ways to
implement the LCU formula for the compensation of Trotter remainders. In the
random-sampling implementation Faehrmann et al. (2021); Wan et al. (2022), we
sample the elementary unitaries $\\{V_{i}\\}$ based on the corresponding
probability $\\{p_{i}\\}$ and embed the realization of the LCU formula into a
Hadamard test Kitaev (1995). In the coherent implementation Berry et al.
(2014, 2015), we introduce $m=\lceil\log_{2}\Gamma\rceil$ ancillary qubits to
realize the superposition of different $\\{V_{i}\\}$. Our major focus is on
the random-sampling implementation, where we can estimate the properties of
the target state $U(t)\ket{\psi_{0}}$ with only one ancillary qubit.

We summarize the gate complexities of the algorithms with random-sampling
implementation in Table 2, while the ones for the coherent implementation in
Table 3.

Table 2: Summary of the random-sampling Trotter-LCU Hamiltonian simulation protocols and corresponding gate complexities. Here, $\tilde{\log}(x):=\log(x)/\log\log(x)$. The direct Taylor-series (DTS), paired Taylor-series (PTS), and nested-commutator (NC) are discussed in Sec. E, F, and G, respecively. Since the result of NC compensation relies on the specific Hamiltonian structure, here we main provide results for the lattice Hamiltonians. Trotter order | Direct Taylor-series compensation | Paired Taylor-series compensation | Nested-commutator compensation
---|---|---|---
$0$th-order | N/A | $\mathcal{O}((\lambda t)^{2}\tilde{\log}(1/\varepsilon))$ | N/A
$1$st-order | $\mathcal{O}((\lambda t)^{2}(L+\tilde{\log}(1/\varepsilon)))$ | $\mathcal{O}((\lambda t)^{1+\frac{1}{3}}(L+\tilde{\log}(1/\varepsilon)))$ | $\mathcal{O}(n^{1+\frac{2}{3}}t^{1+\frac{1}{3}}\varepsilon^{-\frac{1}{3}})$
$2k$th-order | $\mathcal{O}((\lambda t)^{1+\frac{1}{2k}}(L+\tilde{\log}(1/\varepsilon)))$ | $\mathcal{O}((\lambda t)^{1+\frac{1}{4k+1}}(L+\tilde{\log}(1/\varepsilon)))$ | $\mathcal{O}(n^{1+\frac{2}{4k+1}}t^{1+\frac{1}{4k+1}}\varepsilon^{-\frac{1}{4k+1}})$
Table 3: Summary of the coherent Trotter-LCU Hamiltonian simulation protocols and corresponding gate complexities. Here, $\tilde{\log}(x):=\log(x)/\log\log(x)$. Since the result of NC compensation relies on the specific Hamiltonian structure, here we main provide results for the lattice Hamiltonians. Trotter order | Direct Taylor-series compensation | Paired Taylor-series compensation | Nested-commutator compensation
---|---|---|---
$0$th-order |  $\mathcal{O}((\lambda t)L\tilde{\log}(1/\varepsilon))$ Berry et al. (2015) | $\mathcal{O}((\lambda t)L\tilde{\log}(1/\varepsilon))$ | N/A
$1$st-order | $\mathcal{O}((\lambda t)L\tilde{\log}(1/\varepsilon))$ | $\mathcal{O}((\lambda t)L\tilde{\log}(1/\varepsilon))$ | $\mathcal{O}(n^{1+\frac{2}{3}}t^{1+\frac{1}{3}}\varepsilon^{-\frac{1}{3}})$
$2k$th-order | $\mathcal{O}((\lambda t)L\tilde{\log}(1/\varepsilon))$ | $\mathcal{O}((\lambda t)L\tilde{\log}(1/\varepsilon))$ | $\mathcal{O}(n^{1+\frac{2}{4k+1}}t^{1+\frac{1}{4k+1}}\varepsilon^{-\frac{1}{4k+1}})$

In Table 4, we compare the Trotter-LCU algorithms with the current advanced
algorithms. We mainly consider $0$th-order PTSc, $K$th-order PTSc, $K$th-order
NCc algorithms with random sampling implementation of LCU and the $K$th-order
NCc algorithm with coherent implementation of LCU. The $0$th-order PTSc
algorithm is independent of the number of terms $L$ in Hamiltonians, which is
especially useful for many quantum chemistry problem with many terms; the
$K$th-order PTSc algorithm owns the best time dependence of
$\mathcal{O}(t^{1+\frac{1}{2K+1}})$ over all random-sampling based algorithms,
which is suitable for a long-time high-precision simulation; the $K$th-order
NCc algorithm owns almost linear $n$ and $t$ dependence for the lattice
models.

Table 4: Compare Trotter-LCU methods with former algorithms. The implementation hardness refers to whether one needs to implement multi-controlled gates with plenty of ancillary qubits. The algorithms with commutator scaling can use advanced analysis to reduce the gate complexity of $n$-qubit lattice Hamiltonians from $\mathcal{O}(n^{2})$ to $\mathcal{O}(n)$. Algorithm | Implementation hardness | Simulation accuracy | Commutator scaling | Time dependence
---|---|---|---|---
$K$th-order Trotter Suzuki (1990) | Easy | Low, $\mathcal{O}(\varepsilon^{-1/K})$ | Yes, $\mathcal{O}(n^{1+1/K})$ | $\mathcal{O}(t^{1+1/K})$
Post-Trotter (LCU, QSP) Berry et al. (2015); Low and Chuang (2019) | Hard | High, $\mathcal{O}(\tilde{\log}(1/\varepsilon))$ | No, $\mathcal{O}(n^{2})$ | $\mathcal{O}(t)$
$0$th-order PTSc (random) | Easy | High, $\mathcal{O}(\tilde{\log}(1/\varepsilon))$ |  No, $\mathcal{O}(n^{2})$ ($L$ independent) | $\mathcal{O}(t^{2})$
$K$th-order PTSc (random) | Easy | High, $\mathcal{O}(\tilde{\log}(1/\varepsilon))$ | No, $\mathcal{O}(n^{2+\frac{1}{2K+1}})$ | $\mathcal{O}(t^{1+\frac{1}{2K+1}})$
$K$th-order NCc (random) | Easy | Medium, $\mathcal{O}(\varepsilon^{-1/(2K+1)})$ | Yes, $\mathcal{O}(n^{1+\frac{2}{2K+1}})$ | $\mathcal{O}(t^{1+\frac{1}{2K+1}})$
$K$th-order NCc (coherent) | Hard | Medium, $\mathcal{O}(\varepsilon^{-1/(2K+1)})$ | Yes, $\mathcal{O}(n^{1+\frac{2}{2K+1}})$ | $\mathcal{O}(t^{1+\frac{1}{2K+1}})$

##  Appendix D Preliminaries

In this work, we focus on the construction of $U(t):=e^{-iHt}$ of a
Hamiltonian $H$ on $n$ qubits. An $n$-qubit Hamiltonian $H$ is a Hermitian
operator which can be written as

$H=\sum_{l=1}^{L}H_{l}=\sum_{l=1}^{L}\alpha_{l}P_{l}=\lambda\sum_{l=1}^{L}p_{l}P_{l},$
(9)

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}$. For simplicity, we also assume that
$t>0$ in the following discussion. If $t<0$, we will absorb the signs into the
Pauli operators, which will not affect the analysis.

In what follows, we will denote the properties related to the summands of
Hamiltonians $\\{H_{l}\\}_{l=1}^{L}$ to be vectors,

$\displaystyle\vec{\alpha}$
$\displaystyle:=[\alpha,\alpha_{2},...,\alpha_{L}],$ (10)
$\displaystyle\vec{p}$ $\displaystyle:=[p_{1},p_{2},...,p_{L}],$
$\displaystyle\vec{H}$ $\displaystyle:=[H_{1},H_{2},...,H_{L}],$
$\displaystyle\vec{P}$ $\displaystyle:=[P_{1},P_{2},...,P_{L}].$

In Eq. 9, $\lambda=|\vec{\alpha}|=\sum_{l=1}^{L}\alpha_{l}$ is the
$l_{1}$-norm of the coefficient vector $\vec{\alpha}$. Then $\vec{p}$ will be
a probability distribution.

In all the Hamiltonian simulation algorithms discussed in this work, we divide
the real-time evolution $U(t)$ into $\nu$ segments,

$U(t)=e^{-iHt}=(U(x))^{\nu}=\left(e^{-iHx}\right)^{\nu},$ (11)

where $x:=t/\nu$. In this way, we focus on the construction of $U(x)$, which
is usually simpler. Since we assume that $t>0$, $x$ will also be positive.
This will not affect the generality of our results as it can be easily
extended to the case with $t<0$.

Here, we introduce several notation to simplify the summation or
multiplication of similar terms. For a group of index
$\\{a_{1},a_{2},...,a_{m}\\}$ of the same type with length $m$, we denote it
as $a_{1:m}$. The summation of the indices is

$\sum a_{1:m}:=\sum_{i=1}^{m}a_{i}.$ (12)

Moreover, for the multiplication of the quantities
$O_{a_{1}}O_{a_{2}}...O_{a_{m}}$ with the same type of indices, we denote it
as $O^{(m)}_{a_{1:m}}$, i.e.,

$\displaystyle O^{(m)}_{a_{1:m}}$ $\displaystyle:=\prod_{i=1}^{m}O_{a_{i}},$
(13) $\displaystyle\sum_{a_{1:m}}O^{(m)}_{a_{1:m}}$
$\displaystyle:=\sum_{a_{1},a_{2},...,a_{m}}O^{(m)}_{a_{1:m}}.$

When $\\{O_{i}\\}_{i}$ is a group of operators, the value of
$O^{(m)}_{a_{1:m}}$ may be unequal to $O^{(m)}_{a_{m:1}}$ with the inverse
multiplication direction.

Especially, for the indices $l$ related to the Hamiltonian sparsity, i.e., $l$
is enumerated from $1$ to $L$, we adopt the vector notations mentioned in Eq.
10. Moreover, when multiple $l$-related vectors occur in the same summation or
product formula, we will treat them as index-wise operation, i.e.,

$\displaystyle\sum f(\vec{\alpha},\vec{r})$
$\displaystyle=\sum_{l=1}^{L}f(\alpha_{l},r_{l}),$ (14)
$\displaystyle\prod^{\leftarrow}f(\vec{\alpha},\vec{r})$
$\displaystyle=\prod_{l=1}^{L}f(\alpha_{l},r_{l}).$

When the output of the function $f(\vec{\alpha},\vec{r})$ is a number, we will
omit the arrow in the product as the product direction becomes insignificant.

###  D.1 Trotter formulas

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

$\displaystyle
S_{1}(x)=\prod_{l=1}^{L}e^{-ixH_{l}}=e^{-ixH_{L}}...e^{-ixH_{2}}e^{-ixH_{1}}=:\prod^{\leftarrow}e^{-ix\vec{H}}.$
(15)

Here, $\vec{H}$ is defined in Eq. 10. In Eq. 15, 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$. On the other hand,
the Hermitian conjugation of $S_{1}(x)$ is,

$\displaystyle
S_{1}(x)^{\dagger}=\prod_{l=L}^{1}e^{ixH_{l}}=e^{ixH_{1}}e^{ixH_{2}}...e^{ixH_{L}}=:\prod^{\rightarrow}e^{ix\vec{H}}.$
(16)

The second-order Trotter formula is

$\displaystyle S_{2}(x)$
$\displaystyle=D_{2}(-x)^{\dagger}D_{2}(x)=S_{1}(-\frac{x}{2})^{\dagger}S_{1}(\frac{x}{2})$
(17)
$\displaystyle=\prod^{\rightarrow}e^{-i\frac{1}{2}x\vec{H}}\prod^{\leftarrow}e^{-i\frac{1}{2}x\vec{H}},$

where $D_{2}(x):=S_{1}(\frac{x}{2})$. From Eq. 17 we have
$S_{2}(-x)^{\dagger}=S_{2}(x)$.

The general $2k$th-order Trotter formula is

$\displaystyle S_{2k}(x)$ $\displaystyle=D_{2k}(-x)^{\dagger}D_{2k}(x)$ (18)
$\displaystyle=[S_{2k-2}(u_{k}x)]^{2}S_{2k-2}((1-4u_{k})x)[S_{2k-2}(u_{k}x)]^{2},$

with $u_{k}:=1/(4-4^{1/(2k-1)})$ for $k\geq 1$. The operator $D_{2k}(x)$ is
defined recursively,

$D_{2k}(x):=D_{2k-2}((1-4u_{k})x)S_{2k-2}(u_{k}x)^{2}.$ (19)

By induction from Eq. 18 and Eq. 19 we can show that
$S_{2k}(-x)^{\dagger}=S_{2k}(x)$. The operator $D_{2k}(x)$ will take the
following form,

$D_{2k}(x)=\prod^{\leftarrow}e^{-iv_{\kappa}x\vec{H}}\prod^{\rightarrow}e^{-iv_{\kappa-1}x\vec{H}}...\prod^{\rightarrow}e^{-iv_{2}x\vec{H}}\prod^{\leftarrow}e^{-iv_{1}x\vec{H}},$
(20)

where $\kappa=5^{k-1}$, $\\{v_{i}\\}_{i=1}^{\kappa}$ are positive numbers
indicating the segment length, which are determined by Eq. 18 and Eq. 19 and
satisfy $\sum_{i=1}^{\kappa}v_{i}=\frac{1}{2}$.

In what follows, we also denote the zeroth-order Trotter formula to be

$S_{0}(x)=I,$ (21)

for consistency. We denote the multiplicative remainder of the Trotter
formulas as

$\displaystyle V_{0}(x)$ $\displaystyle=U(x),$ (22) $\displaystyle V_{1}(x)$
$\displaystyle=U(x)S_{1}(x)^{\dagger},$ $\displaystyle V_{2k}(x)$
$\displaystyle=U(x)S_{2k}(x)^{\dagger}.$

In Ref. Suzuki (1990), Suzuki proves the following order condition for the
Trotter formulas,

$S_{K}(x)=U(x)+\mathcal{O}(x^{K+1})=e^{-iHx+\mathcal{O}(x^{K+1})},$ (23)

for $K=1$ or even positive $K$. As a result, the remainder $V_{K}(x)$ will
only contain the terms of $x^{q}$ with $q\geq K+1$. We will use the following
order condition of Trotter formulas in the later discussion.

######  Lemma 1 (Order condition of Trotter formulas, Thm 4 in Childs et al.
(2021)).

Let $H$ be an Hermitian operator and $U(x):=e^{-iHx}$ to be the time-evolution
operator with $x\in\mathbb{R}$. For the Trotter formula $S_{K}(x)$ defined in
Eq.Eq. 15 and Eq. 18, where $K=1$ or positive even number, we have

  1. 1.

Additive error: $A_{K}(x):=S_{K}(x)-U(x)=\mathcal{O}(x^{K+1})$,

  2. 2.

Multiplicative error:
$M_{K}(x):=U(x)S_{K}(x)^{\dagger}-I=\mathcal{O}(x^{K+1})$,

  3. 3.

For the exponentiated error $E_{K}(x)$ such that
$S_{K}(x)=\mathcal{T}\exp\left(\int_{0}^{x}d\tau(-iH+E_{K}(\tau))\right)$, we
have $E_{K}(x)=\mathcal{O}(x^{K})$.

###  D.2 Linear combination of unitaries

Instead of decomposing $U(x)$ as a product of elementary unitaries, another
way is to decompose $U(x)$ to a summation of elementary unitaries. We now
provide a formal definition of a LCU formula.

######  Definition 1 (Childs and Wiebe (2012)).

A $(\mu,\varepsilon)$-linear-combination-of-unitary (LCU) formula of an
operator $V$ is defined to be

$\tilde{V}=\sum_{i=0}^{\Gamma-1}c_{i}V_{i}=\mu\sum_{i=0}^{\Gamma-1}\Pr(i)V_{i},$
(24)

such that the spectral norm distance $\|V-\tilde{V}\|\leq\varepsilon$. Here,
$\mu>0$ is the $l_{1}$-norm of the coefficient vector, $\Pr(i)$ is a
probability distribution over different unitaries $U_{i}$, and
$\\{U_{i}\\}_{i=0}^{\Gamma-1}$ is a group of unitaries. Here, we assume
$c_{i}>0$ for all $i$ and absorb the phase into the unitaries $V_{i}$. We call
$\mu$ the $1$-norm of this $(\mu,\varepsilon)$-LCU formula. Moreover, we
define the $1$-norm $\|V\|_{1}$ of an operator $V$ to be its smallest $1$-norm
over all possible $(\mu,0)$-LCU formulas for $V$.

Following the $1$-norm definition in Definition 1, we know that $\|U\|_{1}=1$
for any given unitary $U$. Also, it is easy to verify that this $1$-norm
definition satisfies the following properties,

$\displaystyle\|V\|_{1}$ $\displaystyle\geq 0,$ (25) $\displaystyle\|V\|_{1}$
$\displaystyle=0\Leftrightarrow V=0,$ $\displaystyle\|\alpha V\|_{1}$
$\displaystyle=|\alpha|\|V\|_{1},$ $\displaystyle\|V+W\|_{1}$
$\displaystyle\leq\|V\|_{1}+\|W\|_{1},$ $\displaystyle\|VW\|_{1}$
$\displaystyle\leq\|V\|_{1}\|W\|_{1},$

for all operators $V$ and $W$. Therefore, $\|V\|_{1}$ in Definition 1 is
indeed a well-defined matrix norm.

In this work, we will discuss two ways to implement the LCU formula, shown in
Fig. 3. The first way is to use it in a random-sampling style. In this case,
we are not able to implement the target unitary $U$ deterministically.
Instead, we focus on estimating the observable values $\mathrm{Tr}(U\rho
U^{\dagger}O)$. The second way is to implement it coherently. In this case, we
can implement $U$ and prepare the target state.

Figure 3: Schematic diagram for the Trotter-LCU algorithms. (a) We first split
the real-time evolution $U(t)$ to $\nu$ segments, each with evolution time
$x=t/\nu$. For each segment, we first implement $K$th-order Trotter formula
($K=0,1,2k,k\in\mathbb{N}_{+}$). Then, we compensate the Trotter remainder
$V_{K}(x)$ using LCU formula. (b) For a $(\mu,\varepsilon)$-LCU formula, we
consider two different implementations. For the random-sampling
implementation, we implement the controlled-$V_{i}$ gate where $V_{i}$ is
sampled based on the probability distribution $p_{i}$. For the coherent
implementation, we use three LCU implementation gates $W$ (or $W^{\dagger}$)
defined in Eq. 45 and two reflection gates $R$ defined in Eq. 48. We need
$m:=\lceil\log_{2}\Gamma\rceil$ qubits, where $\Gamma$ is the number of
elementary unitaries in the LCU formula.

We now introduce these two LCU implementation methods in detail. The usage of
random-sampling implementation is introduced in Fig. 4 in detail. We first
prepare the $\ket{+}=\frac{1}{\sqrt{2}}(\ket{0}+\ket{1})$ state on the
ancillary qubit $A$, using it to perform controlled operations

$\displaystyle\ket{0}_{A}\bra{0}\otimes
I_{B}+\ket{1}_{B}\bra{1}\otimes(U_{i})_{B},$ (26)
$\displaystyle\ket{0}_{A}\bra{0}\otimes(U_{j})_{B}+\ket{1}_{B}\bra{1}\otimes
I_{B},$

where $B$ is the target system, $U_{i}$ and $U_{j}$ are independently sampled
based on the LCU formula in Definition 1. If we want to estimate
$\mathrm{Tr}(U\rho U^{\dagger}O)$ for a given observable $O$, then we directly
measure $\mu^{2}X_{A}\otimes O_{B}$ whose expectation value is
$\mathrm{Tr}(\tilde{U}\rho\tilde{U}^{\dagger}O)$ Faehrmann et al. (2021). If
we want to learn many observables efficiently, we can adopt the shadow
estimation circuit Zeng et al. (2022) where a random unitary $V\in\mathbb{G}$
is performed, followed by a computational basis measurement. $\mathbb{G}$ is a
tomographic-complete unitary set. Based on the measurement values $a,\vec{z}$
and the unitary information $V$, we can generate a matrix-valued estimator
$\hat{\sigma}$ for the target state $\sigma=U\rho U^{\dagger}$,

$\hat{\sigma}=\mu^{2}(-1)^{a}\mathcal{M}^{-1}(V^{\dagger}\ket{\vec{z}}\bra{\vec{z}}V).$
(27)

Here, $\mathcal{M}$ is a invertible linear map determined by $\mathbb{G}$,

$\displaystyle\mathcal{M}(O)$
$\displaystyle=\mathop{\mathbb{E}}\limits_{V\in\mathbb{G}}\sum_{\vec{z}\in\mathbb{Z}^{n}_{2}}\braket{\vec{z}}{VOV^{\dagger}}{\vec{z}}\,V^{\dagger}\ket{\vec{z}}\bra{\vec{z}}V.$
(28)

Using the circuit in Fig. 4(b), we can learn the information in the target
state $\sigma$ with extra sample cost of $\mu^{2}$. As long as $\mu$ is a
constant with respect to $t$ and $\varepsilon$, the sample complexity of the
Hamiltonian simulation based on random sampling will be the same as the case
when we prepare the target state $\sigma$ by implementing $U(t)$ on $\rho$.
One can refer to Ref. Zeng et al. (2022) for more details.

Figure 4:  Random-sampling implementation of Hamiltonian simulation based on
$(\mu,\varepsilon)$-LCU formula in Definition 1. $B$ is the original system
and $A$ is an extra ancillary qubit. $U_{i}$ and $U_{j}$ are independently
sampled based on the random sampling LCU formula. (a) Estimation circuit for
the expectation value $\mathrm{Tr}[U\rho U^{\dagger}O]$ of the given
observable $O$ Faehrmann et al. (2021). (b) Circuit to generate the matrix
estimator for the target state $\sigma=U\rho U^{\dagger}$, which can be used
to learn many observables efficiently Zeng et al. (2022). $V$ is a random
unitary chosen from a preset tomographically-complete set $\mathbb{G}$. The
estimator for $\sigma$ is
$\hat{\sigma}=\mu^{2}(-1)^{a}\mathcal{M}^{-1}(V^{\dagger}\ket{\vec{z}}\bra{\vec{z}}V).$

For simplicity, we will focus on the estimation of a single observable $O$ in
this work. Once we have a $(\mu,\varepsilon)$-LCU formula, we can estimate
$\mathrm{Tr}(U\rho U^{\dagger}O)$ with the following performance guarantee.

######  Proposition 1 (Performance of the random-sampling LCU implementation
(c.f., Theorem 2 in Faehrmann et al. (2021))).

For a target unitary $U$ and its $(\mu,\varepsilon)$-LCU formula defined in
Definition 1, if we estimate the value $\braket{O}_{U}:=\mathrm{Tr}(U\rho
U^{\dagger}O)$ with an initial state $\rho$ and observable $O$ using the
circuit in Fig. 4(a) for $N$ times, then the distance between the mean
estimator value $\hat{O}$ and the true value $\braket{O}_{U}$ is bounded by

$|\hat{O}-\braket{O}_{U}|\leq\|O\|(3\varepsilon+\varepsilon_{n}),$ (29)

with successful probability $1-\delta$. Here,
$N=2\mu^{4}\ln(2/\delta)/\varepsilon_{n}^{2}$, $\|O\|$ is the spectral norm of
$O$.

###### Proof.

We first study the estimation bias, i.e., the distance between the expectation
value $\bar{O}$ and the target value $\braket{O}_{U}$. For a specific choice
of $U_{i}$ and $U_{j}$ in Fig. 4(a), the expectation value of the observable
$\mu^{2}\braket{X_{A}\otimes O_{B}}$ is

$\bar{O}_{i,j}=\frac{\mu^{2}}{2}\left(\mathrm{Tr}(U_{i}\rho
U_{j}^{\dagger}O)+\mathrm{Tr}(U_{j}\rho U_{i}^{\dagger}O)\right).$ (30)

If we randomly sample $U_{i}$ and $U_{j}$ independently in Fig. 4(a), then
based on the LCU formula in Eq. 24, the expectation value is then

$\bar{O}:=\mathbb{E}_{i,j}\bar{O}_{i,j}=\mathrm{Tr}(\tilde{U}\rho\tilde{U}^{\dagger}O).$
(31)

We can bound the distance between $\bar{O}$ and $\braket{O}_{U}$ by

$\displaystyle\quad|\bar{O}-\braket{O}_{U}|=|\mathrm{Tr}(\tilde{U}\rho\tilde{U}^{\dagger}O)-\mathrm{Tr}(U\rho
U^{\dagger}O)|$
$\displaystyle\leq|\mathrm{Tr}(\tilde{U}\rho\tilde{U}^{\dagger}O)-\mathrm{Tr}(\tilde{U}\rho
U^{\dagger}O)|+|\mathrm{Tr}(\tilde{U}\rho U^{\dagger}O)-\mathrm{Tr}(U\rho
U^{\dagger}O)|$
$\displaystyle=|\mathrm{Tr}[O\tilde{U}\rho(\tilde{U}-U)^{\dagger}]|+|\mathrm{Tr}[(\tilde{U}-U)\rho
U^{\dagger}O]|$ (32)
$\displaystyle\leq(1+\varepsilon)\varepsilon\|O\|+\varepsilon\|O\|$
$\displaystyle\leq 3\varepsilon\|O\|.$

In the forth line, we use the fact that $\|\rho\|\leq 1$, $\|U\|\leq 1$, and
$\|\tilde{U}\|\leq\|U\|+\|U-\tilde{U}\|\leq 1+\varepsilon$.

Now, we study the statistical error caused by finite $N$. The measured value
$\hat{O}_{i}$ of the observable $\mu^{2}X\otimes O$ is in the range of
$[-\mu^{2}\|O\|,\mu^{2}\|O\|]$. Using the Hoeffding bound in Lemma 6, we have
the following probability tail bound for the mean estimator
$\hat{O}:=\frac{1}{N}\sum_{i=1}^{N}\hat{O}_{i}$,

$\Pr(\left|\hat{O}-\bar{O}\right|\geq\varepsilon_{n}\|O\|)\leq
2\exp\left(-\frac{2N\varepsilon_{n}^{2}}{4\mu^{4}}\right).$ (33)

Combining Eq. 32 and Eq. 33 and setting
$N:=2\mu^{4}\ln(2/\delta)/\varepsilon_{n}^{2}$, we have

$|\hat{O}-\braket{O}_{U}|\leq|\bar{O}-\braket{O}_{U}|+\left|\hat{O}-\bar{O}\right|=\|O\|(3\varepsilon+\varepsilon_{n}),$
(34)

with successful probability $1-\delta$.

∎

From Proposition 1 we can see that, the $1$-norm $\mu$ of the LCU formula
affects the sample comlexity while the accuracy factor $\varepsilon$ introduce
extra bias in the observable estimation.

In the later discussion, we will come up with the product of LCU formulas. We
will use following proposition.

######  Proposition 2 (Product of LCU formulas).

Suppose we have a $(\mu,\varepsilon)$-LCU formula $\tilde{U}$ for a unitary
$U$, i.e.,

$\tilde{U}=\mu\sum_{i}\Pr(i)U_{i},$ (35)

with $\|U-\tilde{U}\|\leq\varepsilon$.

Then the product formula

$\tilde{U}^{\nu}=\mu^{\nu}\sum_{i_{1},i_{2},...,i_{\nu}}\Pr(i_{1})\Pr(i_{2})...\Pr(i_{\nu})U_{1}U_{2}...U_{\nu},\quad\nu\in\mathbb{N},$
(36)

is a $(\mu^{\prime},\varepsilon^{\prime})$-LCU formula for $U^{\nu}$ with

$\displaystyle\mu^{\prime}$ $\displaystyle=\mu^{\nu},$ (37)
$\displaystyle\varepsilon^{\prime}$
$\displaystyle\leq\nu\mu^{\prime}\varepsilon.$

###### Proof.

The $1$-norm is obvious. We now bound the distance $\varepsilon^{\prime}$
between $\tilde{U}^{\nu}$ and $U^{\nu}$.

$\displaystyle\|U^{\nu}-\tilde{U}^{\nu}\|$
$\displaystyle=\left\|\sum_{k=1}^{\nu}U^{k-1}(U-\tilde{U})\tilde{U}^{\nu-k}\right\|$
$\displaystyle\leq\sum_{k=1}^{\nu}\|U^{k-1}\|\|U-\tilde{U}\|\|\tilde{U}^{\nu-k}\|$
$\displaystyle\leq\nu\|U-\tilde{U}\|\sum_{k=1}^{\nu}\max\\{\|U\|,\|\tilde{U}\|\\}^{\nu-1}$
(38) $\displaystyle\leq\nu\varepsilon\mu^{\nu-1}$
$\displaystyle\leq\nu\mu^{\prime}\varepsilon.$

∎

Finally, we remark that, when there are deterministic components in the LCU
formula, we can then simplify the random-sampling implementation by removing
the controlled operation on the deterministic components. Suppose we have the
following LCU formula for an operator $V$,

$\tilde{V}=\mu\sum_{i_{1},i_{2},...,i_{\nu}}\Pr(i_{1},i_{2},...,i_{\nu})V_{i_{\nu}}W_{\nu}...V_{i_{2}}W_{2}V_{i_{1}}W_{1},$
(39)

such that $\|V-\tilde{V}\|\leq\varepsilon$, where $W_{1},W_{2},...,W_{\nu}$
are some fixed unitaries. Then according to Definition 1, Eq. 39 is a
$(\mu,\varepsilon)$-LCU formula of $V$ with the elementary unitaries to be
$V_{\vec{i}}=V_{i_{\nu}}W_{\nu}...V_{i_{2}}W_{2}V_{i_{1}}W_{1}$ for
$\vec{i}=\\{i_{1},i_{2},...,i_{\nu}\\}$. Instead of naively apply the Hadamard
test circuit in Fig. 4(a), where we need to apply the controlled version of
$V_{\vec{i}}$ twice, we can improve the circuit implementation shown in Fig.
5(c). The composite unitary realized by the circuit in Fig. 5(c) is

$\displaystyle\prod_{a=1}^{\nu}\left(\ket{0}_{A}\bra{0}\otimes(V_{j_{a}})_{B}+\ket{1}_{A}\bra{1}\otimes
I_{B}\right)\left(\ket{0}_{A}\bra{0}\otimes
I_{B}+\ket{1}_{A}\bra{1}\otimes(V_{i_{a}})_{B}\right)\left(I_{A}\otimes(W_{a})_{B}\right)$
(40) $\displaystyle=$
$\displaystyle\prod_{a=1}^{\nu}\left(\ket{0}_{A}\bra{0}\otimes(V_{j_{a}}W_{a})_{B}+\ket{1}_{A}\bra{1}\otimes(V_{i_{a}}W_{a})_{B}\right)$
$\displaystyle=$
$\displaystyle\ket{0}_{A}\bra{0}\otimes(V_{j_{\nu}}W_{\nu}...V_{j_{1}}W_{1})_{B}+\ket{1}_{A}\bra{1}\otimes(V_{i_{\nu}}W_{\nu}...V_{i_{1}}W_{1})_{B},$

which is the same as the unitary realized in Fig. 5(b). As a result, these two
implementation are equivalent. The performance guarantee in Proposition 1 also
holds for the improved implementation.

Figure 5:  Improve the random-sampling implementation of Hamiltonian
simulation based on $(\mu,\varepsilon)$-LCU formula in Eq. 39. (a) The LCU
formula of $V$ in Eq. 39 with deterministic components. (b) The naive
application of Hadamard test using the circuit in Fig. 4(a). Here, we need to
apply the controlled version of the whole elementary
$V_{\vec{i}}:=V_{i_{2}}W_{2}V_{i_{1}}W_{1}$ twice. (c) The improved
application of Hadamard test. We remove the control on the deterministic
unitaries and apply them only once.

The gate complexity of the random-sampling implementation of the Trotter-LCU
algorithm is determined by the segment number $\nu$, Trotter order $K$ and the
gate complexity of each elementary gate in the LCU formula. As shown in Fig.
3, to construct controlled-$U(t)$, we split it to $\nu$ segments. In each
segment, we need to implement $K$th-order Trotter circuits and LCU circuits.
The gate complexity of each elementary gate in the LCU circuit is determined
by the truncation order $s_{c}$ of the Taylor-series compensation. Therefore,
the gate complexity of the overall algorithm using $K$th Trotter formula
($K=0,1,2k$) is given by

$\displaystyle N_{K}^{(RS)}=\mathcal{O}(\nu(\kappa_{K}L+s_{c}))$ (41)

where

$\kappa_{K}=\begin{cases}0,&K=0,\\\ 1,&K=1,\\\ 2\cdot
5^{K/2-1},&K=2k,k\in\mathbb{N}_{+}.\end{cases}$ (42)

The second way to implement LCU formula is to use it coherently. To this end,
we need to introduce an ancillary system $A$ with the dimension $\Gamma$ which
costs $\lceil\log_{2}\Gamma\rceil$ qubits. For the LCU lemma defined in Eq.
24, we define the amplitude-encoding unitary $U_{AE}$ (also called the state-
preparation unitary) to be,

$U_{AE}\ket{0}_{A}:=\sum_{i=0}^{\Gamma-1}\sqrt{\Pr(i)}\ket{i},$ (43)

and the select gate to be

$\mathrm{Sel}(\tilde{V}):=\sum_{i=0}^{\Gamma-1}\ket{i}\bra{i}\otimes V_{i}.$
(44)

Then, the following controlled-gate $W_{AB}$ acting on ancillary $A$ and
system $B$ defines a way to realize LCU coherently when we prepare $\ket{0}$
on $A$ and measure it on computational basis to get $\ket{0}$,

$W:=(U_{AE}^{\dagger}\otimes I)\mathrm{Sel}(\tilde{V})(U_{AE}\otimes I).$ (45)

More precisely, if we set the initial state on $B$ to be $\ket{\psi}$, then we
have

$W\ket{0}_{A}\ket{\psi}_{B}=\frac{1}{\mu}\ket{0}_{A}\left(\tilde{V}\ket{\psi}\right)_{B}+\sqrt{1-\frac{1}{\mu^{2}}}\ket{\Phi_{\perp}}_{AB},$
(46)

where $\ket{\Phi_{\perp}}_{AB}$ is a state whose ancillary state is supported
in the subsapce orthogonal to $\ket{0}$.

If we perform computation-basis measurement directly on $A$, the successful
probability to obtain $\tilde{V}\ket{\psi}$ is $\frac{1}{\mu^{2}}$. To boost
the successful probability to nearly deterministic, we can introduce the
amplitude amplification techniques Berry et al. (2014). To this end, we
consider the following isometry

$\mathcal{V}:=-WRW^{\dagger}RW\left(\ket{0}\otimes I\right)_{AB}.$ (47)

Here, $R$ is a reflection over $\ket{0}$ on the system $A$,

$R=I-2\ket{0}_{A}\bra{0}.$ (48)

In what follows, we set $\mu=2$ for simplicity. If $\tilde{V}$ is a unitary,
then we have $(\bra{0}\otimes I)\mathcal{V}=\tilde{V}$. For a general
$\tilde{V}$, we can verify that

$(\bra{0}\otimes
I)\mathcal{V}=\frac{3}{\mu}\tilde{V}-\frac{4}{\mu^{3}}\tilde{V}\tilde{V}^{\dagger}\tilde{V}=\frac{3}{2}\tilde{V}-\frac{1}{2}\tilde{V}\tilde{V}^{\dagger}\tilde{V}.$
(49)

It can be proven that Childs et al. (2018), when $\varepsilon$ is small, the
distance between the resulting operator
$\frac{3}{2}\tilde{V}-\frac{1}{2}\tilde{V}\tilde{V}^{\dagger}\tilde{V}$ will
be closed to the ideal unitary $V$,

$\left\|\frac{3}{2}\tilde{V}-\frac{1}{2}\tilde{V}\tilde{V}^{\dagger}\tilde{V}-V\right\|\leq\frac{\varepsilon^{2}+3\varepsilon+4}{2}\varepsilon.$
(50)

Also, when $\varepsilon$ is small, then the successful probability to project
the ancilla to $\ket{0}$ is high and the resulting operator will be closed to
$V$. We will use the following performance bound.

######  Lemma 2 (Performance of coherent LCU implementation Childs et al.
(2018)).

For a target unitary $U$ and its $(2,\varepsilon)$-LCU formula defined in
Definition 1, if we implement the LCU based on the isometry $\mathcal{V}$
defined in Eq. 47, then if we measure the ancillary system $A$ on the
computational basis and obtain $\ket{0}_{A}$, then the distance between
implemented operation and the target unitary $U$ is bounded by

$\|(\bra{0}\otimes I)\mathcal{V}-U\|\leq\xi.$ (51)

The successful probability to obtain $\ket{0}_{A}$ is larger than
$(1-\xi)^{2}$. Here,

$\xi:=\frac{\varepsilon^{2}+3\varepsilon+4}{2}\varepsilon.$ (52)

In practice, it is not always possible to get a LCU formula for a given $U$
with $\mu=2$. When we have a LCU formula with $\mu<2$, we can always introduce
an extra ancillary qubit $A^{\prime}$ and apply a rotation

$R(\mu)\ket{0}=\frac{\mu}{2}\ket{0}+\sqrt{1-\frac{\mu^{2}}{4}}\ket{1}.$ (53)

Together with the unitary operator $W$ on $A$ and $B$, the composite operation
have the effect

$\displaystyle R(\mu)_{A^{\prime}}\otimes
W_{AB}\ket{00}_{AA^{\prime}}\ket{\psi}_{B}$ (54)
$\displaystyle=\frac{1}{2}\ket{00}_{A^{\prime}A}\tilde{V}\ket{\psi}_{B}+\frac{\sqrt{3}}{2}\ket{\Phi^{\prime}}_{A^{\prime}AB}.$

As a result, we can treat $R(\mu)_{A^{\prime}}\otimes W_{AB}$ as the new $W$
operator to apply the amplitude amplification using the same method introduced
above.

When $\mu>2$ for a LCU formula, we can split the target unitary $U$ to several
segments to reduce $\mu$. For example, for the time evolution $U(t)$, we can
split it to $\nu$ segments $U(x)^{\nu}$ so that $\mu\leq 2$ for each segment.
In this case, suppose we implement $(\mu,\varepsilon)$-LCU formula for each
segment, the overall accuracy is $\nu\xi$ and the successful probability is
$(1-\nu\xi)^{2}$ where $\xi$ is defined in Eq. 52. This can be proven by
simply modifying Lemma 2 and utilize triangle inequalities.

######  Corollary 1 (Coherent LCU implementation with segments, Theorem G. 5
in Childs et al. (2018)).

Suppose we have a $(\mu,\varepsilon)$-LCU formula $\tilde{U}$ for a unitary
$U$ with $\mu\leq 2$. Suppose we want to implement $U^{\nu}$ where $\nu$ is a
positive integer. If we implement the LCU of $U$ based on the isometry
$\mathcal{V}$ defined in Eq. 47 for $\nu$ times, then

$\|\left((\bra{0}\otimes I)\mathcal{V}\right)^{\nu}-U^{\nu}\|\leq\nu\xi.$ (55)

The successful probability to obtain $\ket{0}_{A}$ for all segments is larger
than $(1-\nu\xi)^{2}$. Here, $\xi$ is defined in Eq. 52.

The gate complexity of the coherent implementation of the Trotter-LCU
algorithm is determined by the segment number $\nu$, the Trotter order $K$,
the number of elementary unitaries $\Gamma$ and the gate complexity of each
elementary unitary in the LCU formula. The value $\Gamma$ is related to the
specific compensation method we use, which is usually proportional to $L$ and
the truncation order $s_{c}$. The gate complexity is

$\displaystyle N_{K}^{(C)}=\mathcal{O}(\nu(\kappa_{K}L+\Gamma)),$ (56)

Although we split the evolution $U(t)$ into $\nu$ segments in both random-
sampling and coherent implementation, they follow different analysis
procedures. We summarize and compare their analysis procedure in Fig. 6. In
the random-sampling implementation, we first obtain an overall
$(\mu_{tot},\varepsilon_{tot})$-LCU formula for the whole time evolution
$U(t)$ based on Proposition 2, then analyze the performance based on
Proposition 1 and optimize the value of $\nu$. As a comparison, in the
coherent implementation, we first find proper $x$ value so that we can
construct a LCU formula with $\mu=2$. After that, we calculate the segment
number $\nu=t/x$ and analyze the performance directly based on Corollary 1.
Even starting with the same LCU formula for a given segment, the different
analysis will yield different performance for these two implementation
methods.

Figure 6: Comparison of the performance analysis based on random-sampling and
coherent implementation.

###  D.3 Construct LCU formulas for Trotter remainders

The construction of LCU formulas for the Trotter remainder $V_{K}(x)$ is at
the core of the implementation and performance of the Trotter-LCU algorithms.
We need to consider the following requirements for the LCU formulas,

  1. 1.

The LCU formula should be able to achieve a high accuracy (i.e., a small
$\varepsilon$) with few terms.

  2. 2.

For the random-sampling implementation, the $1$-norm $\mu$ of the LCU formula
should be as small as possible.

  3. 3.

For the coherent implementation, the number of elementary unitaries $\Gamma$
in the LCU formula should be as small as possible.

In this section, we summarize the major idea of three different methods to
construct LCU formulas for $V_{K}(x)$, illustrated in Fig. 7. We mainly
consider different realizations of Taylor-series compensation Berry et al.
(2015). This promises us to achieve a high simulation accuracy $\varepsilon$
with few terms.

For simplicity, let us now assume that the Hamiltonian can be written as
$H=H_{1}+H_{2}$, where $H_{1}$ and $H_{2}$ are two Hermitian operators with
known $(\mu_{1},0)$\- and $(\mu_{2},0)$-LCU decomposition with Pauli
operators, respectively,

$\displaystyle H_{1}=\mu_{1}U_{1}:=\mu_{1}\sum_{i}q_{1}^{(i)}P_{1}^{(i)},$
(57) $\displaystyle
H_{2}=\mu_{2}U_{2}:=\mu_{2}\sum_{i}q_{2}^{(i)}P_{2}^{(i)},$

where $\\{q_{1}^{(i)}\\}_{i}$ and $\\{q_{2}^{(i)}\\}_{i}$ are two probability
distributions; $\\{P_{1}^{(i)}\\}_{i}$ and $\\{P_{2}^{(i)}\\}_{i}$ are two
groups of Pauli operators. We denote $\lambda:=\mu_{1}+\mu_{2}$. Therefore,

$H=\lambda(p_{1}U_{1}+p_{2}U_{2})=\lambda\sum_{l=1,2}p_{l}\sum_{i}q_{l}^{(i)}P_{l}^{(i)},$
(58)

which is a $(\lambda,0)$-LCU formula of $H$. Here, $p_{l}:=\mu_{l}/\lambda$.

For the time being, we also assume that $e^{ixH_{1}}$ and $e^{ixH_{2}}$ can be
easily implemented for any $x\in\mathbb{R}$. We are going to implement the
time-evolution operator

$U(x)=e^{-ixH}=e^{-ix(H_{1}+H_{2})}.$ (59)

Figure 7: Ideas behind three different methods used for LCU formula
construction of $V_{K}(x)$. ($K=1$ or $2k$, $k\in\mathbb{N}_{+}$.) (a) In the
direct Taylor-series compensation, we expand all the exponental terms in
$V_{K}(x)$ independently using Taylor-series expansion. We remove the $1$st-
to $K$th-order expansion terms based on the order condition of the Trotter
formula to reduce the $1$-norm to the order $\mathcal{O}(1+(\lambda
x)^{K+1})$. (b) In the paired Taylor-series compensation, we use the fact that
the leading-order expansion terms $F_{K,K+1}$ to $F_{K,2K+1}$ are all anti-
Hermitian. Therefore, we pair the $0$th-order term with leading order terms
using Euler’s formula to construct the Pauli rotation unitaries. This can
further suppress the $1$-norm of LCU formula to $\mathcal{O}(1+(\lambda
x)^{2K+2})$. (c) In the nested-commutator compensation, we expand the Trotter
error with the nested-commutator form. In this way, we reduce the number of
elementary unitaries in the LCU formula. For the lattice model, the number of
unitaries terms in the leading order will be linear to $n$.

We now summarized the main idea of the three methods to construct LCU formulas
based on Taylor-series compensation. We will take the lowest nontrivial
Trotter order in each case as simple examples.

  1. 1.

Direct Taylor-series compensation (Sec. E). Since $V_{K}(x)$ is a product of
exponential terms, we can expand all the exponential terms in $V_{K}(x)$
independently by Taylor-series compensation. We then expand the polynomials
and classify the terms $F_{K,s}(x)$ based on the power of $x$. Afterwards, we
consider the order condition of Trotter formulas in Lemma 1 to eliminate the
first $K$ non-identity terms in $V_{K}(x)$ with respect to the power of $x$.
This will provide a LCU formula of $V_{K}(x)$.

We take the first-order Trotter remainder as an example. The $1$st-order
Trotter formula and its remainder are

$\displaystyle S_{1}(x)$ $\displaystyle=e^{-ixH_{2}}e^{-ixH_{1}},$ (60)
$\displaystyle V_{1}(x)$
$\displaystyle=U(x)V_{1}^{\dagger}(x)=e^{-ix(H_{1}+H_{2})}e^{ixH_{1}}e^{ixH_{2}}.$

Now, we expand all the exponential terms independently in $V_{1}(x)$ by Taylor
series,

$\displaystyle V_{1}(x)$
$\displaystyle=\sum_{r=0}^{\infty}\frac{(-ix)^{r}}{r!}\left(H_{1}+H_{2}\right)^{r}\sum_{r_{1}=0}^{\infty}\frac{(ix)^{r_{1}}}{r_{1}!}H_{1}^{r_{1}}\sum_{r_{2}=0}^{\infty}\frac{(ix)^{r_{2}}}{r_{2}!}H_{2}^{r_{2}}$
(61) $\displaystyle=e^{2\lambda
x}\sum_{r;r_{1},r_{2}}\mathrm{Poi}\left(r,r_{1},r_{2};\lambda
x,\mu_{1}x,\mu_{2}x\right)(-i)^{r-r_{1}-r_{2}}\left(\sum_{l=1,2}p_{l}U_{l}\right)^{r}U_{1}^{r_{1}}U_{2}^{r_{2}},$

which is a LCU formula with the $1$-norm $e^{2\lambda x}$. Here,
$\mathrm{Poi}(r;x):=e^{-x}\frac{x^{r}}{r!}$ is the Poisson distribution. We
also denote the multiplication of Poisson distribution as,

$\displaystyle\mathrm{Poi}(r_{1},r_{2},...r_{l};x_{1},x_{2},...,x_{l})$
$\displaystyle:=\mathrm{Poi}(r_{1};x_{1})\mathrm{Poi}(r_{2};x_{2})...\mathrm{Poi}(r_{l};x_{l}).$
(62)

We now expand $V_{1}(x)$ and classify all the terms based on power of $x$,

$\displaystyle V_{1}(x)$ $\displaystyle=\sum_{s=0}^{\infty}F_{1,s}(x),$ (63)
$\displaystyle F_{1,s}(x)$
$\displaystyle=\sum_{\begin{subarray}{c}r,r_{1},r_{2}\\\
r+r_{1}+r_{2}=s\end{subarray}}\frac{(-ix)^{r}(ix)^{r_{1}}(ix)^{r_{2}}}{r!r_{1}!r_{2}!}\left(H_{1}+H_{2}\right)^{r}H_{1}^{r_{1}}H_{2}^{r_{2}}.$

The formula of $F_{1,s}(x)$ in Eq. 63 is a LCU formula with $1$-norm
$\frac{(2\lambda x)^{s}}{s!}$.

To suppress the $1$-norm of $V_{2}(x)$ in Eq. 61, we are going to utilize the
order condition of the Trotter formula in Lemma 1. We then have

$F_{1,1}(x)=0.$ (64)

The final LCU formula under direct Taylor-series compensation is then,

$\displaystyle V_{1}(x)$ $\displaystyle=I+\sum_{s=2}^{\infty}F_{1,s}(x),$ (65)

whose $1$-norm is upper bounded by $e^{(2\lambda x)^{2}}$, as proven in
Proposition 4.

The higher-order direct Taylor-series compensation follows the same logic. The
difference is that, we are able to eliminate more terms based on the order
condition of the Trotter formula. For a general $K$th-order Trotter remainder,
from Lemma 1 we have

$F_{K,1}(x)=...=F_{K,K}(x).$ (66)

This will help to further suppress the $1$-norm of the LCU formula.

  2. 2.

Paired Taylor-series compensation (Sec. F). If we carefully check the Taylor-
series expansion formula for $V_{1}(x)$ in Eq. 65, we notice that $F_{1,2}(x)$
and $F_{1,3}(x)$ are both anti-Hermitian. This is clear from the exponential
form of $V_{1}(x)$,

$V_{1}(x)=\exp(i(E_{1,2}x^{2}+E_{1,3}x^{3}+E_{1,4}x^{4}+...)),$ (67)

where $\\{E_{1,s}\\}$ are some Hermitian operators which can be obtained by
applying the Baker-Campbell-Hausdorff (BCH) formula on the definition of
$V_{1}(x)$ in Eq. 60. The first order term $E_{1,1}$ vanishes due to the order
condition in Lemma 4. We have

$\displaystyle F_{1,2}(x)$
$\displaystyle=ix^{2}E_{1,2}=ix^{2}\lambda_{1,2}\sum_{j}p_{1,2}^{(j)}P_{1,2}^{(j)},$
(68) $\displaystyle F_{1,3}(x)$
$\displaystyle=ix^{3}E_{1,3}=ix^{3}\lambda_{1,3}\sum_{j}p_{1,3}^{(j)}P_{1,3}^{(j)},$

by the Taylor-series expansion on Eq. 67. Since $E_{1,2}$ and $E_{1,3}$ are
Hermitian operators, we can naturally expand them by the Hermitian Pauli
operators $\\{P_{1,2}^{(j)}\\}$ and $\\{P_{1,3}^{(j)}\\}$ with real
coefficients. Here, we put all the minus signs in the expansion into the Pauli
operators.

The major idea of the paired Taylor-series compensation is to take the
advantage of the following Euler’s formula for Pauli operator $P$ to further
reduce the $1$-norm of the LCU formula. For a Pauli matrix $P$ and
$y\in\mathbb{R}$,

$I+iyP=\sqrt{1+y^{2}}e^{i\theta(y)P},$ (69)

where $\theta(y):=\tan^{-1}(y)$. We can see that, the $1$-norm of the formula
decrease from $1+y$ to $\sqrt{1+y^{2}}\leq 1+\frac{1}{2}y^{2}$, whose
$y$-order has been doubled.

If we apply Eq. 69 to the $0$th-order term $I$ and the $2$nd-order term
$F_{1,2}(x)$ in Eq. 68 as an example, we have

$\displaystyle I+F_{1,2}(x)$
$\displaystyle=\lambda_{1,2}\sum_{j}p_{1,2}^{j}(I+ix^{2}P_{1,2}^{(j)})$ (70)
$\displaystyle=\lambda_{1,2}\sqrt{1+x^{2}}\sum_{j}p_{1,2}^{j}\exp\left(i\theta(x^{2})P_{1,2}^{(j)}\right).$

After the expansion-order pairing, we can see that the $1$-norm becomes
$1+\mathcal{O}(x^{4})$ whose $x$-order has been doubled.

In practice, however, it is hard to calculate the explicit form of $F_{1,2}$
and $F_{1,3}$ with the Hermitian Pauli operator expansions in general.
Fortunately, in Sec. F we show that one can use the LCU formulas of $F_{1,2}$
and $F_{1,3}$ with the direct expansion form in Eq. 63. By introducing some
redundant anti-Hermitian operators, we are able to implement the LCU in the
random-sampling way with a very simple sampling scheme. We remark that, since
the anti-Hermitian operators cannot be directly paired with $I$, we
intentionally multiply them by an imaginary phase $i$. Since the summation of
all the anti-Hermitian operators in $F_{1,2}$ and $F_{1,3}$ in Eq. 63 are
zeros due to the equivalence of Eq. 63 and Eq. 68, the addtional phase $i$ on
the anti-Hermitian terms will not affect the correctness of the LCU formula.
For the higher-order Trotter with $K\geq 2$, we can perform the similar
operation to suppress the 1-norm of the LCU formula to $e^{(e+c_{k})(2\lambda
x)^{4k+2}}$, where $c_{k}<0.3$ is a small constant number defined in Eq. 202.

  3. 3.

Nested-commutator compensation (Sec. G). In the two methods above, we try to
construct the power-series expansion of $V_{K}(x)$ by first performing Taylor-
series expansion on each exponential term in it and then classify the terms by
the power of $x$. This is direct and simple, but not optimal, since many terms
in this expansion can be cancelled out to further reduce the norm.

One way to further reduce the $1$-norm of $V_{2}(x)$ is to express the $s$th-
order term $F_{2,s}(x)$ using the nested-commutators Childs and Su (2019);
Childs et al. (2021) in the following form,

$\mathrm{ad}_{H_{2}}^{m_{\kappa}}\mathrm{ad}_{H_{1}}^{n_{\kappa}}...\mathrm{ad}_{H_{2}}^{m_{2}}\mathrm{ad}_{H_{1}}^{n_{2}}\mathrm{ad}_{H_{2}}^{m_{1}}\mathrm{ad}_{H_{1}}^{n_{1}}(H_{1}+H_{2}),$
(71)

where $\sum_{j=1}^{\kappa}(m_{j}+n_{j})=s$. For two operators $P$ and $Q$, the
adjoint operator $\mathrm{ad}_{P}Q$ is defined as

$\mathrm{ad}_{P}Q:=[P,Q]=PQ-QP.$ (72)

We also define

$e^{\mathrm{ad}_{P}}Q:=\sum_{s=0}^{\infty}\frac{1}{s!}\mathrm{ad}_{P}^{s}Q=e^{P}Qe^{-P}.$
(73)

The last equality is a standard Lie-Algebraic identity.

Here, we take the $1$st-order Trotter remainder as an example. From Eq. 63 we
know that the second-order expansion term of $V_{1}(x)$ is

$\displaystyle F_{1,2}(x)$
$\displaystyle:=\sum_{\begin{subarray}{c}r,r_{1},r_{2}\\\
r+r_{1}+r_{2}=2\end{subarray}}\frac{(-ix)^{r}(ix)^{r_{1}}(ix)^{r_{2}}}{r!r_{1}!r_{2}!}\left(H_{1}+H_{2}\right)^{r}H_{1}^{r_{1}}H_{2}^{r_{2}}$
(74) $\displaystyle=\frac{(2\lambda
x)^{2}}{2}\sum_{\begin{subarray}{c}r,r_{1},r_{2}\\\
r+r_{1}+r_{2}=2\end{subarray}}\mathrm{Pr}(r,r_{1},r_{2})(-i)^{2r-2}\left(\sum_{l=1,2}p_{l}U_{l}\right)^{r}U_{1}^{r_{1}}U_{2}^{r_{2}}.$

Using Eq. 74 directly, we obtain a LCU formula of $F_{1,2}(x)$ with the
$1$-norm $(2\lambda x)^{2}/2$. However, Eq. 74 can be reduced to

$F_{1,2}(x)=-\frac{x^{2}}{2}[H_{1},H_{2}]=-\frac{x^{2}}{2}\mathrm{ad}_{H_{1}}H_{2}.$
(75)

In this case, if there are many commuted terms in $H_{1}$ and $H_{2}$, the
$1$-norm of $F_{1,2}(x)$ can be further suppressed.

In Sec. G, we derive the explicit formulas of $F_{K,s}$ ($s=K+1,...,2K$) with
the nested commutator form. We show that, for the lattice model Hamiltonian,
this will enable the Trotter-LCU algorithms with an almost-linear system size
dependence.

##  Appendix E Trotter-LCU formula with direct Taylor-series compensation

We start from the simplest LCU construction method based on direct Taylor-
series compensation. We are going to consider the Hamiltonian $H$ with Pauli
decomposition in Eq. 9, $H=\sum_{l=1}^{L}\alpha_{l}P_{l}$. The ideas behind
the direct Taylor-series compensation method are

  1. 1.

Taylor-series expansion on all the exponential terms in $V_{K}(x)$ separately.

  2. 2.

Expanding and recombining the exponential terms based on the power of $x$.

  3. 3.

Eliminate the first $K$ nonzero $x$-order terms based on the order condition
in Lemma 1.

###  E.1 Zeroth-order case

For the zeroth-order case, $S_{0}(x)=I$. The Trotter remainder $V_{0}(x)$ is
just the real-time evolution itself,

$V_{0}(x)=U(x)=e^{-ixH}.$ (76)

We then expand the remainder by Taylor series,

$\displaystyle V_{0}(x)$ $\displaystyle=\sum_{s=0}^{\infty}\frac{(-i\lambda
x)^{s}}{s!}\sum_{l_{1},...,l_{s}}p_{l_{1}}p_{l_{2}}...p_{l_{s}}P_{l_{1}}P_{l_{2}}...P_{l_{s}}$
(77) $\displaystyle=\sum_{s=0}^{\infty}F_{0,s}(x),$

where

$F_{0,s}(x)=\frac{(-i\lambda
x)^{s}}{s!}\sum_{l_{1},...,l_{s}}p_{l_{1}}p_{l_{2}}...p_{l_{s}}P_{l_{1}}P_{l_{2}}...P_{l_{s}}.$
(78)

In practice, to limit the number of gates in $V_{0,s}$ as well as avoiding
sampling over infinitely many integers, we set a truncation integer $s_{c}\geq
1$ to the expansion order $s$. After this truncation, the approximate LCU
formula of $V_{0}(x)$ is

$\tilde{V}_{0}(x)=\tilde{U}_{0}(x)=\sum_{s=0}^{s_{c}}F_{0,s}(x).$ (79)

The $1$-norms of $F_{0,s}(x)$ and $\tilde{V}_{0}(x)$ in Eq. 78 and Eq. 79 are

$\displaystyle\|F_{0,s}(x)\|_{1}$ $\displaystyle=\frac{(\lambda x)^{s}}{s!},$
(80) $\displaystyle\mu_{0}(x)$
$\displaystyle:=\|\tilde{V}_{0}(x)\|_{1}=\sum_{s=0}^{s_{c}}\|F_{0,s}(x)\|_{1}\leq
e^{\lambda x}$

We can write $F_{0,s}(x)$ and $\tilde{V}_{0}(x)$ in the normalized form,

$\displaystyle F_{0,s}(x)$ $\displaystyle=\|F_{0,s}(x)\|_{1}V_{0,s},$ (81)
$\displaystyle V_{0,s}$
$\displaystyle:=(-i)^{s}\sum_{l_{1},...,l_{s}}p_{l_{1}}p_{l_{2}}...p_{l_{s}}P_{l_{1}}P_{l_{2}}...P_{l_{s}},$
$\displaystyle\tilde{V}_{0}(x)$
$\displaystyle=\mu_{0}(x)\sum_{s=0}^{s_{c}}\mathrm{Pr}_{0}(s)V_{0,s},$
$\displaystyle\mathrm{Pr}_{0}(s)$
$\displaystyle:=\frac{\|F_{0,s}(x)\|_{1}}{\mu_{0}(x)},s=0,1,...,s_{c}.$

We have the following proposition to characterize the LCU formula in Eq. 79.

######  Proposition 3 (0th-order Trotter-LCU formula by direct Taylor-series
compensation).

For $x\geq 0$ and $s_{c}\geq 1$, $\tilde{V}_{0}(x)$ in Eq. 79 is a
$(\mu_{0}(x),\varepsilon_{0}(x))$-LCU formula of $V_{0}(x)=U(x)$ with

$\displaystyle\mu_{0}(x)$ $\displaystyle\leq e^{\lambda x},$ (82)
$\displaystyle\varepsilon_{0}(x)$ $\displaystyle\leq\left(\frac{e\lambda
x}{s_{c}+1}\right)^{s_{c}+1}.$

###### Proof.

Since the $1$-norm of $\tilde{V}_{0}(x)$ in Eq. 79 is the same as $V_{0,s}$ in
Eq. 81, we naturally have $\mu_{0}(x)=e^{\lambda x}$.

For the distance bound, we have

$\displaystyle\|\tilde{V}_{0}(x)-V_{0}(x)\|$
$\displaystyle\leq\sum_{s>s_{c}}\|F_{0,s}(x)\|$ (83)
$\displaystyle\leq\sum_{s>s_{c}}\frac{(\lambda x)^{s}}{s!}\|V_{0,s}(x)\|$
$\displaystyle\leq\sum_{s>s_{c}}\frac{(\lambda x)^{s}}{s!}$
$\displaystyle\leq\left(\frac{e\lambda x}{s_{c}+1}\right)^{s_{c}+1}.$

In the second inequality, we use the fact that

$\|V_{0,s}\|\leq\sum_{l_{1},l_{2},...,l_{s}}p_{l_{1}}p_{l_{2}}...p_{l_{s}}\|(-i)^{s}P_{l_{1}}P_{l_{2}}...P_{l_{s}}\|\leq
1.$ (84)

In the forth line of Eq. 83, we apply Corollary 2 in Appendix K. ∎

Now, we consider the Trotter-LCU formula $\tilde{U}_{0}(x)$ in Eq. 79 for both
random-sampling and coherent usages.

For the random-sampling implementation, the overall LCU formula for $U(t)$ is
to repeat $\tilde{V}_{0}(x)=\tilde{U}_{0}(x)$ for $\nu$ times,
$\tilde{U}_{0,tot}(t)=\tilde{U}_{0}(x)^{\nu}$. Using Proposition 2 and 3, we
conclude that $\tilde{U}_{0,tot}(t)$ is a
$(\mu_{0,tot}(t),\varepsilon_{0,tot}(t))$-LCU formula of $U(t)$ with

$\displaystyle\mu_{0,tot}(t)$ $\displaystyle=\mu_{0}(x)^{\nu}\leq e^{\lambda
t},$ (85) $\displaystyle\varepsilon_{0,tot}(t)$
$\displaystyle\leq\nu\mu_{0,tot}(t)\varepsilon_{0}(x)\leq\nu e^{\lambda
t}\left(\frac{e\lambda x}{s_{c}+1}\right)^{s_{c}+1}.$

In this case, the $1$-norm will exponentially increasing no matter what
segment number $\nu$ we choice. Consequently, the $0$th-order Trotter-LCU
formula with direct Taylor-series compensation is not suitable for the random-
sampling usage.

The inefficiency of $0$th-order random-sampling implement can also be
understood from another perspective of classical simulation. Note that, all
the gates in the $0$th-order algorithm are Clifford gates — to implement
controlled-$\tilde{U}_{0}(x)$, we only need to perform many controlled-Pauli
gates. Based on the Gottesman-Knill’s theorem Aaronson and Gottesman (2004),
we know that this circuit can be easily simulated by classical computers. On
the other hand, it is widely believed that the real-time evolution of a
generic Hamiltonian is exponentially hard to simulate on a classical computer
Schuch et al. (2008). As a result, this will let to the exponentially
increasing sampling complexity to the evolution time $t$.

For the coherent implementation, to apply Corollary 1, we need to set

$\mu_{0}(x)=2\Leftarrow x=\frac{\ln 2}{\lambda}.$ (86)

Therefore, the number of segments $\nu$ is then

$\nu=\frac{t}{x}=\frac{\lambda t}{\ln 2}=\mathcal{O}(\lambda t).$ (87)

To determine the truncation order $s_{c}$ based on the overall simulation
accuracy $\varepsilon$, we only need to make sure that the accuracy for each
segment is $\xi:=\varepsilon/\nu$. Then we can solve $s_{c}$ based on the
following constraints

$\displaystyle\epsilon_{0}:=\left(\frac{e\ln 2}{s_{c}+1}\right)^{s_{c}+1},$
(88)
$\displaystyle\frac{\epsilon_{0}^{2}+2\epsilon_{0}+4}{2}\epsilon_{0}\leq\frac{\varepsilon}{\nu}.$

When $\epsilon_{0}$ is small enough such that
$\epsilon_{0}(\epsilon_{0}+2)\leq 2$, it suffices to have

$\displaystyle 3\nu\epsilon_{0}=3\nu\left(\frac{e\ln
2}{s_{c}+1}\right)^{s_{c}+1}\leq\varepsilon$ (89) $\displaystyle\Rightarrow$
$\displaystyle
s_{c}\geq\left\lceil\frac{\ln(\frac{\varepsilon}{3\nu})}{W_{0}\left(\frac{1}{e\ln
2}\ln(\frac{\varepsilon}{3\nu})\right)}-1\right\rceil.$

Here, the second line is derived based on Lemma 8. $W_{0}(y)$ is the principle
branch of the Lambert-W function. From Lemma 9 we know that
$W_{0}(y)=\mathcal{O}(\log(y))$. As a result,
$s_{c}=\mathcal{O}(\log(1/\varepsilon)/\log\log(1/\varepsilon))$.

###  E.2 First-order case

Now, we study the direct Taylor-series compensation of the remainders of
Trotter formulas. We will take $1$st-order case as an example. The Trotter
remainder $V_{1}(x)$ is

$V_{1}(x)=U(x)S_{1}(x)^{\dagger}.$ (90)

From Eq. 77 and Eq. 16 we have

$\displaystyle U(x)$ $\displaystyle=\sum_{r=0}^{\infty}\frac{(\lambda
x)^{r}}{r!}\sum_{l_{1},...,l_{r}}p_{l_{1}}p_{l_{2}}...p_{l_{r}}(-i)^{r}P_{l_{1}}P_{l_{2}}...P_{l_{r}}$
(91) $\displaystyle=e^{\lambda x}\sum_{r=0}^{\infty}\mathrm{Poi}(r;\lambda
x)\sum_{l_{1:r}}p^{(r)}_{l_{1:r}}(-i)^{r}P^{(r)}_{l_{1:r}}$ $\displaystyle
S_{1}(x)^{\dagger}$
$\displaystyle=\prod^{\rightarrow}e^{ix\vec{H}}=e^{ixH_{1}}e^{ixH_{2}}...e^{ixH_{L}}$
$\displaystyle=\left(\prod
e^{\vec{\alpha}x}\right)\sum_{\vec{r}}\mathrm{Poi}(\vec{r};\vec{\alpha}x)\prod^{\rightarrow}(i\vec{P})^{\vec{r}}.$

Here, we adopt the vector notations introduced in Sec. D to simplify the
expressions. In Eq. 91, we also extend the notation of Poisson distribution,

$\mathrm{Poi}(\vec{r};\vec{\alpha}x):=\prod_{l=1}^{L}\mathrm{Poi}(r_{l};\alpha_{l}x).$
(92)

Based on Eq. 91, we then write the Taylor-series expansion of first-order
remainder as follows,

$\displaystyle V_{1}(x)$ $\displaystyle=U(x)S_{1}(x)^{\dagger}$ (93)
$\displaystyle=\sum_{r;\vec{r}}\frac{(\lambda
x)^{r}}{r!}\left(\prod\frac{(\vec{\alpha}x)^{\vec{r}}}{\vec{r}!}\right)\sum_{l_{1:r}}p^{(r)}_{l_{1:r}}(-i)^{r-\sum\vec{r}}P^{(r)}_{l_{1:r}}\prod^{\rightarrow}\vec{P}^{\vec{r}},$
$\displaystyle=e^{2\lambda x}\sum_{r;\vec{r}}\mathrm{Poi}(r,\vec{r};\lambda
x,\vec{\alpha}x)\sum_{l_{1:r}}p^{(r)}_{l_{1:r}}(-i)^{r-\sum\vec{r}}P^{(r)}_{l_{1:r}}\prod^{\rightarrow}\vec{P}^{\vec{r}},$

which is a LCU formula with $1$-norm $e^{2\lambda x}$. Here, $r$ and $\vec{r}$
are two groups of independent variables.

Now, we use the order condition of the Trotter formula in Lemma 1 to eliminate
some terms in Eq. 93, so that the $1$-norm is suppressed. To this end, we
first rewrite Eq. 93 by classifying the terms based on the power of $x$, which
is determined by the value $s=r+\sum\vec{r}$.

$\displaystyle
V_{1}(x)=\sum_{s=0}^{\infty}F_{1,s}(x)=\sum_{s=0}^{\infty}\|F_{1,s}(x)\|_{1}V_{1,s}.$
(94)

Here, $\|F_{1,s}(x)\|_{1}$ is the $1$-norm of the $s$-order expansion formula
$F_{1,s}(x)$. $V_{1,s}$ denotes the normalized LCU formula for the $s$-order
terms. We have

$\displaystyle F_{1,s}(x)$
$\displaystyle:=\sum_{r+\sum\vec{r}=s}\frac{(\lambda
x)^{r}}{r!}\left(\prod\frac{(\vec{\alpha}x)^{\vec{r}}}{\vec{r}!}\right)\sum_{l_{1:r}}p^{(r)}_{l_{1:r}}(-i)^{2r-s}P^{(r)}_{l_{1:r}}\prod^{\rightarrow}\vec{P}^{\vec{r}},$
(95) $\displaystyle V_{1,s}$
$\displaystyle:=\sum_{r;\vec{r}}\Pr(r,\vec{r}|s)\sum_{l_{1:r}}p^{(r)}_{l_{1:r}}(-i)^{2r-s}P^{(r)}_{l_{1:r}}\prod^{\rightarrow}\vec{P}^{\vec{r}}.$

In the expression of $V_{1}^{(s)}$, we use $r-\sum\vec{r}=r-(s-r)=2r-s$. The
conditional probability $\Pr(r,\vec{r}|s)$ indicates the probability to sample
$r$ and $\vec{r}$ when their summation $s$ is given,

$\displaystyle\Pr(r,\vec{r}|s)$
$\displaystyle=\frac{\Pr\\{(\hat{r}=r,\hat{\vec{r}}=\vec{r})\cap(\hat{s}=s)\\}}{\Pr\\{\hat{s}=s\\}}$
(96)
$\displaystyle=\frac{s!}{r!\prod\vec{r}!}\left(\frac{1}{2}\right)^{r}\prod\left(\frac{\vec{p}}{2}\right)^{\vec{r}}$
$\displaystyle=\mathrm{Mul}(\\{r,\vec{r}\\};\\{\frac{1}{2},\frac{\vec{p}}{2}\\};s),$

which is a multinomial distribution $\mathrm{Mul}(\cdot;\cdot;s)$ with $s$
trials and $(L+1)$-outcomes. In each trial, the $(L+1)$ outcomes
$\\{r;r_{1},r_{2},...,r_{L}\\}$ occur with the corresponding probability
$\\{\frac{1}{2},\frac{1}{2}p_{1},\frac{1}{2}p_{2},...,\frac{1}{2}p_{L}\\}$.
Recall that $\vec{p}$ is defined in Eq. 10.

We first estimate the normalization cost $\|F_{1,s}(x)\|_{1}$ for different
$s$ orders,

$\displaystyle\|F_{1,s}(x)\|_{1}$
$\displaystyle=\sum_{s;\vec{r}}\mathbbm{1}\left[r+\sum\vec{r}=s\right]\frac{(\lambda
x)^{r}}{r!}\left(\prod\frac{(\vec{\alpha}x)^{\vec{r}}}{\vec{r}!}\right)$ (97)
$\displaystyle=\frac{(\lambda x+\sum\vec{\alpha}x)}{s!}=\frac{(2\lambda
x)^{s}}{s!}:=\eta_{s}.$

In the second inequality, we use Corollary 4. We denote

$\eta_{s}=\eta_{s}(2\lambda x):=\frac{(2\lambda x)^{s}}{s!}.$ (98)

We will frequently use $\eta_{s}$ in the following discussion.

From Eq. 97 we can see that, similar to the expansion of $V_{0}(x)$, the
$1$-norm of the expansion terms with different orders of $x$ in $V_{1}(x)$
follow the Possion distribution. This motivates us to eliminate the low-order
terms such as $F_{1,1}(x)$. Based on Lemma 1, we have

$F_{1}^{(1)}(x)=0.$ (99)

As a result, we can directly remove the term $F_{1}^{(1)}(x)$ in Eq. 94. The
resulting formula is,

$\displaystyle V_{1}(x)$ $\displaystyle=I+\sum_{s=2}^{\infty}\eta_{s}V_{1,s}.$
(100)

In practice, we introduce a truncation $s_{c}\geq 1$ in the expansion order of
$x$. The truncated formula is

$\displaystyle\tilde{V}_{1}(x)$
$\displaystyle=I+\sum_{s=2}^{s_{c}}\eta_{s}V_{1,s}$ (101)
$\displaystyle=\mu_{1}(x)\left(\mathrm{Pr}_{1}(0)I+\sum_{s=2}^{s_{c}}\mathrm{Pr}_{1}(s)V_{1,s}\right),$

where

$\displaystyle\mu_{1}(x)$
$\displaystyle:=\eta_{0}+\sum_{s=2}^{s_{c}}\eta_{s}\leq e^{2\lambda
x}-(2\lambda x)\leq e^{(2\lambda x)^{2}},$ (102)
$\displaystyle\mathrm{Pr}_{1}(s)$
$\displaystyle:=\frac{\eta_{s}}{\mu_{1}(x)}=\frac{1}{\mu_{1}(x)}\frac{(2\lambda
x)^{s}}{s!},s=0;2,3,...,s_{c}.$ $\displaystyle V_{1,s}$
$\displaystyle=\sum_{r;\vec{r}}\mathrm{Pr}(r,\vec{r}|s)\sum_{l_{1:r}}p^{(r)}_{l_{1:r}}(-i)^{2r-s}P^{(r)}_{l_{1:r}}\prod^{\rightarrow}\vec{P}^{\vec{r}}.$

Here, the conditional probability
$\mathrm{Pr}(r,\vec{r}|s):=\mathrm{Mul}(\\{r,\vec{r}\\};\\{\frac{1}{2},\frac{\vec{p}}{2}\\};s)$.

Combined with the deterministic first-order Trotter formula, the overall LCU
formula for $U(x)$ is

$\displaystyle\tilde{U}_{1}(x)=\tilde{V}_{1}(x)S_{1}(x).$ (103)

We have the following proposition to characterize the LCU formula in Eq. 101
and Eq. 103.

######  Proposition 4 (1st-order Trotter-LCU formula by direct Taylor-series
compensation).

For $x>0$ and $s_{c}\geq 1$, $\tilde{V}_{1}(x)$ in Eq. 101 is a
$(\mu_{1}(x),\varepsilon_{1}(x))$-LCU formula of $V_{1}(x)$ with

$\displaystyle\mu_{1}(x)$ $\displaystyle\leq e^{(2\lambda x)^{2}},$ (104)
$\displaystyle\varepsilon_{1}(x)$ $\displaystyle\leq\left(\frac{2e\lambda
x}{s_{c}+1}\right)^{s_{c}+1}.$

As a result, $\tilde{U}_{1}(x)$ in Eq. 103 is a
$(\mu_{1}(x),\varepsilon_{1}(x))$-LCU formula of $U(x)$.

###### Proof.

From $e^{x}-x\leq e^{x^{2}}$ for $x>0$, we know that

$\mu_{1}(x)\leq e^{2\lambda x}-2\lambda x\leq e^{(2\lambda x)^{2}}.$ (105)

For the distance bound, from Eq. 100 and Eq. 101 we have

$\displaystyle\|V_{1}(x)-\tilde{V}_{1}(x)\|$
$\displaystyle\leq\sum_{s>s_{c}}\eta_{s}\|V_{1,s}\|$
$\displaystyle=\sum_{s>s_{c}}\frac{(2\lambda x)^{s}}{s!}$ (106)
$\displaystyle\leq\left(\frac{2e\lambda x}{s_{c}+1}\right)^{s_{c}+1}$

In the second line, we use the fact that $\|V_{1,s}\|\leq 1$. In the third
line, we apply Corollary 2. ∎

Now, we consider the Trotter-LCU formula $\tilde{U}_{1}(x)$ in Eq. 103 for
both random-sampling and coherent usages.

For the random-sampling implementation, the overall LCU formula for $U(t)$ is
to repeat $\tilde{U}_{1}(x)$ for $\nu$ times,
$\tilde{U}_{1,tot}(t)=\tilde{U}_{1}(x)^{\nu}$. Using Proposition 2 and 4, we
conclude that $\tilde{U}_{1}(t)$ is a
$(\mu_{1,tot}(t),\varepsilon_{1,tot}(t))$-LCU formula of $U(t)$ with

$\displaystyle\mu_{1,tot}(t)$ $\displaystyle=\mu_{1}(x)^{\nu}\leq
e^{\frac{(2\lambda t)^{2}}{\nu}},$ (107) $\displaystyle\varepsilon_{1,tot}(t)$
$\displaystyle\leq\nu\mu_{1,tot}(t)\varepsilon_{1}(x)\leq\nu
e^{\frac{(2\lambda t)^{2}}{\nu}}\left(\frac{2e\lambda
x}{s_{c}+1}\right)^{s_{c}+1}.$

To realize a $(\mu,\varepsilon)$-LCU formula for $U(t)$, we only need to set
the segment number $\nu$ and the truncation order $s_{c}$ to be

$\displaystyle\nu$ $\displaystyle\geq\frac{(2\lambda t)^{2}}{\ln\mu},$ (108)
$\displaystyle s_{c}$
$\displaystyle\geq\left\lceil\frac{\ln\left(\frac{\mu(2\lambda
t)^{2}}{\ln(\mu)\varepsilon}\right)}{W_{0}\left(\frac{2\lambda
t}{e\ln(\mu)}\ln\left(\frac{\mu(2\lambda
t)^{2}}{\ln(\mu)\varepsilon}\right)\right)}-1\right\rceil,$

To derive the bound for $s_{c}$ in the second line, we use Lemma 8. Based on
Eq. 41, the gate complexity of the $1$st-order Trotter-LCU algorithm is then

$\mathcal{O}(\nu(L+s_{c}))=\mathcal{O}\left((\lambda
t)^{2}(L+\frac{\log(1/\varepsilon)}{\log\log(1/\varepsilon)})\right).$ (109)

From Eq. 108 we can see that, since we eliminate the first-order term in
$V_{1}(x)$ based on the order condition in Lemma 1, we suppress the 1-norm.
This enables us to consider an efficient random-sampling implementation of
Trotter-LCU formula.

For the coherent implementation of 1st-order DTSc Trotter-LCU, to apply
Corollary 1, we need to set

$\mu_{1}(x)\leq 2\Leftarrow x=\frac{\sqrt{\ln 2}}{2\lambda}.$ (110)

Therefore, the number of segments $\nu$ is then

$\nu=\frac{t}{x}=\frac{2\lambda t}{\sqrt{\ln 2}}=\mathcal{O}(\lambda t).$
(111)

As a result, the smaller value of $\mu(x)$ for each segment is not helpful to
reduce the number of segments $\nu$ for the coherent implementation of LCU.
The truncation order analysis for $s_{c}$ is similar to the one in Sec. E.1.
We also have $s_{c}=\mathcal{O}(\log(1/\varepsilon)/\log\log(1/\varepsilon))$.

###  E.3 Second- and higher-order case

Following similar methods, we can obtain the LCU formulas for $2$nd- and
higher-order Trotter remainders.

For example, the $2$nd-order Trotter remainder $V_{2}(x)$ is

$\displaystyle V_{2}(x)$ $\displaystyle=U(x)S_{2}(x)^{\dagger}=U(x)$ (112)
$\displaystyle=\sum_{r,\vec{r},\vec{r}^{\prime}}\frac{(\lambda
x)^{r}}{r!}\prod\frac{\left(\vec{\alpha}\frac{x}{2}\right)^{\vec{r}}}{\vec{r}!}\prod\frac{\left(\vec{\alpha}\frac{x}{2}\right)^{\vec{r}^{\prime}}}{\vec{r}^{\prime}!}\sum_{l_{1:r}}p_{l_{1:r}}^{(r)}(-i)^{r-\sum(\vec{r}+\vec{r}^{\prime})}P^{(r)}_{l_{1:r}}\prod^{\rightarrow}\vec{P}^{\vec{r}}\prod^{\leftarrow}\vec{P}^{\vec{r}^{\prime}}.$

We then expand $V_{2}(x)$ based on the power of $x$, given by
$s:=r+\sum(\vec{r}+\vec{r}^{\prime})$,

$\displaystyle
V_{2}(x)=\sum_{s=0}^{\infty}F_{2,s}(x)=\sum_{s=0}^{\infty}\|F_{2,s}(x)\|_{1}V_{2,s},$
(113)

where $\|F_{2,s}(x)\|_{1}$ is the 1-norm of the $s$-order expansion formula
$F_{2,s}(x)$. $V_{2,s}$ denotes the normalized LCU formula for the $s$-order
terms. Similar to the derivations in Eq. 96 and Eq. 97, we have

$\displaystyle\|F_{2,s}(x)\|_{1}$
$\displaystyle=\sum_{\begin{subarray}{c}r,\vec{r},\vec{r}^{\prime}\\\
r+\sum(\vec{r}+\vec{r}^{\prime})=s\end{subarray}}\frac{(\lambda
x)^{r}}{r!}\prod\frac{\left(\vec{\alpha}\frac{x}{2}\right)^{\vec{r}}}{\vec{r}!}\prod\frac{\left(\vec{\alpha}\frac{x}{2}\right)^{\vec{r}^{\prime}}}{\vec{r}^{\prime}!}=\frac{(2\lambda
x)^{s}}{s!}:=\eta_{s},$ (114) $\displaystyle V_{2,s}$
$\displaystyle=\sum_{r,\vec{r},\vec{r}^{\prime}}\Pr(r,\vec{r},\vec{r}^{\prime}|s)\sum_{l_{1:r}}p_{l_{1:r}}^{(r)}(-i)^{2r-s}P^{(r)}_{l_{1:r}}\prod^{\rightarrow}\vec{P}^{\vec{r}}\prod^{\leftarrow}\vec{P}^{\vec{r}^{\prime}}.$

Here, the conditional probability $\Pr(r,\vec{r},\vec{r}^{\prime}|s)$ is the
multinomial distribution
$\mathrm{Mul}(\\{r,\vec{r},\vec{r}^{\prime}\\};\\{\frac{1}{2},\frac{\vec{p}}{4},\frac{\vec{p}}{4}\\};s)$
following similar derivation in Eq. 96. In the expression of $V_{2,s}$, we use
$r-\sum(\vec{r}+\vec{r}^{\prime})=r-(s-r)=2r-s$.

Based on Lemma 1 we know that

$F_{2,1}(x)=F_{2,2}(x)=0.$ (115)

Therefore, we can reduce the expansion formula of $V_{2}(x)$ to

$\displaystyle V_{2}(x)=I+\sum_{s=3}^{\infty}\eta_{s}V_{2,s}.$ (116)

In practice, we introduce a truncation $s_{c}\geq 3$ in the expansion order.
We remark that, if $s_{c}<3$, we do not need to introduce the LCU formula to
compensate the Trotter error. The truncated formula is

$\displaystyle\tilde{V}_{2}(x)$
$\displaystyle=I+\sum_{s=3}^{s_{c}}\eta_{s}V_{2,s}$ (117)
$\displaystyle=\mu_{2}(x)\left(\mathrm{Pr}_{2}(0)I+\sum_{s=3}^{s_{c}}\mathrm{Pr}_{2}(s)V_{2,s}\right),$

where

$\displaystyle\mu_{2}(x)$
$\displaystyle:=\eta_{0}+\sum_{s=3}^{s_{c}}\eta_{s},$ (118)
$\displaystyle\mathrm{Pr}_{2}(s)$
$\displaystyle:=\frac{\eta_{s}}{\mu_{2}(x)}=\frac{1}{\mu_{2}(x)}\frac{(2\lambda
x)^{s}}{s!},s=0;3,4,...,s_{c}.$

Combined with the deterministic second-order Trotter formula, the overall LCU
formula for $U(x)$ is

$\displaystyle\tilde{U}_{2}(x)=\tilde{V}_{2}(x)S_{2}(x).$ (119)

The derivation of LCU formula for the $2k$th-order Trotter remainder
$V_{2k}(x)$ is similar to the second order case. From Eq. 22, the remainder of
the general $2k$th-order Trotter formula is,

$V_{2k}(x)=U(x)S_{2k}(x)^{\dagger}=U(x)D_{2k}(x)^{\dagger}D_{2k}(-x),$ (120)

which can be expanded as

$\displaystyle V_{2k}(x)$
$\displaystyle=\sum_{\begin{subarray}{c}r;\vec{r}_{1},...,\vec{r}_{\kappa};\\\
\vec{r}_{1}^{\prime},...,\vec{r}_{\kappa}^{\prime}\end{subarray}}\frac{(\lambda
x)^{r}}{r!}\prod\frac{\left(\vec{\alpha}\nu_{1}x\right)^{\vec{r}_{1}}}{\vec{r}_{1}!}...\prod\frac{\left(\vec{\alpha}\nu_{\kappa}x\right)^{\vec{r}_{\kappa}}}{\vec{r}_{\kappa}!}\,\prod\frac{\left(\vec{\alpha}\nu_{\kappa}x\right)^{\vec{r}^{\prime}_{\kappa}}}{\vec{r}^{\prime}_{\kappa}!}...\prod\frac{\left(\vec{\alpha}\nu_{1}x\right)^{\vec{r}^{\prime}_{1}}}{\vec{r}^{\prime}_{1}!}$
(121)
$\displaystyle\quad\quad\sum_{l_{1:r}}p_{l_{1:r}}^{(r)}(-i)^{r-\sum(\vec{r}_{1}+...+\vec{r}_{\kappa}+\vec{r}_{1}^{\prime}+...+\vec{r}_{\kappa}^{\prime})}\left(P^{(r)}_{l_{1:r}}\right)\prod^{\rightarrow}\vec{P}^{\vec{r}_{1}}\prod^{\leftarrow}\vec{P}^{\vec{r}_{2}}...\prod^{\rightarrow}\vec{P}^{\vec{r}_{\kappa}}\prod^{\leftarrow}\vec{P}^{\vec{r}^{\prime}_{\kappa}}...\prod^{\rightarrow}\vec{P}^{\vec{r}^{\prime}_{2}}\prod^{\leftarrow}\vec{P}^{\vec{r}^{\prime}_{1}}.$

Here, we use $\vec{r}_{i}=[r_{i,1},r_{i,2},...,r_{i,L}]$ to denote a vector
with $L$ elements. In what follows, we denote

$\vec{r}_{1:\kappa}=\\{\vec{r}_{1},\vec{r}_{2},...,\vec{r}_{\kappa}\\}.$ (122)

As a result, there are $\kappa L$ summands $r_{i,j}$ in $\vec{r}_{1:\kappa}$.
Recall that $\kappa=5^{k-1}$, $\\{v_{i}\\}_{i=1}^{\kappa}$ are positive
numbers indicating the segment length, which are determined by Eq. 18 and Eq.
19 which satisfy $\sum_{i=1}^{\kappa}v_{i}=\frac{1}{2}$.

Expanding $V_{2k}(x)$ based on the power of $x$, given by
$s:=r+\sum(\vec{r}_{1:\kappa}+\vec{r}^{\prime}_{1:\kappa})$, we have

$\displaystyle V_{2k}(x)$
$\displaystyle=\sum_{s=0}^{\infty}F_{2k,s}(x)=\sum_{s=0}^{\infty}\|F_{2k,s}(x)\|_{1}V_{2k,s},$
(123) $\displaystyle\|F_{2k,s}(x)\|_{1}$
$\displaystyle=\sum_{r;\vec{r}_{1:\kappa};\vec{r}^{\prime}_{1:\kappa}}\frac{(\lambda
x)^{r}}{r!}\prod\frac{\left(\vec{\alpha}\nu_{1}x\right)^{\vec{r}_{1}}}{\vec{r}_{1}!}...\prod\frac{\left(\vec{\alpha}\nu_{\kappa}x\right)^{\vec{r}_{\kappa}}}{\vec{r}_{\kappa}!}\prod\frac{\left(\vec{\alpha}\nu_{\kappa}x\right)^{\vec{r}^{\prime}_{\kappa}}}{\vec{r}^{\prime}_{\kappa}!}...\prod\frac{\left(\vec{\alpha}\nu_{1}x\right)^{\vec{r}^{\prime}_{1}}}{\vec{r}^{\prime}_{1}!}=\eta_{s},$
$\displaystyle V_{2k,s}$
$\displaystyle=\sum_{r,\vec{r}_{1:\kappa},\vec{r}^{\prime}_{1:\kappa}}\Pr(r,\vec{r}_{1:\kappa},\vec{r}^{\prime}_{1:\kappa}|s)\sum_{l_{1:r}}p_{l_{1:r}}^{(r)}(-i)^{2r-s}\left(P^{(r)}_{l_{1:r}}\right)\prod^{\rightarrow}\vec{P}^{\vec{r}_{1}}\prod^{\leftarrow}\vec{P}^{\vec{r}_{2}}...\prod^{\rightarrow}\vec{P}^{\vec{r}_{\kappa}}\prod^{\leftarrow}\vec{P}^{\vec{r}^{\prime}_{\kappa}}...\prod^{\rightarrow}\vec{P}^{\vec{r}^{\prime}_{2}}\prod^{\leftarrow}\vec{P}^{\vec{r}^{\prime}_{1}}.$

Where
$\Pr(r,\vec{r}_{1:\kappa},\vec{r}^{\prime}_{1:\kappa}|s):=\mathrm{Mul}(\\{r,\vec{r}_{1:\kappa},\vec{r}^{\prime}_{1:\kappa}\\};\\{\frac{1}{2},\nu_{1:\kappa}\frac{\vec{p}}{2},\nu_{1:\kappa}\frac{\vec{p}}{2}\\};s)$.
In the expression of $V_{2k,s}$, we use
$r-\sum(\vec{r}_{1:\kappa}+\vec{r}^{\prime}_{1:\kappa})=r-(s-r)=2r-s$. We use
the notation
$\nu_{1:k}\vec{p}=\\{\nu_{1}\vec{p},\nu_{2}\vec{p},...,\nu_{\kappa}\vec{p}\\}$.

Based on Lemma 1 we know that

$F_{2k}^{(1)}(x)=F_{2k}^{(2)}(x)=...=F_{2k}^{(2k)}=0.$ (124)

Therefore, we can reduce the expansion formula of $V_{2k}(x)$ to

$\displaystyle V_{2k}(x)=I+\sum_{s=2k+1}^{\infty}\eta_{s}V_{2k,s}.$ (125)

In practice, we introduce a truncation $s_{c}\geq 2k+1$ in the expansion
order. The truncated formula is then

$\displaystyle\tilde{V}_{2k}(x)$
$\displaystyle=I+\sum_{s=2k+1}^{s_{c}}\eta_{s}V_{2k,s}$ (126)
$\displaystyle=\mu_{2k}(x)\left(\mathrm{Pr}_{2k}(s=0)I+\sum_{s=2k+1}^{s_{c}}\mathrm{Pr}_{2k}(s)V_{2k,s}\right).$

where

$\displaystyle\mu_{2k}(x)$
$\displaystyle:=\eta_{0}+\sum_{s=2k+1}^{s_{c}}\eta_{s}V_{2k,s}\leq e^{2\lambda
x}-\sum_{s=1}^{2k}\frac{(2\lambda x)^{s}}{s!},$ (127)
$\displaystyle\mathrm{Pr}_{2k}(s)$
$\displaystyle:=\frac{\eta_{s}}{\mu_{2k}(x)}=\frac{1}{\mu_{2k}(x)}\frac{(2\lambda
x)^{s}}{s!}.$

Combined with the deterministic second-order Trotter formula, the overall LCU
formula for $U(x)$ is

$\displaystyle\tilde{U}_{2k}(x)=\tilde{V}_{2k}(x)S_{2k}(x).$ (128)

We have the following proposition to characterize the LCU formula for
$V_{2k}(x)$ and $U_{2k}(x)$ in Eq. 126 and Eq. 128.

######  Proposition 5 ($2k$th-order Trotter-LCU formula by direct Taylor-
series compensation).

For $0<x<(2k+1)/(2\lambda)$ and $s_{c}\geq 2k+1$, $\tilde{V}_{2k}(x)$ in Eq.
126 is a $(\mu_{2k}(x),\varepsilon_{2k}(x))$-LCU formula of $V_{2k}(x)$ with

$\displaystyle\mu_{2k}(x)$ $\displaystyle\leq e^{e(2\lambda x)^{2k+1}},$ (129)
$\displaystyle\varepsilon_{2k}(x)$ $\displaystyle\leq\left(\frac{2e\lambda
x}{s_{c}+1}\right)^{s_{c}+1}.$

As a result, $\tilde{U}_{2k}(x)$ in Eq. 128 is a
$(\mu_{2k}(x),\varepsilon_{2k}(x))$-LCU formula of $U(x)$.

###### Proof.

Here, we focus on the case with $k=1$, i.e., the second-order case. The proof
for general $k$ can be done in the same way.

From Eq. 118 and Corollary 3, we know that

$\mu_{2}(x)\leq e^{2\lambda x}-\sum_{s=1}^{2}\frac{(2\lambda x)^{s}}{s!}\leq
e^{e(2\lambda x)^{3}}.$ (130)

For the distance bound, from Eq. 116 and Eq. 117 we have

$\displaystyle\|V_{2}(x)-\tilde{V}_{2}(x)\|$
$\displaystyle\leq\sum_{s>s_{c}}\eta_{s}\|V_{2,s}\|$
$\displaystyle=\sum_{s>s_{c}}\frac{(2\lambda x)^{s}}{s!}$ (131)
$\displaystyle\leq\left(\frac{2e\lambda x}{s_{c}+1}\right)^{s_{c}+1}$

In the second line, we use the fact that $\|V_{2,s}\|\leq 1$. In the forth
line, we apply Corollary 2.

∎

Now, we consider the Trotter-LCU formula $\tilde{U}_{2k}(x)$ in Eq. 128 for
both random-sampling and coherent usages.

For the random-sampling implementation, the overall LCU formula for $U(t)$ is
to repeat the sampling of $\tilde{U}_{2k}(x)$ for $\nu$ times,
$\tilde{U}_{2k,tot}(t)=\tilde{U}_{2k}(x)^{\nu}$. Using Proposition 2 and 5,
when $0<x<\frac{2k+1}{2\lambda}$ and $s_{c}\geq 1$, we conclude that
$\tilde{U}_{2k,tot}(t)$ is a $(\mu_{2k,tot}(t),\epsilon_{2k,tot}(t))$-LCU
formula of $U(t)$ with

$\displaystyle\mu_{2k,tot}(t)$ $\displaystyle=\mu_{2k}(x)^{\nu}\leq
e^{e\frac{(2\lambda t)^{2k+1}}{\nu^{2k}}},$ (132)
$\displaystyle\varepsilon_{2k,tot}(t)$
$\displaystyle\leq\nu\mu_{2k,tot}(t)\varepsilon_{2k}(x)\leq\nu
e^{e\frac{(2\lambda t)^{2k+1}}{\nu^{2k}}}\left(\frac{2e\lambda
x}{s_{c}+1}\right)^{s_{c}+1}.$

To realize a $(\mu,\varepsilon)$-LCU formula for $U(t)$, we only need to set
the segment number $\nu$ and the truncation order $s_{c}$ to be

$\displaystyle\nu$ $\displaystyle\geq\max\left\\{\nu_{2k}(t),\frac{2\lambda
t}{2k+1}\right\\},$ (133) $\displaystyle s_{c}$
$\displaystyle\geq\max\left\\{\left\lceil\frac{\ln\left(\frac{\mu}{\varepsilon}\nu_{2k}(t)\right)}{W_{0}\left(\frac{1}{2e\lambda
t}\nu_{2k}(t)\ln\left(\frac{\mu}{\varepsilon}\nu_{2k}(t)\right)\right)}-1\right\rceil,2k+1\right\\},$

where $\nu_{2k}(t):=\left(\frac{2e\lambda
t}{\ln\mu}\right)^{\frac{1}{2k}}2\lambda t$. To derive the bound for $s_{c}$
in the second line, we use Lemma 8. Based on Eq. 41, the gate complexity of
the $2k$th-order Trotter-LCU algorithm is then

$\mathcal{O}(\nu(L+s_{c}))=\mathcal{O}\left((\lambda
t)^{1+\frac{1}{2k}}(L+\frac{\log(1/\varepsilon)}{\log\log(1/\varepsilon)})\right).$
(134)

We illustrate the random-sampling implementation of the DTSc formula in Eq.
103 and Eq. 128 in Fig. 8. In Fig. 8(b), we mainly set $K=1$ as an example due
to its simplicity. For the higher-order case, the only difference is that the
sampling of $\vec{r}$ and $\vec{r}^{\prime}$ need to be the variables
$\vec{r}_{1:\kappa}$ and $\vec{r}_{1:\kappa}^{\prime}$ on multiple stages. The
sampling probability
$\Pr(r,\vec{r}_{1:\kappa},\vec{r}^{\prime}_{1:\kappa}|s):=\mathrm{Mul}(\\{r,\vec{r}_{1:\kappa},\vec{r}^{\prime}_{1:\kappa}\\};\\{\frac{1}{2},\nu_{1:\kappa}\frac{\vec{p}}{2},\nu_{1:\kappa}\frac{\vec{p}}{2}\\};s)$
is given in Eq. 123.

Figure 8: (a) Circuit for the random-sampling implementation of the Trotter-
LCU formula. $V_{i}$ and $V_{j}$ are two random unitaries sampled
independently based on different LCU formula for $V_{K}(x)$. (b) The sampling
procedure of $V_{i}$ or $V_{j}$ in the DTSc algorithm. We set $K=1$ as
example.

For the coherent implementation of DTSc formula in Eq. 128, to apply Corollary
1, we need to set

$\mu_{2k}(x)\leq 2\Leftarrow x=\frac{1}{2\lambda}\left(\frac{\ln
2}{e}\right)^{\frac{1}{2k+1}}.$ (135)

Therefore, the number of segments $\nu$ is then

$\nu=\frac{t}{x}=2\lambda t\left(\frac{e}{\ln
2}\right)^{\frac{1}{2k+1}}=\mathcal{O}(\lambda t).$ (136)

Similar to the first-order case, the smaller value of $\mu(x)$ for each
segment is not helpful to reduce the number of segments $\nu$ for the coherent
implementation of LCU. The truncation order analysis for $s_{c}$ is similar to
the one in Sec. E.1. We also have
$s_{c}=\mathcal{O}(\log(1/\varepsilon)/\log\log(1/\varepsilon))$.

##  Appendix F Trotter-LCU formula with paired Taylor-series compensation

In this section, we discuss how to improve the Taylor-series compensation
methods introduced in Sec. E to further suppress the $1$-norm of the LCU
formulas. This will remarkably improve the $t$-dependence of the algorithms
based on random-sampling implementation. The main idea behind this improvement
is to apply Euler’s formula for Pauli operators to pair the Pauli operator
with $I$ to suppress the $1$-norm.

We remark that, similar “paired” Taylor-series compensation idea is considered
in Ref. Wan et al. (2022) for the zeroth-order case. It is also proposed in
Ref. Yang et al. (2021) for the case with Trotter formulas. Our contributions
here are: 1) we propose an explicit LCU formula with easy experimental
implementation and simple sampling methods; 2) based on the explicit LCU
formula, we provide rigorous performance guarantee and gate complexity
analysis.

###  F.1 Zeroth-order case

We start from the zeroth-order case, where the Trotter remainder
$V_{0}(x)=U(x)=e^{-ixH}$.

In Eq. 77, we have expanded the $V_{0}(x)$ based on Taylor series,

$\displaystyle V_{0}(x)$ $\displaystyle=\sum_{s=0}^{\infty}F_{0,s}(x),$ (137)
$\displaystyle F_{0,s}(x)$ $\displaystyle=\frac{(-i\lambda
x)^{s}}{s!}\sum_{l_{1:s}}p_{l_{1:s}}^{(s)}P_{l_{1:s}}^{(s)}.$

The $1$-norm of $V_{0}(x)$ is

$\|V_{0}(x)\|_{1}=\sum_{s=0}^{\infty}\|F_{0,s}(x)\|_{1}=1+(\lambda
x)+\frac{1}{2!}(\lambda x)^{2}+\frac{1}{3!}(\lambda x)^{3}+...=e^{\lambda x}.$
(138)

In the discussion in Sec. E.1, we have seen that the $1$-norm $e^{\lambda x}$
is too large for an efficient random-sampling implementation of the zeroth-
order algorithm. When $\lambda x$ is a small value, the major contribution to
$\|V_{0}(x)\|_{1}$ comes from the low-order terms of $F_{0,s}(x)$.

Note that the first-order term $F_{0,1}(x)$ is anti-Hermitian. To suppress
$\|V_{0}(x)\|_{1}$, we rewrite $V_{0}(x)$ as follows,

$\displaystyle V_{0}(x)$ $\displaystyle=\sum_{s=0}^{\infty}F_{0,s}(x)$ (139)
$\displaystyle=I-i\lambda
x\sum_{l=1}^{L}p_{l}P_{l}+\sum_{s=2}^{\infty}F_{0,s}(x)$
$\displaystyle=\sum_{l=1}^{L}p_{l}(I-i\lambda
xP_{l})+\sum_{s=2}^{\infty}F_{0,s}(x).$

Then, we can apply Eq. 69 on Eq. 139 to convert the first-order term to Pauli
rotation unitaries

$\displaystyle V_{0}(x)=\sqrt{1+(\lambda
x)^{2}}\sum_{l=1}^{L}p_{l}e^{i\theta_{0}P_{l}}+\sum_{s=2}^{\infty}F_{0,s}(x),$
(140)

where $\theta_{0}:=\tan^{-1}(-\lambda x)$.

In practice, we introduce a truncation $s_{c}\geq 2$ to the expansion order of
$x$. After this truncation, the approximated LCU formula of $V_{0}^{(p)}(x)$
is

$\displaystyle\tilde{V}_{0}^{(p)}(x)$
$\displaystyle:=\tilde{U}_{0}^{(p)}(x)=\sqrt{1+(\lambda
x)^{2}}\sum_{l=1}^{L}p_{l}e^{i\theta_{0}P_{l}}+\sum_{s=2}^{s_{c}}F_{0,s}(x)$
(141)
$\displaystyle=\mu_{0}^{(p)}(x)\left(\mathrm{Pr}_{0}^{(p)}(1)V_{0,1}^{(p)}+\sum_{s=2}^{s_{c}}\mathrm{Pr}_{0}^{(p)}(s)V_{0,s}^{(p)}\right),$

where

$\displaystyle\mu_{0}^{(p)}(x)$ $\displaystyle:=\sqrt{1+(\lambda
x)^{2}}+\sum_{s=2}^{s_{c}}\|F_{0,s}(x)\|_{1}$
$\displaystyle\leq\sqrt{1+(\lambda x)^{2}}+\left(e^{\lambda x}-1-\lambda
x\right),$ $\displaystyle\mathrm{Pr}_{0}^{(p)}(s)$
$\displaystyle:=\frac{1}{\mu_{0}^{(p)}(x)}\begin{cases}\sqrt{1+(\lambda
x)^{2}},&s=1,\\\ \|F_{0,s}(x)\|_{1}=\frac{(\lambda
x)^{s}}{s!},&s=2,3,...,s_{c},\end{cases}$ (142) $\displaystyle V_{0,s}^{(p)}$
$\displaystyle:=\begin{cases}\sum_{l}p_{l}e^{i\theta_{0}P_{l}},&s=1,\\\
\sum_{l_{1:s}}p_{l_{1:s}}^{(s)}(-i)^{s}P_{l_{1:s}}^{(s)},&s=2,3,...,s_{c}.\end{cases}$

After “pairing” the terms with $s=0$ and $s=1$, we obtain Eq. 141 which is a
new LCU formula with the $1$-norm $\mu_{0}^{(p)}(x)$.

We have the following proposition to characterize the LCU formula in Eq. 141.

######  Proposition 6 (0th-order Trotter-LCU formula by paired Taylor-series
compensation).

For $x\geq 0$ and $s_{c}\geq 2$, $\tilde{V}_{0}^{(p)}(x)$ in Eq. 141 is a
$(\mu_{0}^{(p)}(x),\varepsilon_{0}^{(p)}(x))$-LCU formula of $V_{0}(x)=U(x)$
with

$\displaystyle\mu_{0}^{(p)}(x)$ $\displaystyle\leq e^{\frac{3}{2}(\lambda
x)^{2}},$ (143) $\displaystyle\varepsilon_{0}^{(p)}(x)$
$\displaystyle\leq\left(\frac{e\lambda x}{s_{c}+1}\right)^{s_{c}+1}.$

###### Proof.

For the normalization factor, we have

$\displaystyle\mu_{0}^{(p)}(x)$ $\displaystyle\leq\sqrt{1+(\lambda
x)^{2}}+\left(e^{\lambda x}-1-\lambda x\right)$ $\displaystyle\leq
1+\frac{1}{2}(\lambda x)^{2}+\left(e^{\lambda x}-1-\lambda x\right)$
$\displaystyle=e^{\lambda x}-\lambda x+\frac{1}{2}(\lambda x)^{2}$ (144)
$\displaystyle\leq e^{(\lambda x)^{2}}+\frac{1}{2}(\lambda x)^{2}$
$\displaystyle\leq e^{\frac{3}{2}(\lambda x)^{2}}.$

The forth inequality is due to $e^{x}-x\leq e^{x^{2}}$ for $x\in\mathbb{R}$.

The proof of the distance bound $\varepsilon_{0}^{(p)}(x)$ is the same as the
one of $\varepsilon_{0}(x)$ in Proposition 3. ∎

Now, we consider the Trotter-LCU formula $\tilde{U}_{0}^{(p)}(x)$ in Eq. 141
for both random-sampling and coherent usages.

For the random-sampling implementation, the overall LCU formula for $U(t)$ is
to repeat $\tilde{U}_{0}^{(p)}(x)$ for $\nu$ times,
$\tilde{U}_{0,tot}^{(p)}(t)=\tilde{U}_{0}^{(p)}(x)^{\nu}$. Using Proposition 2
and 6, we conclude that $\tilde{U}_{0,tot}^{(p)}(t)$ is a
$(\mu_{0,tot}^{(p)}(t),\varepsilon_{0,tot}^{(p)}(t))$-LCU formula of $U(t)$
with

$\displaystyle\mu_{0,tot}^{(p)}(t)$ $\displaystyle=\mu_{0}^{(p)}(x)^{\nu}\leq
e^{\frac{3}{2}\frac{(\lambda t)^{2}}{\nu}},$ (145)
$\displaystyle\varepsilon_{0,tot}^{(p)}(t)$
$\displaystyle\leq\nu\mu_{0,tot}^{(p)}(t)\varepsilon_{0}^{(p)}(x)\leq\nu
e^{\frac{3}{2}\frac{(\lambda t)^{2}}{\nu}}\left(\frac{e\lambda
x}{s_{c}+1}\right)^{s_{c}+1}.$

To realize a $(\mu,\varepsilon)$-LCU formula for $U(t)$, we only need to set
the segment number $\nu$ and the truncation order $s_{c}$ to be

$\displaystyle\nu$ $\displaystyle\geq\frac{3}{2}\frac{(\lambda
t)^{2}}{\ln\mu},$ (146) $\displaystyle s_{c}$
$\displaystyle\geq\left\lceil\frac{\ln\left(\frac{3\mu(\lambda
t)^{2}}{2\ln(\mu)\varepsilon}\right)}{W_{0}\left(\frac{3\lambda
t}{2e\ln(\mu)}\ln\left(\frac{3\mu(\lambda
t)^{2}}{2\ln(\mu)\varepsilon}\right)\right)}-1\right\rceil,$

To derive the bound for $s_{c}$ in the second line, we use Lemma 8. Based on
Eq. 41, the gate complexity of the $0$th-order Trotter-LCU algorithm is then

$\mathcal{O}(\nu s_{c})=\mathcal{O}\left((\lambda
t)^{2}\frac{\log(1/\varepsilon)}{\log\log(1/\varepsilon)}\right).$ (147)

We remark that, comparing to other Trotter-LCU algorithms, this $0$th-order
random-sampling algorithm has a unique property that it is independent of the
number of terms $L$ in the Hamiltonian decomposition. This will make it
especially useful in the quantum chemistry tasks Campbell (2019); Wan et al.
(2022).

From Eq. 146 we can see that, different from introducing Trotter formulas in
Sec. E.2, we can also eliminate the first-order term of $x$ in the LCU formula
by constructing Pauli rotation unitaries. This enable us to consider an
efficient random-sampling implementation of Trotter-LCU formula.

For the coherent implementation, to apply Corollary 1, we need to set

$\mu_{0}^{(p)}(x)\leq 2\Leftarrow x=\frac{1}{\lambda}\sqrt{\frac{2}{3}\ln 2}.$
(148)

Therefore, the number of segments $\nu$ is then

$\nu=\frac{t}{x}=\frac{\sqrt{3}\lambda t}{\sqrt{2\ln 2}}=\mathcal{O}(\lambda
t).$ (149)

As a result, the smaller value of $\mu(x)$ for each segment is not helpful to
reduce the number of segments $\nu$ for the coherent implementation of LCU.
The truncation order analysis for $s_{c}$ is similar to the one in Sec. E.1.
We also have $s_{c}=\mathcal{O}(\log(1/\varepsilon)/\log\log(1/\varepsilon))$.

We remark that, in Ref. Wan et al. (2022), Wan et al. propose an alternative
way to construct LCU based on the paired Taylor-series compensation. Instead
of pairing the first two terms with $s=0$ and $1$, they consider to pair all
the expansion orders by grouping the $s$-th and $(s+1)$-th expansion terms
with even $s$,

$\displaystyle V_{0}(x)$
$\displaystyle=\sum_{s=0}^{\infty}F_{0,s}(x)=\sum_{s:\mathrm{even}}\frac{(-ix)^{s}}{s!}H^{s}\left(I-i\frac{x}{s+1}H\right)$
(150) $\displaystyle=\sum_{s:\mathrm{even}}\frac{(\lambda
x)^{s}}{s!}(-iH)^{s}\sum_{l^{\prime}}p_{l^{\prime}}\left(I-i\frac{x}{s+1}P_{l^{\prime}}\right).$

Applying Eq. 69 to Eq. 150, we have,

$\displaystyle V_{0}(x)$ $\displaystyle=\sum_{s:\mathrm{even}}\frac{(\lambda
x)^{s}}{s!}(-iH)^{s}\sum_{l^{\prime}}p_{l^{\prime}}\left(I-i\frac{\lambda
x}{s+1}P_{l^{\prime}}\right)$ (151)
$\displaystyle=\sum_{s:\mathrm{even}}\frac{(\lambda
x)^{s}}{s!}\sqrt{1+\left(\frac{\lambda
x}{s+1}\right)^{2}}\sum_{l_{1},...,l_{s};l^{\prime}}p_{l_{1}}...p_{l_{s}}p_{l^{\prime}}(-i)^{s}\left(P_{l_{1}}...P_{l_{s}}W_{s;l^{\prime}}\right)$

where $W_{s;l^{\prime}}:=e^{i\theta_{s}P_{l}}$ is a Pauli rotation gate with
the rotation angle $\theta_{s}:=\tan^{-1}(-\lambda x/(s+1))$. In Ref. Wan et
al. (2022), the authors prove that the $1$-norm of $V_{0}(x)$ in Eq. 151
satisfies

$\|V_{0}(x)\|_{1}\leq e^{(\lambda x)^{2}},$ (152)

which owns the same $x$-order as $\mu_{0}^{(p)}(x)$ in Eq. 145 but with an
improved overhead. This full-order-pairing approach takes advantage of the
simple structure of $V_{0}(x)$, since the $(s+1)$-order term $F_{0,s+1}(x)$
differs from the neighboring term $F_{0,s}(x)$ only by an operator of $H$. In
what follows, we are going to consider the Trotter remainders $V_{K}(x)$ which
cannot be written as a single exponential term. To this end, we will mainly
consider the “leading-order” pairing method Yang et al. (2021), where the
meaning of leading orders will be clarified soon.

###  F.2 First-order case

Now, we consider the $1$st-order Trotter remainder $V_{1}(x)$ in Eq. 100. From
Eq. 100 and Eq. 102 we can directly expand $V_{1}(x)$ in the following form,

$\displaystyle V_{1}(x)$ $\displaystyle=I+\sum_{s=2}^{\infty}F_{1,s}(x),$
(153) $\displaystyle F_{1,s}(x)$
$\displaystyle=\sum_{r+\sum\vec{r}=s}\frac{(\lambda
x)^{r}}{r!}\left(\prod\frac{(\vec{\alpha}x)^{\vec{r}}}{\vec{r}!}\right)\sum_{l_{1:r}}p^{(r)}_{l_{1:r}}(-i)^{2r-s}P^{(r)}_{l_{1:r}}\prod^{\rightarrow}\vec{P}^{\vec{r}}$
$\displaystyle=\eta_{s}\sum_{r,\vec{r}}\Pr(r,\vec{r}|s)\sum_{l_{1:r}}p_{l_{1:r}}^{(r)}(-i)^{2r-s}P^{(r)}_{l_{1:r}}\prod^{\rightarrow}\vec{P}^{\vec{r}}=\eta_{s}V_{1,s}(x).$

Recall that $\eta_{s}=(2\lambda x)^{s}/s!$. Here, the $s$th-order expansion
term $F_{1,s}(x)$ is written as a summation of Pauli operators
$P^{(r)}_{l_{1:r}}\prod^{\rightarrow}\vec{P}^{\vec{r}}$. The conditional
probability $\Pr(r,\vec{r}|s)$ is the multinomial distribution
$\text{Mul}(r,\vec{r};\\{\frac{1}{2},\frac{\vec{p}}{2}\\};s)$ defined in Eq.
96.

For the convenience of later discussion, we simplify the notation of
$F_{1,s}(x)$ as follows,

$\displaystyle
F_{1,s}(x)=\eta_{s}\sum_{r,\gamma}\Pr(r,\gamma|s)P_{1}(r,\gamma),$ (154)

where $\gamma$ is used to denote all the expansion variables
$\\{\vec{r},l_{1:r}\\}$ besides $r$. $\Pr(r,\gamma|s)$ and $P_{1}(r,\gamma)$
are then defined as

$\displaystyle\Pr(r,\gamma|s)$
$\displaystyle:=\Pr(r,\vec{r}|s)p^{(r)}_{l_{1:r}},$ (155) $\displaystyle
P_{1}(r,\gamma)$
$\displaystyle:=(-i)^{2r-s}P^{(r)}_{l_{1:r}}\prod^{\rightarrow}\vec{P}^{\vec{r}}.$

Note that, the multiplication of Pauli operators may generate imaginary global
phase, making the operator anti-Hermitian. If we want to pair the expansion
term $F_{1,s}(x)$ with $I$ using the Euler’s formula in Eq. 69, we have to
make sure that the Pauli operators $P_{1}(r,\gamma)$ with $r+\vec{r}=s$ are
all with imaginary coefficients $i$ or $-i$. This, however, does not holds in
general. One can classify the Pauli decomposition terms of $F_{1,s}$ in Eq.
153 into Hermitian and anti-Hermitian types,

$\displaystyle
F_{1,s}(x)=\eta_{s}\sum_{\\{r,\gamma\\}\in\text{Her}}\Pr(r,\gamma|s)P_{1,s}(r,\gamma)+i\eta_{s}\sum_{\\{r,\gamma\\}\in\text{anti-
Her}}\Pr(r,\gamma|s)(-i)P_{1,s}(r,\gamma),$ (156)

where $\\{r,\gamma\\}\in\text{Her}$ and $\\{r,\gamma\\}\in\text{anti-Her}$,
respectively, indicates the set of $\\{r,\gamma\\}$ such that
$P_{1}(r,\gamma)$ is Hermitian and anti-Hermitian. When
$\\{r,\gamma\\}\in\text{Her}$, the corresponding Pauli operator owns a real
coefficient, which cannot be paired with $I$ based on the Euler’s formula.

On the other hand, we can write the exponential form of $V_{1}(x)$ by applying
the BCH formula on the definition of $V_{1}(x)$ in Eq. 22,

$V_{1}(x)=\exp\left(i(E_{1,2}x^{2}+E_{1,3}x^{3}+E_{1,4}x^{4}+...)\right),$
(157)

where $\\{E_{1,s}\\}$ are some Hermitian operators. The first order term
$E_{1,1}$ vanishes due to the order condition in Lemma 4. From Eq. 157 we have

$\displaystyle F_{1,2}(x)$ $\displaystyle=ix^{2}E_{1,2},$ (158) $\displaystyle
F_{1,3}(x)$ $\displaystyle=ix^{3}E_{1,3}$

by the Taylor-series expansion on Eq. 157.

Comparing Eq. 153 and Eq. 158, we can see that

$\displaystyle\sum_{\\{r,\gamma\\}\in\text{Her}}\Pr(r,\gamma|s)P_{1,s}(r,\gamma),\quad
s=2,3.$ (159)

This is because $F_{1,2}(x)$ and $F_{1,3}(x)$ are anti-Hermitian from Eq. 158.

From Eq. 156 and Eq. 159 we can modify the form of $F_{1,s}(x)$ ($s=2$ or $3$)
as follows,

$\displaystyle F_{1,s}(x)$
$\displaystyle=\eta_{s}\sum_{\\{r,\gamma\\}\in\text{Her}}\Pr(r,\gamma|s)P_{1,s}(r,\gamma)+i\sum_{\\{r,\gamma\\}\in\text{anti-
Her}}\Pr(r,\gamma|s)(-i)P_{1,s}(r,\gamma),$ (160)
$\displaystyle=i\eta_{s}\sum_{\\{r,\gamma\\}\in\text{Her}}\Pr(r,\gamma|s)P_{1,s}(r,\gamma)+i\sum_{\\{r,\gamma\\}\in\text{anti-
Her}}\Pr(r,\gamma|s)(-i)P_{1,s}(r,\gamma),$
$\displaystyle=i\eta_{s}\sum_{r,\gamma}\Pr(r,\gamma|s)(-i)^{\mathbbm{1}[P_{1,s}(r,\gamma):\text{anti-
Her}]}P_{1,s}(r,\gamma).$

In Eq. 160, we intentionally add an extra phase $i$ on the terms with
Hermitian Pauli operators, which has no effect on $F_{1,s}(x)$ as they own
zero summation value. In this way, all the Pauli expansion terms in
$F_{1,2}(x)$ and $F_{1,3}(x)$ are with imaginary coefficient, which can be
paired with $I$ using the Euler’s formula in Eq. 69. We call the $2$nd- and
$3$rd-order terms the leading orders. Why we still keep the Hermitian terms
with zero summation value? The major reason is that we want to simplify the
sampling procedure for the variables $r,\vec{r}$ and $l_{1},...,l_{r}$ in the
random-sampling implementation: when all the Hermitian terms are kept, the
expansion variables $r,\vec{r}$ and $l_{1},...,l_{r}$ can be sampled from
$\Pr(r,\gamma)$ defined in Eq. 155. In this way, the variables
$\\{r,\vec{r}\\}$ follows a simple multinomial distribution
$\text{Mul}(r,\vec{r};\\{\frac{1}{2},\frac{\vec{p}}{2}\\};s)$ defined in Eq.
96 while $\\{l_{1},l_{2},...,l_{r}\\}$ are sampled identically and
independently from $\vec{p}$.

Now we are going to construct the LCU formula with the expansion-order
pairing. Similar to Eq. 94, we now expand $V_{1}(x)$ as

$\displaystyle V_{1}(x)$
$\displaystyle=I+\sum_{s=2}^{3}F_{1,s}(x)+\sum_{s>3}F_{1,s}(x),$ (161)
$\displaystyle F_{1,s}(x)$
$\displaystyle=i\eta_{s}V_{1,s}^{\prime}=i\eta_{s}\sum_{r,\gamma}\Pr(r,\gamma|s)(-i)^{\mathbbm{1}[P_{1,s}(r,\gamma):\text{anti-
Her}]}P_{1,s}(r,\gamma),\quad s=2,3$ $\displaystyle F_{1,s}(x)$
$\displaystyle=\eta_{s}V_{1,s}=\eta_{s}\sum_{r,\gamma}\Pr(r,\gamma|s)(-i)^{2r-s}P_{1,s}(r,\gamma),\quad
s\geq 4.$

As is introduced in Eq. 160, we modify the phase of the Hermitian terms in
$F_{1,2}(x)$ and $F_{1,3}(x)$ for the later usage of pairing. In this way, the
operator

$P_{1,s}^{\prime}(r,\gamma):=(-i)^{\mathbbm{1}[P_{1,s}(r,\gamma):\text{anti-
Her}]}P_{1,s}(r,\gamma),\quad s=2,3$ (162)