# High-precision and low-depth eigenstate property estimation:
theory and resource estimation

Jinzhao Sun<EMAIL_ADDRESS>Blackett Laboratory, Imperial College
London, London SW7 2AZ, United Kingdom  Computer Laboratory, University of
Cambridge, Cambridge CB3 0FD, United Kingdom  Pei Zeng<EMAIL_ADDRESS>Pritzker School of Molecular Engineering, The University of Chicago, Illinois
60637, USA  Tom Gur<EMAIL_ADDRESS>Computer Laboratory, University of
Cambridge, Cambridge CB3 0FD, United Kingdom  M. S. Kim<EMAIL_ADDRESS>Blackett Laboratory, Imperial College London, London SW7 2AZ, United Kingdom

###### Abstract

Estimating the eigenstate properties of quantum many-body systems is a long-
standing, challenging problem for both classical and quantum computing. For
the task of eigenstate preparation, quantum signal processing (QSP) has
established near-optimal query complexity
$\mathcal{O}(\Delta^{-1}\log(\varepsilon^{-1}))$ by querying the block
encoding of the Hamiltonian $H$ where $\Delta$ is the energy gap and
$\varepsilon$ is the target precision. However, QSP is challenging for both
near-term noisy quantum computers and early fault-tolerant quantum computers
(FTQC), which are limited by the number of logical qubits and circuit depth.
To date, early FTQC algorithms have focused on querying the perfect time
evolution $e^{-iHt}$. It remains uncertain whether early FTQC algorithms can
maintain good asymptotic scaling at the gate level. Moreover, when considering
qubit connectivity, the circuit depth of existing FTQC algorithms may scale
suboptimally with system size. Here, we present a full-stack design of a
random sampling algorithm for estimating the eigenenergy and the observable
expectations on the eigenstates, which can achieve high precision and good
system size scaling. The gate complexity has a logarithmic dependence on the
inverse precision ${\mathcal{O}}(\log^{1+o(1)}(1/\varepsilon))$ for generic
Hamiltonians, which cannot be achieved by methods using Trottersiation to
realise $e^{-iHt}$ like in QETU. For $n$-qubit lattice Hamiltonians, our
method achieves near-optimal system size dependence with the gate complexity
$\mathcal{O}(n^{1+o(1)})$. When restricting the qubit connectivity to a linear
nearest-neighbour architecture, The method shows advantages in circuit depth,
with $\mathcal{O}(n^{o(1)})$ for lattice models and $\mathcal{O}(n^{2+o(1)})$
for electronic structure problems. We compare the resource requirements (CNOT
gates, T gates and qubit numbers) by state-of-the-art methods, including phase
estimation, QSP, and QETU, in tasks of lattice models and molecular problems.
We provide a toolbox for analysing resource requirements for these advanced
methods when compiled at the gate level.

##  I Introduction

Estimating the properties of the ground state and excited states of quantum
many-body systems is a long-standing problem of fundamental interest, which
has applications in condensed matter physics, quantum chemistry and material
science bharti2021noisy ; mcardle2020quantum ; bauer2020quantum ;
dalzell2023quantum . However, this is a challenging task for both classical
computing and quantum computing. Indeed, estimating the ground state energy of
a two-local Hamiltonian is a QMA-complete problem kempe2006complexity  (see
numerical evidence in lee2023evaluating ). Nevertheless, quantum computing is
a promising paradigm for addressing this problem under certain assumptions,
which include a non-vanishing energy gap $\Delta$ and a non-vanishing overlap
of the initial state and the target state. There has been considerable
progress in the development of quantum computing algorithms (see
dalzell2023quantum  for a review), which include quantum phase estimation
(QPE) kitaev1995quantum ; higgins2007entanglement ; knill2007optimal ;
rall2021faster ; meister2022resource ; martyn2021grand ; ding2023even ,
spectral filter algorithms (e.g., spectral filter by a linear combination of
unitaries (LCU) Ge19 ; keen2021quantum ; an2023linear , random sampling
yang2021accelerated ; wan2021randomized ; lin2021heisenberglimited ;
zeng2021universal ; lu2021algorithms ; zhang2022computing ; huo2021shallow ;
wang2023quantum ; ding2024quantum ; he2022quantum ; wang2023faster ;
wang2024qubit , quantum signal processing (QSP) lin2020near ; lin2020optimal ,
quantum singular value transformation (QSVT) gilyen2019quantum ;
martyn2021grand ), dissipation-based algorithms cubitt2023dissipative ;
ding2023single ; chen2023local , variational algorithms McClean2016theory ;
cerezo2021variational ; liu2024variational , and others huggins2022unbiasing ;
feng2024escaping ; kovalsky2023self ; motta2020determining ; xu2023quantum ;
zhang2022quantumQMC ; hejazi2023adiabatic . For the task of eigenstate
preparation, spectral filter-based methods provide a rigorous and
deterministic solution with clear assumptions of the initial overlap and
energy gap, which will be the major focus of this work.

In spectral filter methods, QSP has achieved near-optimal query complexity
$\mathcal{O}(\Delta^{-1}\log(\varepsilon^{-1}))$ for target precision
$\varepsilon$, which is a widely accepted as one of the state-of-the-art
algorithms lin2020near . While it is favourable in the long term, QSP hinges
on querying the block encoding of the Hamiltonian $H$, which is challenging
for near-term noisy intermediate-scale quantum (NISQ) or early fault-tolerant
quantum computing (FTQC) applications. In the early FTQC regime
lin2021heisenberglimited ; katabarwa2023early , where the number of logical
qubits is limited, minimising controlled operations and achieving low circuit
depth are essential goals. Considering physical devices, qubit connectivity
must be taken into account, as depth and controlled operations are often
interrelated. Compiling a nonlocal controlled gate into nearest-neighbour
gates devulapalli2022quantum  can significantly increase both the required
qubit count and the overall circuit depth. Therefore, it is desirable to
design quantum algorithms tailored for NISQ and early FTQC, which satisfy the
constraints of fewer qubit numbers and low circuit depth.

Considering the feature of early FTQC lin2021heisenberglimited , there has
been considerable progress on quantum algorithms for ground-state property
estimation wang2023faster ; lin2021heisenberglimited ; zeng2021universal ;
lu2021algorithms ; zhang2022computing ; huo2021shallow ; he2022quantum ;
wang2023quantum ; ding2024quantum ; ding2023even ; wang2023faster . In 2021,
Lin and Tong proposed a spectral filter algorithm based on random sampling
lin2021heisenberglimited , which achieved the Heisenberg limit for ground-
state energy estimation. Using similar techniques, Zeng et al.
zeng2021universal  and Zhang et al. zhang2022computing  showed that the time
complexity can achieve logarithmic in precision $\log(1/\varepsilon)$ for the
ground state property estimation problem 111Note that Ref. zeng2021universal
extends to eigenstate property estimation. . On the other hand, following the
spirit of QSP, Dong et al. proposed the quantum eigenvalue transformation of
unitaries (QETU) algorithm dong2022ground  to prepare the ground state. QETU
achieves near-optimal asymptotic scaling but avoids querying the block-
encoding of the Hamiltonian by querying time evolution, making it appealing to
early fault-tolerant quantum computers. The above algorithms assume the usage
of only one ancillary qubit and achieve good asymptotic sample and time
complexity, which are competitive for the application of noisy intermediate-
scale quantum and early FTQC. However, these algorithms assume perfect and
efficient queries to the time evolution operator $e^{-iHt}$. It remains
unclear whether the good properties of the above early FTQC algorithms can be
preserved if we further expand $e^{-iHt}$ into elementary gates. For example,
if $e^{-iHt}$ is implemented using Trotterisation methods, this will eliminate
the advantage of logarithmic scaling in precision. On the other hand, if we
introduce advanced block-encoding-based methods to realise $e^{-iHt}$, it
needs many ancillary qubits and nonlocal controlled gates, which violates the
spirit of early FTQC. Moreover, the system-size dependence of the algorithms
is rarely discussed in existing works, as it highly depends on the detailed
circuit-level implementation of $e^{-iHt}$ as well as the qubit connectivity
of the device. This raises the question when considering the gate complexity
and qubit connectivity, whether we can design high-precision and low-depth
algorithms for NISQ and early FTQC applications.

In this work, we provide an efficient quantum algorithm from end to end for
the task of eigenstate property estimation by taking advantage of the recent
advances in spectral filters and Hamiltonian simulations. We consider
realising a non-unitary spectral filter through real-time evolution,
accompanied by a high-precision Hamiltonian simulation scheme. To do so, we
build up a framework based on random composite LCU formulae for analysing the
gate complexity for the task of eigenstate property estimation. It contains
several hierarchies. We consider realising the non-unitary operator by random
sampling real-time evolution $e^{-iHt}$ with different time lengths $t$
lu2021algorithms ; zeng2021universal ; ding2024quantum ; wang2023quantum ;
zhang2022computing ; huo2021shallow ; he2022quantum ; lu2021algorithms ;
wang2023faster . To preserve high precision when implementing time evolution,
we choose to compensate the Trotter error zeng2022simple  by using a composite
LCU formula, as opposed to only implementing the Trotter formulae in
conventional Trotter methods. The first contribution of this paper is that the
gate complexity for each circuit (at a single run) is shown to be
$\mathcal{O}(\operatorname{poly}\log(\varepsilon^{-1}))$, outperforming the
QPE-based method and matching the result by QSP. We note that QETU can achieve
high-precision ground state preparation with one ancillary qubit. However,
QETU hinges on a coherent implementation of $e^{-iHt}$, which rules out any
random sampling method, and thus, it is not straightforward to achieve
logarithmic dependence on inverse precision. Furthermore, we can achieve near-
optimal system size scaling for lattice Hamiltonians by exploiting the
commutation relations. In this work, we mainly discuss the use of $2k$th-order
Trotter formula and then compensate the Trotter error. For $n$-qubit lattice
Hamiltonians, the total gate count is shown to be
$\mathcal{O}(n^{1+\frac{2}{4k+1}}\varepsilon^{-\frac{1}{4k+1}})$ where $2k$ is
the order of the Trotter formula used in this work. Our result shows better
system size dependence compared with QSP, and a better precision scaling
compared with QETU, as can be found in Table 1.

In addition to addressing gate complexity for generic Hamiltonians, the second
contribution of our work is the advantages in circuit depth for various
physical problems. These include lattice models and second-quantised plane-
wave electronic structures with $n$ spin orbitals and $L=\mathcal{O}(n^{2})$
terms in the Hamiltonian. Notably, the 2D Fermi-Hubbard model falls within the
problem class. In theory, for problems that conserve certain symmetries, we
demonstrate the ability to achieve high-precision and low-depth eigenstate
property estimation without the need for ancillary qubits. To accomplish this,
we design new Trotter and Trotter-error compensation circuits that maintain
the system’s symmetries. Thanks to this ancilla-free feature, when restricting
the qubit connectivity, the random sampling method shows advantages in circuit
depth over existing best practices. When restricted to nearest-neighbour
architecture, for lattice models, the circuit depth scales as
$d=\mathcal{O}(n^{\frac{2}{4k+1}})$, while the circuit depth for QETU combined
with Trotter scales as $d_{\mathrm{QETU}}=\mathcal{O}(n^{1+\frac{1}{2k}})$.
For second-quantised plane-wave electronic structures, the circuit depth is
$d=\mathcal{O}(n^{2+\frac{2}{4k+1}})$, while naive implementation of the QETU-
Trotter method requires $d_{\mathrm{QETU}}=\Omega(n^{3+\frac{1}{2k}})$.

To address the application in the NISQ and FTQC eras, we present the actual
gate count and number of qubits required for solving practical problems
accordingly. In the existing resource estimation works kim2022fault ;
reiher2017elucidating ; goings2022reliably ; wecker2014gate ;
campbell2019random , the energy estimation or eigenstate preparation for
practical problems, such as second-quantised quantum chemistry problems with
$L=\mathcal{O}(n^{4})$ terms lee2021even ; babbush2018encoding ;
berry2019qubitization  and condensed-phase electrons kivlichan2020improved ;
campbell2021early , is predominantly based on phase estimation. A typical
strategy is to encode the eigenspectra of the Hamiltonian in a unitary for
phase estimation by the evolution $e^{-iHt}$, which is synthesised by
Trotterisation kivlichan2020improved ; campbell2021early , or a qubitised
quantum walk babbush2018encoding ; lee2021even  with eigenspectrum
proportionally to $e^{\pm i\arccos(H/\lambda)}$, where $\lambda$ is a
parameter related to the norm of the Hamiltonian. However, due to the cost of
phase estimation, the circuit depth will inevitably be polynomial in the
precision, which is not optimal for eigenstate energy and property estimation.

The third contribution is to present resource estimations with the method
developed in this work, which has better asymptotic scaling in precision and
system size in theory, and also has a smaller circuit compilation overhead in
practice. In the NISQ applications, the major bottleneck is the number of two-
qubit gates for the noisy quantum computer. The major overhead for the error-
corrected quantum computer is the T gate count babbush2018encoding , which
requires more gates to perform error correction than Clifford gates.
Therefore, we study the resource requirement in terms of CNOT and T gates for
ground-state property estimation for representative physical models in
condensed matter and chemistry. We provide a systematic comparison with the
existing advanced methods, which include QPE combined with Trotter formulae or
qubitised quantum walk (as have been used in lee2021even ; babbush2018encoding
; goings2022reliably ), QSP lin2020near  and QETU. The CNOT gate cost for a
20-site Heisenberg model is about $3\times 10^{5}$ while the T gate cost is
about $6\times 10^{6}$. A notable contribution of this work is that for
different query models, we provide a useful toolbox for researchers to analyse
the individual costs for elementary units (like block encoding and controlled
real-time evolution), and thus enable comparison across different eigenstate
preparation methods with various initial conditions. We hope that our
framework (decomposing the task into the elementary operations) and the
toolbox for analysing the cost for each elementary operation can be useful as
a building block for resource estimations for other quantum algorithms based
on querying the block encoding of $H$ or the time evolution.

Figure 1:  (a) The random composite LCU formulae components of $g_{\tau}$.
(a1) The spectral filter is decomposed into a combination of real-time
evolution $U(t_{i})=e^{-iHt_{i}}$, as described by Eq. 1. (a2) The real-time
evolution is divided into $\nu$ segments. Each segment with time duration
$t_{i}/\nu$ is realised by a Trotter formula $S$ and the Trotter remainder
$V$. The Trotter remainder $V$ can be written into an LCU form in Eq. 3. (a3)
A random composite LCU formula of $g_{\tau}$. The summands of the circuit are
shown here (each circuit instance is described by Eq. 35). In the one-ancilla
scheme $W_{i_{q}}(q\leq\nu)$ is composed of a Pauli rotation operator
$R=\exp(-i\theta_{i_{q}}P_{i_{q}})$ and a tensor product of single-qubit Pauli
operators, where $P_{i_{q}}$ is drawn from a certain probability distribution.
The maximum gate count in $W_{i_{q}}$ is $\operatorname{wt_{m}}{H}+n$, where
$\operatorname{wt_{m}}{H}$ is the largest weight of the Hamiltonian terms. The
controlled $W$ operation is shown in (c3). In the ancilla-free measurement
scheme illustrated in (c2), $W_{i_{q}}$ is chosen to be symmetry-conserved
operators to conserve the symmetry of the real-time evolution. (b1) The
spectral weight of the effectively imaginary-time-evolved state after applying
the spectral filter $g_{\tau}(H-E_{0})$ where $E_{0}$ is assumed to be known.
The $y$ axis represents the spectral weight
$\braket{u_{i}}{g_{\tau}(H-E_{0})}{\psi_{0}}$, where $|u_{i}\rangle$ is the
eigenstate with energy $E_{i}$. By increasing the imaginary time $\tau$, only
the spectral weight on the ground state will be preserved. The contributions
from the other eigenstates are suppressed exponentially with $\tau$. (b2) The
eigenenergy search using
$D_{\tau}(\omega)=\braket{\psi_{0}}{g_{\tau}^{2}(H-\omega)}{\psi_{0}}$ where
$D_{\tau}(\omega)$ is defined in Eq. 11. By increasing the imaginary time
$\tau$, the peaks become sharper (from the blue line to the red line) and the
true eigenenergies can be resolved by finding the peaks. The key quantity
encountered in (b1) and (b2) is $\braket{U^{\dagger}(t_{j})OU(t_{i})}$ in Eq.
36. The component in $\braket{U^{\dagger}(t_{j})OU(t_{i})}$ can be measured
using the circuit in (c). There are two extreme cases: Case I: $t_{j}=0$, and
Case II $t_{i}=t_{j}$. When $t_{j}=0$, the circuit reduces to the one in the
green box. When $t_{j}=t_{i}$, the circuit reduces to the one in the orange
box. Case I requires more resources than other cases. The resources in other
cases with $t_{j}<t_{i}$ will be between Case I and Case II. (c2) Ancilla-free
measurement scheme for symmetry-conserved cases. $U_{p}$ is the unitary
operator for state preparation $\ket{\psi_{0}}=U_{p}\ket{0^{\otimes n}}$ or
$\frac{1}{\sqrt{2}}(\ket{\psi_{\textrm{Ref}}}+\ket{\psi_{0}})=U_{p}\ket{0^{\otimes
n}}$ which is discussed in subsubsection III.4.2. (c3) For general cases in
(c1), the controlled $S$ operation takes the form of $\exp(-it_{i}P/\nu)$ with
Pauli operator $P$, which can be compiled into controlled single-qubit Pauli
rotation gates, CNOT gates and single-qubit Pauli gates.

The rest of this paper is organised as follows. In Sec. II, we outline the
algorithm and overview the main results in this work. In Sec. III, we show the
construction of a spectral filter with the random-sampling composite LCU
formula. We show the methods for eigenenergy estimation and eigenstate
property estimation, and provide the circuit depth and gate complexity for the
two tasks. In Sec. IV, we present the resource cost in terms of the ancillary
qubits, CNOT gates and non-Clifford gates (T gates and single-qubit Pauli
rotation gates) for various physical systems, including the lattice models and
quantum chemistry problems.

##  II Framework and main results

In this section, we introduce the framework of eigenenergy and eigenstate
property estimation with the random-sampling composite linear-combination-of-
unitary (LCU) formulae. Let us consider an $n$-qubit gapped quantum system
with the Hamiltonian $H$. The eigenstate $\ket{u_{j}}$ and the corresponding
eigenenergy $E_{j}$ of the Hamiltonian satisfy,
$H\ket{u_{j}}=E_{j}\ket{u_{j}},$ for $j\leq 2^{n}$. The tasks concerned in
this work are (1) to estimate the eigenenergy $E_{j}$, and (2) to estimate an
eigenstate property, characterised by an observable expectation on the target
eigenstate $\braket{u_{j}}{O}{u_{j}}$. The problems with the assumptions and
the notations used in this work are stated below.

######  Problem 1 (Eigenstate property estimation).

Consider a gapped $n$-qubit Hamiltonian with a Pauli decomposition
$H=\sum_{l=1}^{L}\alpha_{l}P_{l}:=\lambda\sum_{l=1}^{L}\tilde{\alpha}_{l}P_{l}$
where $P_{l}$ is a Pauli operator, $\lambda:=\sum_{l}|\alpha_{l}|$, and
$\tilde{\alpha}:=\alpha_{l}/\lambda$. The task is to estimate the expectation
value of an observable $O$ on the $j$th eigenstate of $H$, $\ket{u_{j}}$. We
assume that we have an estimate of $E_{j}$, $\hat{E}_{j}$ with a small
estimation error $\kappa:=|\hat{E}_{j}-E_{j}|$. We assume that we have a good
initial state $\ket{\psi_{0}}$ which has a nonvanishing overlap with the
target eigenstate,
$\eta:=|\braket{\psi_{0}}{u_{j}}|^{2}=\Omega(1/\operatorname{poly}(n))$. We
assume a nonvanishing energy gap $\Delta:=\min(E_{j+1}-E_{j},E_{j}-E_{j-1})$.

Suppose the observable has a Pauli decomposition as
$O=\sum_{l=1}^{L_{o}}o_{l}P_{l}$ with Pauli operators $P_{l}$ and positive
coefficients $o_{l}$, and $\|O\|_{1}:=\sum_{l=1}^{L_{o}}|o_{l}|$. The aim is
to find an estimator $\hat{v}$, such that it is close to
$\braket{u_{j}}{O}{u_{j}}$ with at least probability $1-\vartheta$, i.e.,
$\Pr(|\hat{v}-\braket{u_{j}}{O}{u_{j}}|\leq\varepsilon)\geq 1-\vartheta$.

Note that in Problem 1, we follow the convention in zeng2021universal , which
has a slight difference in the definition of the initial overlap $\eta$
compared to Refs. lin2020near ; dong2022ground .

######  Problem 2 (Eigenenergy estimation).

Under the same definitions and assumptions introduced in Problem 1, the aim is
to find an eigenenergy estimator $\hat{E_{j}}$, such that it is $\kappa$-close
to $E_{j}$ with at least probability $1-\vartheta$, i.e.,
$\Pr(|\hat{E_{j}}-E_{j}|\leq\kappa)>1-\vartheta$.

The design of the eigenstate property estimation algorithm is sketched in
Algorithm 1 and illustrated in Fig. 1. To access the physical properties of
eigenstates, a natural idea is to apply a spectral filter to the initial state
that projects out the contributions from the other unwanted eigenstates, as
illustrated in Fig. 1(b1). While the spectral filter $g$ is non-unitary by
construction, we can effectively realise it using LCU techniques, either by
coherent Ge19 ; keen2021quantum  or random-sampling implementation
zeng2021universal ; lu2021algorithms ; zhang2022computing ; huo2021shallow ;
he2022quantum ; wang2023quantum ; yang2021accelerated . The decomposition of
$g$ into real-time evolution $U(t_{i})=e^{-iHt_{i}}$ with time length $t_{i}$
is shown in Fig. 1(a1). We further realise each $U(t_{i})$ by another random
LCU formula. Within each time segment, the LCU formula involves both the
Trotter formula $S$ and the Trotter-error-compensation term $V$, as shown in
Fig. 1(a2). Overall, it forms a composite LCU formula illustrated in Fig.
1(a3), which involves the summation and product of individual LCU components.
The hierarchy of different LCU components is illustrated in Fig. 2. The error
propagation will be analysed in Sec. III and Appendix in order to prove the
main theorems. The quantum circuit realisation is shown in Fig. 1(c),
including both one-ancilla and ancilla-free schemes. A more detailed
Hamiltonian-dependent circuit compilation method will be discussed in
subsection III.4. Below, we elaborate on the composite LCU formula for
eigenstate property estimation.

The spectral filter $g$ is a non-unitary operator defined on the $n$-qubit
system, which is usually a function of the target Hamiltonian $H$. Choices for
spectral filters include the imaginary-time evolution operator
$g_{\tau}(H)=e^{-\tau H}$ or the Gaussian operator
$g_{\tau}(H)=e^{-\tau^{2}H^{2}}$. For example, by applying
$g_{\tau}(H)=e^{-\tau|H-E_{0}|}$ to an initial state $\ket{\psi_{0}}$, the
state becomes
$\ket{\psi(\tau)}=c(\tau)\sum_{i}c_{i}e^{-|E_{i}-E_{i}|}\ket{u_{i}}$ where
$c(\tau)$ is the normalisation factor. Provided the assumptions in Problem 1,
i.e., a nonvanishing energy gap $\Delta=E_{1}-E_{0}$ and a nonvanishing
$c_{0}=\eta^{-1/2}$, the spectral weight of the unwanted excited eigenstates
is suppressed with increasing $\tau$. In the large $\tau$ limit, the
unnormalised state becomes
$\lim_{\tau\rightarrow\infty}|\psi(\tau)\rangle=c_{0}\ket{u_{0}}$, as
illustrated in Fig. 1(b1).

We note that the spectral filter method has been well-established in the
existing literature, in particular lin2020near ; zeng2021universal ;
zhang2022computing ; wang2023quantum . In this paper, we introduce a full-
stack approach to eigenstate property and energy estimation by random-sampling
composite LCU decomposition of a Gaussian operator, and we refer to this full-
stack random-sampling LCU approach as our method for simplicity.

Now, let us discuss the construction of the composite LCU formulae. A
$(\mu,\varepsilon)$-random-sampling linear-combination-of-unitary (RLCU)
formula of a general operator $g$ is defined to be

$\tilde{g}=\mu\sum_{i}\Pr(i)U_{i},$ (1)

such that the spectral norm distance $\|g-\tilde{g}\|\leq\varepsilon$, as
illustrated in Fig. 1(a1). Here, $\mu>0$ is the normalisation factor, $\Pr(i)$
is a probability distribution associated with an instance specified by $i$,
and $\\{U_{i}\\}_{i}$ is a group of unitaries. Eq. 1 can be extended to a
continuous form

$\tilde{g}=\mu\int_{-x_{c}}^{+x_{c}}dxp(x)U(x)$ (2)

where we require that $p(x)~{}(x\in[x_{c},x_{c}])$ is a well-defined
probability distribution.

As discussed above, a natural choice for $U_{i}$ is the real-time evolution
$U_{i}=U(t_{i}):=e^{-iHt_{i}}$ with time length $t_{i}$ because quantum
systems in nature evolve governed by the Hamiltonian; for instance, analogue
quantum simulators can directly perform $U(t_{i})$. The Gaussian spectral
filter can be decomposed on the basis of real-time evolution
$g(H)=\int^{\infty}_{-\infty}dp(t)e^{-itH}$ with a Gaussian probability
distribution $p(t)=\frac{1}{\sqrt{2\pi}}e^{-t^{2}/4}$. We could truncate the
time range and formally express it as $g=\sum_{i}\Pr(i)e^{-iHt_{i}}.$ Though
the Gaussian spectral filter is represented in an integral form, it has a
well-defined probability distribution and can be well-characterised by Eq. 1.
Therefore, we will mainly discuss the property of composite LCU form in the
discretised form, while we use the continuous form when dealing with the
Gaussian spectral filter. While most early FTQC algorithms assume that real-
time evolution can be implemented perfectly, $U(t_{i})$ in general is not
directly implementable. An established way to implement $U(t_{i})$ without
involving other oracles is through Trotter formulae: when the Hamiltonian can
be decomposed into Pauli operators, each $U(t_{i})$ can be implemented using
elementary quantum gates without querying any other oracles and introducing
additional gates. However, by using Trotter formulae, the gate complexity will
be polynomial in the inverse accuracy childs2021theory . Therefore, any
ground-state property estimation protocol based on Trotterisation will lose
the advantage of high precision in the first place. The question is how to
preserve the logarithmic dependence on the precision.

Given a Trotter formula $S(t_{i})$, the Trotter error can be expressed as
$V(t_{i})=U(t_{i})S^{\dagger}(t_{i})$. The problem with the Trotter method is
that the remainder of a $k$th Trotter formula (i.e., Trotter error) is non-
negligible, which is polynomial in the order of $k$. To preserve the high
precision property, we choose to implement the Trotter remainder $V$ as well,
as opposed to implementing $S$ only in conventional Trotter methods, which is
illustrated in Fig. 1(a2,a3). By doing so, we can implement a Trotter-error-
free $U(t_{i})$. Specifically, we consider dividing the time evolution
operator $U(t_{i})$ into $\nu$ segments, which can be written as
$U(t_{i})=\left(S(t_{i}/\nu)V(t_{i}/\nu)\right)^{\nu}$ where $S(t_{i}/\nu)$ is
a Trotter formula and $V(t_{i}/\nu)$ is the corresponding Trotter error within
duration $t_{i}/\nu$. The key idea of this Trotter error compensation is that
the Trotter error can be further decomposed into easy-to-implement unitary
operators

$V=\mu\sum_{j}\Pr(j)W_{j}$ (3)

where $\Pr(j)$ is the probability associated with the unitary operator $W_{j}$
in the decomposition of $V$, and $\mu$ is the normalisation factor.
Importantly, the decomposition in Eq. 3 can be explicitly derived by
$V=US^{\dagger}$ using the Taylor expansion. By doing so, we can compensate
the Trotter error by sampling $W_{j}$ according to $\Pr(j)$ in Eq. 3. This
effectively realises a Trotter-error-compensated real-time dynamics, with the
individual terms illustrated in Fig. 1(a3). In general cases, the elementary
operator $W_{j}$ in Eq. 3 can be chosen as Pauli operators or the
exponentiation of Pauli operators, as commonly used in Trotter and LCU methods
childs2018toward ; zeng2022simple . However, choosing Pauli operators as the
bases may break the symmetry of the real-time evolution. To enable ancilla-
free measurement, we choose $W_{j}$ to be symmetry-conserved operator.
Specifically, to conserve the symmetry of either particle number or total
spins, $W_{j}$ is chosen to be SWAP, a tensor product of Pauli-$Z$ operators,
and their exponentiation, which will be discussed in subsection III.3.

The spectral filter can now be formally rewritten as

$\tilde{g}=\sum_{i}\Pr(i)\left(V(\frac{t_{i}}{\nu})S(\frac{t_{i}}{\nu})\right)^{\nu}$
(4)

with the aim of being close to the original spectral filter. The spectral
filter consists of other LCU formulae in it. We could see that Eq. 4 is a
modified version of the original LCU formula given by Eq. 1, which involves
hierarchical formulations, as illustrated in Fig. 1(a). Specifically, it
adopts a composite structure that integrates both the products of individual
LCU components and the summations of LCU components. This work will analyse
the propagation of errors and the change of normalisation factors within this
composite LCU framework, outlined in Fig. 2.

We will need to introduce a few tools to analyse the properties of the
composite form of LCU formulae.

######  Proposition 1 (Composite LCU formula).

Suppose $g_{1}$ takes a $(\mu_{1},\varepsilon_{1})$-LCU formula of $g$,
$g_{1}=\mu_{1}\sum_{i}\Pr(i)U(t_{i})$ in Eq. 1. Suppose that each of the
summands $U(t)$ is real-time evolution and is divided into $\nu$ segments,
with each segment $U(\delta t)$, $\delta t=t/\nu$. Suppose that $U(\delta t)$
has a $(\mu_{2},\varepsilon_{2})$-LCU formula, $\tilde{U}(\delta
t)=\mu_{2}\sum_{\ell}\Pr(\ell)P_{\ell}$. Then the formula

$g_{2}=\mu_{1}\sum_{i}\Pr(i)(\tilde{U}(t_{i}/\nu))^{\nu}$ (5)

is a $(\mu,\varepsilon)$-LCU formula of $g$, with $\mu:=\mu_{1}\mu_{2}^{\nu}$,
and $\varepsilon\leq\varepsilon_{1}+\mu_{1}\mu\varepsilon_{2}$.

The proof of Proposition 1 is shown in Appendix A.1. Leveraging Proposition 1,
we could calculate the error in the composite LCU form provided the error of
each individual LCU approximation.

The key steps of the eigenstate property estimation algorithm and the circuit
implementation are summarised in Algorithm 1 below.

1:Construct a composite LCU formula as described in Proposition 1. We first
expand a Gaussian spectral filter $g$ into a combination of real-time
evolution $U(t_{i})$ by Eq. 2. We decompose
$U(t_{i})=(S_{2k}(t_{i}/\nu)V_{2k}(t_{i}\nu))^{\nu}$ into Pauli operators,
which takes a composite LCU form, with the LCU form of $V$ in Eq. 3 and the
explicit form in Eq. 16 in Sec. III). Here $S_{2k}$ is the $2k$th-order
Trotter formula and $V_{2k}$ is the Trotter remainder. We consider effectively
realising both the Trotter $S_{2k}$ and the Trotter remainder $V_{2k}$ to
enable high-precision simulation (where Trotter error is compensated).

2:Realise the composite LCU formula by random sampling. We first generate
$t_{i}$ and $t_{j}$ (at the level of real-time evolution) and $\vec{i}$ and
$\vec{j}$ (at the level of gates, see Eq. 35) by sampling from the
corresponding probability distributions, which gives us an explicit form of
the circuit in each run (see Eq. 35 and Eq. 36 for example). The circuit for
evaluating $N(O)$ and $D$ (in Eq. 30) is shown in Fig. 1(c1,c2). The ancilla-
free scheme is discussed in subsection III.3.

3:Estimate the observable expectation value by
$\braket{O}=\frac{N_{\tau}(O)}{D_{\tau}}$, and eigenenergy by
$E_{j}=\operatorname*{arg\,max}_{\omega}D_{\tau}(\omega)$. The estimator of
$N_{\tau}(O)$ and $D_{\tau}$ is given by Eq. 40. The observable expectation on
the eigenstate is given by Eq. 11. As illustrated in Fig. 1(b) The eigenenergy
could be determined by $\max_{\omega}D_{\tau}(\omega)$ where the search
process is evaluated by pure classical computation. The eigenenergy can be
estimated with $\hat{D}_{\tau}(E)$ where the search process with different $E$
can be evaluated by pure classical computation.

Algorithm 1 An overview of the random-sampling algorithm for eigenstate property and eigenenergy estimation described in Problem 1 and Problem 2.  Methods | Gate complexity | Gate complexity (lattice models) | Extra qubits
---|---|---|---
QPE + Trotter ($2k$th-order) | ${\mathcal{O}}(L\tilde{\Delta}^{-(1+\frac{1}{2k})}\varepsilon^{-(1+\frac{1}{k})})$ | ${\mathcal{O}}(n^{1+\frac{1}{2k}}\Delta^{-(1+\frac{1}{2k})}\varepsilon^{-(1+\frac{1}{k})})$ | $\log(\varepsilon^{-1})+\log(\Delta^{-1})$
QPE + QW lee2021even  | $\tilde{\mathcal{O}}(L\tilde{\Delta}^{-1}\varepsilon^{-1})$ | ${\mathcal{O}}(n^{2}\Delta^{-1}\varepsilon^{-1})$ | $\log(L)+\log(\varepsilon^{-1})+\log(\Delta^{-1})$
QSP lin2020near  | ${\mathcal{O}}(L\tilde{\Delta}^{-1}\log\varepsilon^{-1})$ | ${\mathcal{O}}(n^{2}\Delta^{-1}\log\varepsilon^{-1})$ | $\log(L)+\log\log(\varepsilon^{-1})+\log(\Delta^{-1})$
QETU dong2022ground  | $\tilde{\mathcal{O}}(L\tilde{\Delta}^{-(1+\frac{1}{2k})}\varepsilon^{-\frac{1}{2k}})$ | ${\mathcal{O}}(n^{2}\Delta^{-(1+\frac{1}{2k})}\log\varepsilon^{-1})$ | 1
This work ($2k$th-order) | $\mathcal{O}(L(\tilde{\Delta}^{-1}\log\varepsilon^{-1})^{(1+\frac{1}{4k+1})})$ | $\tilde{\mathcal{O}}(n^{1+\frac{2}{4k+1}}\Delta^{-(1+\frac{1}{4k+1})}\varepsilon^{-\frac{1}{4k+1}})$ | 0 or 1
(zeroth-order) | $\mathcal{O}(\tilde{\Delta}^{-2}\log^{2}\varepsilon^{-1})$ | $\mathcal{O}(n^{2}\Delta^{-2}\log^{2}\varepsilon^{-1})$ | 0 or 1
Table 1: Comparison of observable estimation on the eigenstate of a Hamiltonian (Problem 1). Here we compare the gate complexity of the algorithm in a single coherent run. The results in the second column in the table are based on Algorithm 1 and Theorem 1 using the $2k$th-order Trotter error compensation. Here, $\tilde{\Delta}:=\lambda/\Delta$. The results in the third column in the table are based on Theorem 2 when the $2k$th-order Trotter formula is used. The eigenenergy is assumed to be known a priori and the unknown eigenenergy case is analysed in subsection B.5. The dependence on $\eta$ is not included in the table since it only appears in sample complexity for our method. Similar to other random-sampling spectral filter methods (see zhang2022computing ; lin2021heisenberglimited ), the sample complexity with respect to $\eta$ is $\mathcal{O}(\eta^{-2})$ for which the optimal scaling is $\mathcal{O}(\eta^{-1/2})$ achieved by QSP and QETU with amplitude amplification.  Methods | Gate complexity | Gate complexity (lattice models) | Extra qubits
---|---|---|---
QPE + Trotter ($2k$th-order) | ${\mathcal{O}}(L\lambda^{1+\frac{1}{2k}}\varepsilon^{-(1+\frac{1}{k})})$ | ${\mathcal{O}}(n^{1+\frac{1}{2k}}\varepsilon^{-(1+\frac{1}{k})})$ | $\log(\varepsilon^{-1})+\log(\Delta^{-1})$
QPE + QW lee2021even  | $\tilde{\mathcal{O}}(L\lambda\varepsilon^{-1})$ | $\tilde{\mathcal{O}}(n^{2}\varepsilon^{-1})$ | $\log(L)+\log(\varepsilon^{-1})+\log(\Delta^{-1})$
QSP lin2020near  | $\tilde{\mathcal{O}}(L\lambda\varepsilon^{-1})$ | $\tilde{\mathcal{O}}(n^{2}\varepsilon^{-1})$ | $\log(L)+\log\log(\varepsilon^{-1})+\log(\Delta^{-1})$
QETU dong2022ground  | $\tilde{\mathcal{O}}(L\lambda^{1+\frac{1}{2k}}\varepsilon^{-(1+\frac{1}{k})})$ | $\tilde{\mathcal{O}}(n^{1+\frac{1}{2k}}\varepsilon^{-(1+\frac{1}{2k})})$ | 1
This work ($2k$th-order) | $\tilde{\mathcal{O}}(L\lambda^{1+\frac{1}{4k+1}}\varepsilon^{-(1+\frac{1}{4k+1})})$ | $\tilde{\mathcal{O}}(n^{1+\frac{2}{4k+1}}\varepsilon^{-(1+\frac{1}{4k+1})})$ | 0 or 1
(zeroth-order) and wan2021randomized  | $\tilde{\mathcal{O}}(\lambda^{2}\varepsilon^{-2})$ | $\tilde{\mathcal{O}}(n^{2}\varepsilon^{-2})$ | 0 or 1
Table 2: Gate complexity in eigenenergy estimation up to estimation error
$\varepsilon$ (Problem 2). Note that this table focuses on the total gate
complexity (i.e., sample complexity is included). The dependence on $\eta$ is
not included in the table, as similarly discussed in Table 1.

###  II.1 Overview of the results

Our result with random composite LCU can nearly match and, in some cases,
outperform the previous best methods for gate complexity with respect to
$\Delta$, $\varepsilon$, $\lambda$ and $n$, though the sample complexity
becomes worse in $\eta$.  In this section, we provide an overview of the main
results of this work. The eigenenergy is first assumed to be known a priori.
This is the case for linear algebra tasks, in which the true solution
corresponds to the ground state of a constructed Hamiltonian with the
eigenenergy being zero. The task with an unknown eigenenergy will be discussed
in subsection B.5. Before diving into the derivations, we first summarise the
main theoretical results for generic Hamiltonians.

######  Theorem 1 (Eigenstate property estimation for generic Hamiltonians).

Suppose we use the method in Algorithm 1 and the conditions and assumptions in
Problem 1 hold. Observable estimation (Problem 1): To achieve the error of
observable’s expectation on the eigenstate $\ket{u_{j}}$ within $\varepsilon$,
the gate complexity in a single circuit is
$\mathcal{O}\left(5^{k-1}L(\lambda\Delta^{-1}\ln(\eta^{-1}\varepsilon^{-1}))^{1+\frac{1}{4k+1}}\right)$
when the number of samples is
$N_{s}=\mathcal{O}\left({\eta^{-2}\varepsilon^{-2}}\|O\|_{1}^{2}\ln(1/\vartheta)\right)$
with a success probability at least $1-\vartheta$.

Eigenenergy estimation (Problem 2): To achieve the eigenenergy estimation
error within $\kappa$ using Algorithm 1, the gate complexity in a single
circuit is
$\mathcal{O}\left(5^{k-1}L(\lambda\kappa^{-1}\ln(\eta^{-1}))^{1+\frac{1}{4k+1}}\right)$,
with number of samples
$N_{s}=\mathcal{O}\left(\eta^{-2}\ln(1/\vartheta)\right)$ with a success
probability at least $1-\vartheta$, approaching to the Heisenberg limit.
Alternatively, by using the methods proposed in wang2023quantum  and Algorithm
1, the gate complexity in a single circuit is
$\mathcal{O}\left(5^{k-1}L(\lambda\Delta^{-1}\ln(\eta^{-1}\kappa^{-1}))^{1+\frac{1}{4k+1}}\right)$
at the cost of more number of samples
$N_{s}=\mathcal{O}\left(\eta^{-2}\Delta^{4}\kappa^{-4}(\ln({\kappa^{-2}\eta^{-1}}))^{2}\ln(1/\vartheta))\right)$
with a success probability at least $1-\vartheta$.

Comparisons of the state-of-the-art ground-state property estimation methods
and ground-state energy estimation are shown in Table 1 and Table 2. We note
that in this work, we focus on the application in NISQ or early FTQC, in which
circuit depth is more concerned than sample complexity. As in random sampling
methods, $\eta$ only appears in sample complexity and thus $\eta$ dependence
is not included in Table 1. As we shall see, the random sampling method has
the advantage in system size scaling when qubit connectivity is restricted.
However, random sampling methods will require more samples much worse than
coherent methods. We leave the comparison of random-sampling methods with
coherent methods in the Discussion.

As shown in Table 1, our method can achieve polylogarithmic dependence on
inverse precision, outperforming the QPE-based method and matching the result
by QSP. QETU can achieve near-optimal ground state preparation with a single
ancillary qubit by querying real-time evolution. However, QETU hinges on a
coherent implementation of $e^{-iHt}$. That means the Hamiltonian simulation
needs to run coherently, ruling out any random sampling method, and thus it is
not straightforward to achieve polylogarithmic dependence on inverse
precision.

An advantage of our method over QSP is that the commutator relations can be
utilised such that it may achieve a better system size dependence. For
Heisenberg models, the gate complexity of our method is
$\mathcal{O}(n^{1+\frac{2}{4k+1}}\Delta^{-(1+\frac{1}{4k+1})}\varepsilon^{-\frac{1}{4k+1}})$.
In contrast, the solution given by QSP lin2020near  is
$\mathcal{O}(n^{2}\Delta^{-1}\log\varepsilon^{-1})$. We summarise the results
for the gate complexity for different physical Hamiltonians. These examples
include (1) spin-lattice Hamiltonians with nearest-neighbour interactions
(e.g., Ising models and Heisenberg models), (2) a second-quantised plane-wave
electronic structure with $n$ spin orbitals and $L=\mathcal{O}(n^{2})$ terms
in the Hamiltonian, which also include the 2D Fermi-Hubbard model, and (3) a
general second-quantised quantum chemistry problem with $L=\mathcal{O}(n^{4})$
terms.

It is worth comparing our method with the Trotter-based methods. The
implementation of time evolution by Trotter methods will undermine the optimal
scaling with respect to $\lambda$, $\Delta$ and $\varepsilon$. Since
controlled operations are necessary in QETU (even for the control-free
method), it becomes suboptimal when considering the qubit connectivity. Our
RLCU method has the following advantages: (1) it has better scaling with
respect to $\lambda$, $\Delta$ and $\varepsilon$, though the dependence on
$\eta$ is much worse than other advanced methods. (2) When restricting the
qubit connectivity to nearest neighbours, our method has improved in circuit
depth. For the Heisenberg model, QETU has depth
$\mathcal{O}(n^{1+\frac{1}{2k}}\varepsilon^{-\frac{1}{2k}})$, while our method
has depth $\mathcal{O}(n^{\frac{2}{4k+1}}\varepsilon^{-\frac{1}{4k+1}})$,
which shows advantages in the scaling of $n$ and $\varepsilon$.

The result is summarised in Theorem 2 and a comparison with other methods is
shown in the third column of Table 1.

######  Theorem 2 (Eigenstate property estimation for lattice Hamiltonians).

For $n$-qubit Heisenberg Hamiltonians, to estimate the observable on the
eigenstate with a precision $\varepsilon$ and a success probability
$1-\vartheta$, in a single run, the gate complexity is
$\mathcal{O}(n^{1+\frac{2}{4k+1}}\Delta^{-(1+\frac{1}{4k+1})}\varepsilon^{-\frac{1}{4k+1}}\log(\vartheta^{-1}))$.
The circuit depth when compiled on qubits with nearest-neighbour geometry is
$\mathcal{O}(n^{\frac{2}{4k+1}}\Delta^{-(1+\frac{1}{4k+1})}\varepsilon^{-\frac{1}{4k+1}}\log(\vartheta^{-1}))$.

The proof idea for Theorem 1 and Theorem 2 is illustrated in Fig. 2 in Sec.
III. The proof and a more detailed version of Theorem 1 can be found in
Theorem 3 in Appendix C.

Hamiltonians | Ancilla-free method (NN) | QETU (NN) dong2022ground  | arbitrary
---|---|---|---
1D lattice models | $d=\mathcal{O}(\mathbf{n^{\frac{2}{4k+1}}}\Delta^{-(1+\frac{1}{4k+1})}\varepsilon^{-\frac{1}{4k+1}})$ | $d=\mathcal{O}(n^{1+\frac{1}{2k}}\Delta^{-(1+\frac{1}{2k})}\varepsilon^{-\frac{1}{2k}})$ | $d=\mathcal{O}(n^{\frac{2}{4k+1}}\Delta^{-(1+\frac{1}{4k+1})}\varepsilon^{-\frac{1}{4k+1}})$
| $g=\mathcal{O}(n^{1+\frac{2}{4k+1}}\Delta^{-(1+\frac{1}{4k+1})}\varepsilon^{-\frac{1}{4k+1}})$ | $g=\mathcal{O}(n^{1+\frac{1}{2k}}\Delta^{-(1+\frac{1}{2k})}\varepsilon^{-\frac{1}{2k}})$ | $g=\mathcal{O}(n^{1+\frac{2}{4k+1}}\Delta^{-(1+\frac{1}{4k+1})}\varepsilon^{-\frac{1}{4k+1}})$
Electronic structure | $d=\mathcal{O}(\mathbf{n^{2+\frac{2}{4k+1}}}\Delta^{-(1+\frac{1}{4k+1})}\varepsilon^{-\frac{1}{4k+1}})$ | $d=\Omega(n^{3+\frac{1}{2k}}\Delta^{-(1+\frac{1}{2k})}\varepsilon^{-\frac{1}{2k}})$ | $d=\mathcal{O}(n^{2+\frac{2}{4k+1}}\Delta^{-(1+\frac{1}{4k+1})}\varepsilon^{-\frac{1}{4k+1}})$
(Eq. 6) | $g=\mathcal{O}(n^{3+\frac{2}{4k+1}}\Delta^{-(1+\frac{1}{4k+1})}\varepsilon^{-\frac{1}{4k+1}})$ | $g=\mathcal{O}(n^{3+\frac{1}{2k}}\Delta^{-(1+\frac{1}{2k})}\varepsilon^{-\frac{1}{2k}})$ | $g=\tilde{\mathcal{O}}(n^{2+\frac{2}{4k+1}}\Delta^{-(1+\frac{1}{4k+1})}\varepsilon^{-\frac{1}{4k+1}})$
Table 3: Gate complexity with respect to the maximum simulation time $t$,
target precision $\varepsilon$, and system size $n$ for different
Hamiltonians. In the second and third columns, the qubit connectivity is
restricted to a linear nearest-neighbour (NN) architecture. The fourth column
is the result of our method when there is no restriction on qubit
connectivity. In this table, the commutation relation of the Hamiltonian terms
is used to improve the system-size scaling. It is worth mentioning that the
gate complexity of electronic structure problems studied in this table can be
$g=\tilde{\mathcal{O}}(n(\lambda\Delta^{-1}\log\varepsilon^{-1}))^{1+\frac{1}{4k+1}}$
which is nearly linear in $n$ and logarithmic in inverse precision.

The result can be extended to the simulation of molecules. In the plane-wave
dual basis, the second quantised electronic-structure Hamiltonian has the form
(see babbush2018encoding ; kivlichan2020improved ; childs2021theory )

$H=\hat{T}+\hat{V}=\sum_{pq}T_{pq}\hat{a}_{p}^{\dagger}a_{q}+\sum_{p}U_{p}\hat{n}_{p}+\sum_{p\neq
q}V_{pq}\hat{n}_{p}\hat{n}_{q}$ (6)

where $\hat{T}$ and $\hat{V}$ represent the kinetic and potential terms of the
fermionic Hamiltonian, respectively, $\hat{a}$ and $\hat{a}^{\dagger}$ are
fermionic creation and annihilation operators and $\hat{n}_{p}$ is the number
operator for the corresponding spin-orbital, and the total number of terms
$L=\mathcal{O}(n^{2})$. Note that operators are typically indicated by using
the hat notation. In this work, to avoid confusion with the estimator, we will
only use the hat symbol to denote the fermionic operators and the estimators.
Except for the potential terms $\hat{V}$ in molecular Hamiltonians, $V$
denotes the Trotter remainder in this work.

To estimate the eigenstate property of the Hamiltonian described by Eq. 6, the
circuit depth scales as $\mathcal{O}(n^{2+\frac{2}{4k+1}})$. When restricted
to nearest-neighbour (NN) architecture, the circuit depth by QETU scales as
$\Omega(n^{3+\frac{1}{2k}})$. The circuit depth results for different physical
Hamiltonians are displayed in Table 3. Here, we note that QSP and methods
based on full-order pairing wan2021randomized  have worse system-size scaling,
and thus they are not included in Table 3.

The primary contributions of this work are summarised as follows.

  1. 1.

We provide end-to-end gate complexity analysis of the eigenstate property and
energy estimation task, and can achieve nearly logarithmic scaling on inverse
precision, with improved scaling on $\Delta$ and $\lambda$ compared to oracle-
free methods. The comparison with other methods is shown in Table 1 and Table
2. For electronic structure problems specified in Eq. 6 which is usually
compared in the existing literature, the gate complexity for eigenstate
property estimation and eigenenergy estimation can be
$\tilde{\mathcal{O}}(n(\lambda\Delta^{-1}\log\varepsilon^{-1})^{1+\frac{1}{4k+1}})$,
which is nearly linear in $n$ (excluding the $\log n$ factor) and logarithmic
in inverse precision. Note that compared to reiher2017elucidating ;
lee2021even ; babbush2018encoding ; babbush2018low , this result is for the
gate count in a single-run circuit with much worse sample complexity.

  2. 2.

The random sampling approach can exploit the commutation relation of the
target Hamiltonian terms to reduce the gate complexity, outperforming QSP for
the 1D lattice model. Moreover, our approach favours a linear nearest-
neighbour (NN) architecture. When restricted to nearest-neighbour
architecture, for the 1D lattice model, the circuit depth scales
$d=\mathcal{O}(n^{\frac{2}{4k+1}})$, while the circuit depth in the QETU-
Trotter method scales as $d_{\mathrm{QETU}}=\mathcal{O}(n^{1+\frac{1}{2k}})$.
For electronic structure problems specified in Eq. 6, when considering the
commutation relation and restricting to a NN architecture, the circuit depth
is $d=\mathcal{O}(n^{2+\frac{2}{4k+1}})$, while naive implementation of QETU +
Trotter requires $d=\Omega(n^{3+\frac{1}{2k}})$. A side product is that
controlled $e^{-i\theta{H}}$ can be implemented by a linear depth circuit
$d=\mathcal{O}(n)$ summarised in Proposition 4, comparable to the control-free
simulation of electronic problems in kivlichan2018quantum ; babbush2018low .
This result could be used as a subroutine in other quantum algorithms.

  3. 3.

We analyse the actual gate count for typical problems and compare it with the
state-of-the-art methods (see subsection III.4 for the analytic results and
Sec. IV for numerical results).

Figure 2:  The diagram illustrates the hierarchy of various components within
the composite LCU. The figure shows how different sub-sections in Sec. III
connect. This visualisation aids in understanding how the target problem
(Problem 1 and Problem 2) can be achieved, specifically, to achieve
$(\varepsilon,\vartheta)$-eigenstate property estimation. Arrows in the figure
indicate direct connections between components. The key parameters relevant to
the construction of LCU are highlighted.

##  III Methods

In subsection III.1, we first construct the LCU formula for the spectral
filter outlined in Fig. 1 and Algorithm 1. The hierarchy of different
components in the composite LCU is shown in Fig. 2. subsection III.2 shows how
errors in observable expectations can be bounded given the approximation error
of the composite LCU formula. subsection III.3 shows the circuit realisation
of the composite LCU form by presenting an explicit form of each circuit
instance. An ancilla-free composite LCU formula and measurement scheme is
presented for symmetry-conserved cases. subsection III.4 analyses the circuit
depth and gate complexity for general cases and Hamiltonian-specific cases.
Fig. 2 illustrates the proof idea of Theorem 1 and Theorem 2, and shows the
connections of different subsections.

###  III.1 Composite LCU form of the spectral filter

Now, we show how to construct the composite LCU form of the spectral filter.
We mainly consider the Gaussian filter $g_{\tau}(H)=g(\tau H)$ but present the
construction of $g$ in a general form, and hence the results can be readily
applied to the investigation of other functions. For the target Hamiltonian
$H$, suppose we want to prepare the $j$th eigenstate $\ket{u_{j}}$ with
eigenenergy $E_{j}$. We consider the spectral filter $g(\tau(H-\omega))$ with
$\omega=E_{j}$, which can be expressed as follows,

$g(\tau(H-\omega))=\sum_{i=0}^{N-1}g(\tau(E_{i}-\omega))\ket{u_{i}}\bra{u_{i}}.$
(7)

For an input state $\ket{\psi_{0}}=\sum_{i}c_{i}\ket{u_{i}}$, the state after
applying the spectral filter $g(\tau(H-\omega))$ at a finite $\tau$ is given
by

$\ket{\psi(\tau)}=\frac{g(\tau(H-\omega))\ket{\psi_{0}}}{\|g(\tau(H-\omega))\ket{\psi_{0}}\|}=\frac{\sum_{i}g(\tau(E_{i}-\omega))c_{i}\ket{u_{i}}}{\sum_{i}|c_{i}|^{2}g(\tau(E_{i}-\omega))^{2}}.$
(8)

It is easy to see that the Gaussian function $g(\tau(E_{i}-\omega))$ decreases
exponentially with $\tau$ and $E_{i}-\omega$. When taking $\omega=E_{j}$, the
amplitudes of the normalised state $\ket{\psi(\tau)}$ concentrate to the
eigenstates with energy $E_{j}$, and the evolved state asymptotically
approximates the eigenstate $\ket{u_{j}}$ with nonzero
$|\braket{\psi_{0}}{u_{j}}|^{2}\neq 0$ for sufficiently large $\tau$ as
$\lim_{\tau\to\infty}g(\tau(H-E_{j}))\ket{\psi_{0}}\propto\ket{u_{j}}.$

There are two important parameters $\tau$ and $\omega$ in the spectral filter
$g(\cdot)$, in which $\tau$ is a time scaling factor, and $\omega=E_{j}$
indicates a shift of the original function. We shall see that $\tau$ indicates
the timescale for the spectral filtering procedure, and larger $\tau$ will
cool the state closer to the target eigenstate. The shift $\omega$ plays an
important role in searching the eigenenergies (see Fig. 1(b)) and also the
eigenstate property estimation. From the above equation, readers may wonder if
we can still get the concentration around the eigenstate $\ket{u_{j}}$ with
energy $\omega$ if we do not know the value of $E_{j}$ a priori. Note that
$\omega$ only appears in classical post-processing and will not be involved in
quantum measurement. In the implementation, we will classically tune $\omega$
as a classical variable to find the true eigenenergy while it will not
increase any quantum resource cost.

The Gaussian spectral filter can be decomposed into the basis of real-time
evolution as

$\displaystyle
g(\tau(H-\omega))=c\int^{\infty}_{-\infty}\textrm{d}p(x)e^{ix(\omega-H)}$ (9)

with $p(x)$ being a Gaussian distribution. The quantum state given by Eq. 8
becomes

$\displaystyle\ket{\psi(\tau)}\propto
c\int^{\infty}_{-\infty}\textrm{d}p(x)e^{i\tau x\omega}\ket{\phi(x\tau)},$
(10)

which is now a superposition of real-time evolved states
$\ket{\phi(x\tau)}=e^{-ix\tau H}\ket{\psi_{0}}$ with probability distribution
$\textrm{d}p(x)=p(x)\textrm{d}x$.

Instead of preparing the quantum state of Eq. 10, we focus on the goal of
obtaining arbitrary observable expectation values. Specifically, we aim to
measure any observable $O$ of the evolved state $\ket{\psi(\tau)}$, i.e.,

$\braket{O}_{\psi(\tau)}=\braket{\psi(\tau)}{O}{\psi(\tau)}=\frac{N_{\tau}(O)}{D_{\tau}(\omega)},$
(11)

where

$\displaystyle D_{\tau}(\omega)$
$\displaystyle=\braket{\psi_{0}}{g^{2}(\tau(H-\omega))}{\psi_{0}},$ (12)
$\displaystyle N_{\tau}(O)$
$\displaystyle=\braket{\psi_{0}}{g(\tau(H-\omega))Og(\tau(H-\omega))}{\psi_{0}}$

and we have omitted $\omega$ in $N_{\tau}$ to simplify the notation. The
denominator is

$\displaystyle D_{\tau}(\omega)$
$\displaystyle=c^{2}\int_{-\infty}^{\infty}\tilde{p}(x)e^{i\tau
x\omega}\braket{\psi_{0}}{e^{-i\tau xH}}{\psi_{0}}$ (13)

with $\textrm{d}\tilde{p}(x)=\tilde{p}(x)\textrm{d}y$ and
$\tilde{p}(x)=\frac{1}{2}\int^{\infty}_{-\infty}p(\frac{z+x}{2})p(\frac{z-x}{2})\textrm{d}z$.
The expectation values of the denominator and the numerator can be expressed
as $v=\operatorname{Tr}(Og_{\tau}(H-\omega)\rho_{0}g_{\tau}(H-\omega))$, and
the denominator is by taking $O=I$.

For the numerator $N_{\tau}(O)$, we can efficiently obtain it by sampling the
distribution $\textrm{d}p(x,x^{\prime})$ and then estimating the mean value
$\mathbb{E}_{x,x^{\prime}}\braket{\phi(x^{\prime}\tau)}{O}{\phi(x\tau)}$,
where each term can be measured by the Hadamard test circuit or using the
ancilla-free measurement strategy detailed in subsection III.3. We can
similarly obtain the denominator $D_{\tau}$ by estimating
$\mathbb{E}_{x}\braket{\psi_{0}}{e^{-ix\tau H}}{\psi_{0}}$ with probability
$\textrm{d}\tilde{p}(x)$. Therefore, we only need a quantum computer to
efficiently estimate the values of state overlaps like
$\braket{\psi_{0}}{U}{\psi_{0}}$ with $U$ being either $e^{ix^{\prime}\tau
H}Oe^{-ix\tau H}$ or $e^{-ix\tau H}$. In practice, we also consider a cutoff
of the integral from $[-\infty,\infty]$ to $[-x_{c},x_{c}]$ to avoid infinite
integration.

The unnormalised eigenstate can be effectively realised by applying a spectral
filter $g_{\tau\rightarrow\infty}(H-E_{j})$ to an initial state, which holds
when the imaginary-time becomes infinity. When $\tau\rightarrow\infty$, one
could check that

$\braket{O}=\frac{N_{\tau\rightarrow\infty}(O)}{D_{\tau\rightarrow\infty}}$

and thus we arrive at the ideal observable expectation. Note that the
denominator is nonvanishing $D_{\tau\rightarrow\infty}(E_{j})=\eta$ under the
assumption. The observable when considering a finite $\tau$ is estimated by

$\braket{O}_{\tau}=\frac{N_{\tau}(O)}{D_{\tau}}.$ (14)

Furthermore, to realise the spectral filter with unitary operations, we need
to avoid the infinite time length and truncate the time in real-time evolution
to a given threshold. Below we discuss the LCU form of $g$ when considering
the cutoff.

The Gaussian spectral filter takes an explicit composite LCU form of

$\displaystyle g_{\tau}(H-\omega)$
$\displaystyle=c\int_{-x_{c}}^{+x_{c}}dxp(x)e^{ix\tau\omega}e^{-i\tau xH}$
(15)

with truncation $[-x_{c},x_{c}]$. The integrand is a real-time evolution with
total real-time $\tau x$. Suppose that we further use LCU formula to implement
$e^{-i\tau xH}$, which takes the form of

$e^{-i\tau
xH}=\mu(x\tau)\sum_{\vec{r}\in\mathcal{K}_{x}}\Pr(\vec{r},x\tau,\nu(x\tau))U_{\vec{r}}$
(16)

where $\vec{r}$ specifies the unitary $U_{\vec{r}}$ involved in the LCU
formula of $e^{-i\tau xH}$, $U_{\vec{r}}$ is a unitary operator, and
$\Pr(\vec{r},x\tau,\nu(x\tau))$ represents the normalised decomposition
coefficients of $U_{\vec{r}}$. Then, we have

$\displaystyle g_{\tau}(H-\omega)$ (17)
$\displaystyle=c(\mu)\int_{-x_{c}}^{+x_{c}}dxp_{\mu}(x)e^{i\theta_{x}}e^{ix\tau\omega}\sum_{\vec{r}\in\mathcal{K}_{x}}\Pr(\vec{r},x\tau,\nu(x\tau))U_{\vec{r}}$

where we defined $c(\mu):=c\int_{-x_{c}}^{+x_{c}}p(x)\mu(x\tau)dx$ and
$p_{\mu}(x)=p(x)\mu(x\tau)/c(\mu)$ in a similar vein to Proposition 5.

The central question is how to find the operator $U_{\vec{r}}$ and the
corresponding probability distribution, with which the LCU decomposition in
Eq. 16 can be determined. The overall idea is to divide the time evolution
into $\nu$ segments, and for each segment, we realise both the Trotter
evolution and the Trotter error, as proposed in zeng2022simple . For evolution
time $t$, let us denote the real-time evolution $U(t):=e^{-iHt}=U_{m}^{\nu}$
where $U(m)=e^{imH}$ could be expressed as

$\displaystyle{U}(m)={V}_{2k}(m)S_{2k}(m),$ (18)

which consists of a deterministic $2k$th-order Trotter formula $S_{2k}$ and
the Trotter remainder ${V}_{2k}$. The LCU formula of ${V}_{2k}$ can be
expressed as

$\tilde{V}_{2k}(m)=\mu(m)\sum_{r}\Pr(r,m,\nu(m))W_{r}$ (19)

where $W_{k}$ is a unitary operator with the error of the formula
$\varepsilon_{2k}(m)$. The overall LCU formula of $U(t)$ is a product of each
individual formula

$U(t)=(U(m))^{\nu}=\mu(t)\sum_{\vec{r}\in\mathcal{K}_{t}}\Pr(\vec{r},t,\nu(t))U_{\vec{r}}$
(20)

where $\vec{r}:=(r_{1},r_{2},...,r_{\nu})$ with each $U_{r_{i}}$ being sampled
from the distribution of $\Pr(r,m,\nu(m))$ introduced in Eq. 19 and we denote
$\Pr(\vec{r},t,\nu(t)):=\prod_{i\leq\nu}\Pr(r_{i},m,\nu(m))$ and
$U_{\vec{r}}:=\prod_{i\leq\nu}W_{r_{i}}S_{2k}$. With some derivation, one can
prove that $\mu=\mu(t)^{\nu}$ and
$\varepsilon_{2k}(t)\leq\nu\mu(m)\varepsilon_{2k}(m)$. Eq. 20 shows how
$U_{\vec{r}}$ can be sampled and thus how $g$ can be realised. This is
illustrated in Fig. 1(a3); in the figure, the index by $i$ is related to the
probability $\Pr(i)$ when decomposing $g$ into real-time evolution $U(t_{i})$.

Using BCH formula, the $2k$th-order remainder
$V_{2k}(m):=U(m)S_{2k}(m)^{\dagger}$ has an explicit form of $V_{2k}(m)$ as
$V_{2k}(m)=\exp(i\sum_{s=1}^{\infty}E_{2k,s}m^{s})$ with Hermitian operators
$E_{2k,s}$. Since the $2k$th-order Trotter error is compensated, we know that
$E_{2k,s}=0$ for $s\leq 2k$. Expanding $V_{2k}(m)$, we have
$V_{2k}(m)=\sum_{s=0}^{\infty}F_{2k,s}(m)$ where we group the terms by the
order of $m$, and $F_{1,s}$ denotes the $s$-order expansion term of
$V_{2k}(m)$ with $m^{s}$.

Below, we briefly discuss the error in LCU form when considering truncation of
$s$. The truncation error is found to have a quick decrease with an increasing
truncation order $s_{c}$ zeng2022simple . Considering of a finite truncation
and the vanishing of $F_{2k,s\leq 2k}(m)=0$, we can rewrite $V_{2k}^{(s)}(m)$
as

$V_{2k}^{(s)}(m)=I+\sum_{s=2k+1}^{s_{c}}F_{2k,s}(m).$ (21)

Given the truncation, the LCU formula for $U(m)$ is

$\displaystyle{U}_{2k}^{(s_{c})}(m)={V}_{2k,s}(m)S_{2k}(m),$ (22)

which consists of a deterministic second-order Trotter formula and the Trotter
error compensation term. The overall LCU formula for $U(t)$ is to repeat the
sampling of $\tilde{U}_{2k}(m)$ for $\nu$ times,
${U}_{2k}^{(s_{c})}(t)=\left({U}_{2k}^{(s_{c})}(m)\right)^{\nu}$ with
$\nu=t/m$.

Eq. 17 is a composite LCU formula; more specifically, it is a $(c(\mu),0)$ LCU
formula. When we consider a finite $x_{c},s_{c}$, it will introduce some
errors. When considering a finite $s_{c}$, the spectral filter becomes

$g_{\tau,x_{c},s_{c}}(H-\omega)=c\int_{-x_{c}}^{x_{c}}dxp(x)e^{ix\tau\omega}{U}_{2k}^{(s_{c})}(x\tau,\nu(x\tau))$
(23)

It is easy to check that $\|g_{\tau,x_{c},s_{c}}(H-\omega)\|\leq
c(1+\varepsilon_{s_{c}})$. The operator distance between $g_{\tau,x_{c}}$ and
$g_{\tau,x_{c},s_{c}}$ due to finite $s_{c}$ is

$\displaystyle\|g_{\tau,x_{c}}-g_{\tau,x_{c},s_{c}}\|$ (24) $\displaystyle\leq
c\int_{-x_{c}}^{x_{c}}dxp(x)e^{ix\tau\omega}\|U(x\tau)-{U}_{2k}^{(s_{c})}(x\tau,\nu(x\tau))\|\leq\varepsilon_{s_{c}}$

when $\|U-{U}_{2k}^{(s_{c})}(x\tau,\nu(x\tau))\|\leq\varepsilon_{s_{c}}/c$. In
order to achieve an additive error of the approximation

$\|U(t)-{U}_{2k}^{(s_{c})}(t))\|\leq\varepsilon_{s_{c}}$ (25)

Lemma 1 in the Appendix gives the required segment numbers $\nu$ given the
truncation order $s_{c}$.

With the above result, the total error of the spectral filter $g$ can be
bounded by each individual error term. This process is illustrated within the
orange box of Fig. 2. Specifically, we can bound the error in spectral filter
construction by comparing the operator distance between

$\|g_{\tau,x_{c},s_{c}}-g_{\tau,x_{c},s_{c}\rightarrow\infty}\|,$ (26)

which gives an upper bound for the numerator and the denominator. Using the
triangular inequality, the operator error between $g_{\tau}$ and
$g_{\tau,x_{c},s_{c}}$ defined in Eq. 95 is

$\|g_{\tau}-g_{\tau,x_{c},s_{c}}\|\leq\varepsilon_{x_{c}}+\varepsilon_{s_{c}}.$
(27)

In Theorem 3 in Appendix B, we show that when the segment number is of the
order

$\nu=\mathcal{O}((\lambda\tau x_{c})^{1+\frac{1}{4k+1}})$ (28)

it suffices to ensure the LCU construction error up to $\varepsilon$.
Similarly, Theorem 5 in Appendix Appendix B shows the required segment number
for eigenenergy estimation. Together with Proposition 3 in Appendix subsection
III.4, one can prove Theorem 1. For the lattice model’s eigenstate property
estimation, the segment number can be further reduced to
$\nu=\mathcal{O}(n^{\frac{2}{4k+1}}(\tau
x_{c})^{1+\frac{1}{4k+1}}\varepsilon^{-\frac{1}{4k+1}})$. With this result,
one can prove Theorem 2.

###  III.2 Observable estimation

Here, we discuss the eigenstate property estimation when the state is filtered
under the constructed LCU formula discussed above. We first discuss how to
measure the expectation value of observable $O$. Denote the initial state on
which the spectral operator acts as $\rho$. The expectation value of
observable $O$ on the normalised state can be expressed by

$\braket{O}=\frac{N(O)}{D}=\frac{\operatorname{Tr}(g\rho
g^{\dagger}O)}{\operatorname{Tr}(g\rho g^{\dagger})}.$ (29)

We have

$N(O)=\mu^{2}\sum_{ij}\Pr(i)\Pr(j)e^{i(t_{i}-t_{j})\omega}\operatorname{Tr}(U(t_{i})\rho
U^{\dagger}(t_{j})O).$ (30)

Since $U(t_{i})$, $U(t_{j})$ are independent, readers may group it and take
the symmetric estimator as $\operatorname{Tr}(U(t_{i})\rho
U(t_{j})^{\dagger}O)+\operatorname{Tr}(U(t_{i})\rho U(t_{j})^{\dagger}O)$ to
save measurement resources. However, this symmetric estimator is not
applicable to dealing with general cases. In the later discussion, we shall
see that there is an additional phase factor associated with $U(t_{i})$ such
that the imaginary part of $\operatorname{Tr}(U(t_{i})\rho
U(t_{j})^{\dagger}O)$ matter for the final result. For the sake of
consistency, we will focus on the discussion on the asymmetric estimator.

Denote the expectation value of the estimator of the numerator over
measurement outcomes with $U(t_{i})$ and $U(t_{j})$,
$\operatorname{Tr}(U(t_{i})\rho U(t_{j})^{\dagger}O)$, as $\bar{N}_{ij}(O)$.
The expectation value of $N(O)$ is given by sampling over the distribution

$\bar{N}(O)=\mu^{2}\mathbb{E}_{ij}\bar{N}_{ij}(O)$ (31)

As random sampling is used in our method, we cannot prepare the state
deterministically. However, in terms of the expectation of observables, our
scheme owns a similar function to the other schemes.

When considering finite gate complexity and sample complexity, the eigenstate
property is estimated by

$\hat{O}_{\tau,x_{c},s_{c}}=\frac{\hat{N}_{\tau,x_{c},s_{c}}(O)}{\hat{D}_{\tau,x_{c},s_{c}}}.$
(32)

The selection of $\tau,x_{c},s_{c}$ can be determined by analysing the error
of $\hat{O}_{\tau,x_{c},s_{c}}$ compared to the ideal value, which can be
analysed using the operator distance detailed in the above section.
Proposition 2 shows the error in observable estimation when we are given a
$(\mu,\varepsilon)$ random sampling formula of $g$. The proof is shown in the
Appendix.

######  Proposition 2 (Observable estimation using the composite LCU formula).

For a target operator $g$ and its $(\mu,\varepsilon)$-random-sampling formula
defined in Eq. 1, if we estimate the value on the unnormalised state
$N_{g}(O):=\operatorname{Tr}(g\rho g^{\dagger}O)$ with an initial state $\rho$
and observable $O$, then the distance between the mean estimator value
$\hat{O}$ and the true value $N_{g}(O)$ is bounded by

$\varepsilon_{N}:=|\hat{N}_{\tilde{g}}(O)-N_{g}(O)|\leq\|O\|(2\mu^{2}\varepsilon+\varepsilon_{n}),$
(33)

with a success probability $1-\vartheta$. Here, we use the variant of the
Hadamard test circuit for $N=\mu^{4}\ln(2/\vartheta)/\varepsilon_{n}^{2}$
times and $\|O\|$ is the spectral norm of $O$. The error for the denominator
is bounded as $\varepsilon_{D}:=|\hat{D}_{\tilde{g}}-D_{g}|\leq
2\mu^{2}\varepsilon+\varepsilon_{n}$. Given the error $\varepsilon_{D}$ and
$\varepsilon_{N}$, the error for the normalised observable expectation
$\braket{O}_{g}=N_{g}(O)/D_{g}$ is

$\begin{aligned}
|\frac{\hat{N}_{\tilde{g}}(O)}{\hat{D}_{\tilde{g}}}-\frac{N_{g}(O)}{D_{g}}|\leq\frac{1}{D_{g}}((\braket{O}+1)\varepsilon_{D}+\varepsilon_{N})\end{aligned}.$
(34)

Given Proposition 2, the error of the observables can be bounded given the
RLCU form. We remark that the random sampling algorithm is not a deterministic
state preparation that cannot prepare the target state $\sigma$. However, when
we focus on the property estimation of the target state, the estimator owns
similar performance to the normal Hamiltonian simulation methods: the sample
complexities of the former algorithms and RCLU algorithm to learn the
observable properties of the state are similar, as long as the norm of the
estimator is a constant. This is guaranteed by Proposition 2.

Using these results of error propagation in the LCU form, we can determine the
required segment number $\nu_{c}$ in order to achieve the desired accuracy as
shown in Eq. 28. Theorem 3 in Appendix Appendix B and Theorem 5 in Appendix
Appendix C present the required segment number $\nu_{c}$ and the gate count
(with the prefactors included) for achieving desired accuracy for two tasks,
which will be used to estimate the actual gate numbers needed for physical
Hamiltonians.

The total gate count is determined by multiplying $\nu_{c}$ and the gate count
for realising $S_{2k}$ and $V_{2k}$ within each segment, denoted by $g$
hereafter. The gate count for lattice models is straightforward
$g=\mathcal{O}(n)$. For electronic structure problem in Eq. 6, a naive
realisation of $S_{2k}$ and $V_{2k}$ will result in the cost
$g=\mathcal{O}(n^{3})$. However, we could bring down the scaling by
considering the properties of fermions. The kinetic operator is quadratic and
thus can be diagonalised by an efficient circuit transformation, either by
Givens rotations or fermionic fast Fourier transform (FFFT)
kivlichan2020improved ; kivlichan2018quantum ; babbush2018low . By using the
results in low2018hamiltonian , $e^{-im\hat{T}}$ and $e^{-im\hat{V}}$ can be
implemented with $\mathcal{O}(n\log n)$ gates with $m$ being the time length.
The compensation term will cost $\mathcal{O}(n)$ gates. Therefore, we have
$g=\mathcal{O}(n\log n)$ and the total gate count for electronic structure
problem is $g_{\rm
tot}=\tilde{\mathcal{O}}(n(\lambda\Delta^{-1})^{1+\frac{1}{4k+1}}\log\varepsilon^{-1})$.

###  III.3 Estimation by quantum circuits: general strategy and ancilla-free
strategy

In this section, we discuss the circuit realisation for Algorithm 1. For
symmetry-conserved cases, we introduce an ancilla-free LCU construction
strategy and the corresponding measurement strategy. To estimate an
observable, a quantum computer is only used for estimating the quantity
involving real-time evolution $\braket{U^{\dagger}(t_{j})OU(t_{i})},$ where
$U({t_{i}})=e^{-iHt_{i}}$. As that $U(t_{i})$ is implemented by a Trotter-LCU
formula in Eq. 18, involving the Trotter operator $S_{2k}$ and the Trotter
error compensation term $V_{2k}^{(s)}$, the unitary can be expressed in an
explicit LCU form as

$U(t_{i})=\mu\sum_{\vec{i}}\Pr(\vec{i})\prod_{q=1}^{\nu}W_{i_{q}}S$ (35)

where $\vec{i}=(i_{1},i_{2},...,i_{\nu})$. The observable expectation is

$\braket{U^{\dagger}(t_{j})OU(t_{i})}=\mu^{2}\sum_{\vec{i},\vec{j}}\Pr(\vec{i})\Pr(\vec{j})\braket{\vec{\mathbf{j}}}{O}{\vec{\mathbf{i}}}$
(36)

where we denote
$\ket{\vec{\mathbf{i}}}:=\prod_{q=1}^{\nu}W_{i_{q}}S\ket{\psi_{0}}$.

To measure $\braket{\vec{\mathbf{j}}}{O}{\vec{\mathbf{i}}}$, we can use the
Hadamard test circuit with an extra ancillary qubit. For Hamiltonians with
certain symmetries, we do not need this ancillary qubit and thus do not
require any controlled operation. Below, we first introduce the general cases
by introducing an ancillary qubit. Then we discuss the ancilla-free composite
LCU formula and the corresponding measurement schemes.

####  III.3.1 General cases

There are two extreme cases: Case I $t_{i}=0$ (or $t_{j}=0$) and Case II
$t_{i}=t_{j}$. For $t_{j}=0$, each component in Eq. 36
$\braket{\psi_{0}}{O}{\vec{\mathbf{i}}}$ can be measured using a Hadamard test
circuit, which is the circuit depicted in the green box in Fig. 1(c).

For $0<t_{j}<t_{i}$, suppose that $U(t_{j})$ is divided into $m$ segments and
$U(t_{i})$ is divided into $\nu$ segments with the first $m$ segments being
set the same as that of $U(t_{i})$. The quantum circuit that can measure
$\braket{\vec{\mathbf{j}}}{O}{\vec{\mathbf{i}}}$ is shown in Fig. 1. After
post-selecting $\ket{+}$ on the ancillary qubit (but before measurements on
$B$), the circuit on system $B$ outputs

$\frac{1}{2}\left(\prod_{q=1}^{m}W_{i_{q}}S+\prod_{q=1}^{\nu}W_{i_{q}}S\right)\ket{\psi_{0}}$
(37)

The case of $t_{i}=t_{j}$ is a special case in Fig. 1(c1) by only implementing
the circuits in the orange box.

Here, we can use a single-shot measurement strategy to estimate $N(O)$ in Eq.
30. Suppose that the sequence of the unitary operators has been obtained from
sampling. To measure $\braket{\vec{\mathbf{j}}}{O}{\vec{\mathbf{i}}}$, we can
use the circuit in Fig. 1(c1). We first prepare the state $\ket{\psi_{0}}$ and
an extra ancillary qubit prepared on $|+\rangle$. Afterwards, we perform a C-U
gate from ancillary to $\ket{\psi_{0}}$. If we directly measure the ancillary
qubit on the X-basis, the outcome $a$ will be $0$ with a probability of
$\Pr(a=0)=\frac{1}{2}(1+\operatorname{Re}(\braket{\vec{\mathbf{j}}}{O}{\vec{\mathbf{i}}})$
and 1 with a probability of
$\Pr(a=1)=\frac{1}{2}(1-\operatorname{Re}(\braket{\vec{\mathbf{j}}}{O}{\vec{\mathbf{i}}})$.
Then we add a phase gate $S={R_{z}}(\frac{\pi}{2})$ before measurement. The
outcome $b$ will be $0$ with a probability of
$\Pr(b=0)=\frac{1}{2}(1+\operatorname{Im}(\braket{\vec{\mathbf{j}}}{O}{\vec{\mathbf{i}}})$
and 1 with a probability of
$\Pr(b=1)=\frac{1}{2}(1-\operatorname{Im}(\braket{\vec{\mathbf{j}}}{O}{\vec{\mathbf{i}}})$.

We consider the estimator as

$\hat{d}=(-1)^{a}+i(-1)^{b}$ (38)

based on the measurement outcome $a$ and $b$ obtained from the above circuit.
$\hat{d}$ is an unbiased estimator of
$\braket{\vec{\mathbf{j}}}{O}{\vec{\mathbf{i}}}$. One can verify that

$\mathbb{E}_{a}(-1)^{a}=\operatorname{Re}(\braket{\vec{\mathbf{j}}}{O}{\vec{\mathbf{i}}}),~{}\mathbb{E}_{b}(-1)^{b}=\operatorname{Im}(\braket{\vec{\mathbf{j}}}{O}{\vec{\mathbf{i}}})$
(39)

and thus
$\mathbb{E}_{a,b}\hat{d}=\braket{\vec{\mathbf{j}}}{O}{\vec{\mathbf{i}}}$.

We take an estimator as

$\hat{v}=c^{2}(\mu)e^{i\omega(t_{i}-t_{j})}\hat{d}$ (40)

The estimator $\hat{v}$ defined in Eq. 40 is an unbiased estimator of
$\braket{\psi_{0}}{g_{\tau}(H-\omega)Og_{\tau}(H-\omega)\psi_{0}}$,

$\mathbb{E}_{t_{i},t_{j},\vec{i},\vec{j}}\mathbb{E}_{a,b}\hat{v}=\braket{\psi_{0}}{g_{\tau}(H-\omega)Og_{\tau}(H-\omega)\psi_{0}}.$
(41)

By analysing the error dependence on the finite $\tau$, $x_{c}$, segment
number $\nu$ in Eq. 28, we arrive at the result of property estimation in
Theorem 1 (see detailed description in Theorem 3 and Theorem 5). More detailed
proof can be found in Appendix Appendix B and Appendix C.

####  III.3.2 Ancilla-free composite LCU formula and the corresponding
measurement scheme for symmetry conserved systems

In the above section, given a sampled configuration
$(t_{i},t_{j},\vec{i},\vec{j})$, the real and imaginary part of
$\braket{\vec{\mathbf{j}}}{O}{\vec{\mathbf{i}}}$ can be obtained by the
circuit in Fig. 1(c). In cases of the Heisenberg model and electronic
structure problems (Eq. 6), the target problem has certain symmetries
$\mathcal{S}$ satisfying $[H,\mathcal{S}]=0$. For symmetry-conserved systems
with $[U,\mathcal{S}]=0$, it is possible to estimate the expectation value of
a unitary $\braket{\psi_{0}}{U}{\psi_{0}}$ without ancilla, as proposed in
o2021error ; xu2023quantum ; cortes2022quantum . However, the issue here is
that as the unitary operator is realised by a Trotter-LCU expansion in Eq. 35,
at least there exists a Pauli operator $W_{i_{q}}$ that does not commute with
$\mathcal{S}$. In this section, we first briefly introduce how to measure
$\braket{\psi_{0}}{U}{\psi_{0}}$ without ancilla, where the unitary is either
$U=e^{-iHt}$ or $U=e^{-iHt_{1}}\hat{O}e^{-iHt_{2}}$. Then we show how to
design the unitary in Eq. 35 such that $[U,\mathcal{S}]=0$.

The expectation value $\braket{\psi_{0}}{U}{\psi_{0}}$ can be expressed as
$\braket{\psi_{0}}{U}{\psi_{0}}=re^{i\theta}$. If $\ket{\psi_{0}}$ is a
product state, or can be prepared as $\ket{\psi_{0}}=U_{p}\ket{0}^{\otimes
n}$, the amplitude of the expectation value
$r=|\braket{\psi_{0}}{U}{\psi_{0}}|$ can be obtained by measuring
$U_{p}^{\dagger}UU_{p}\ket{0}^{\otimes n}$ in the computational basis. The
next step is to obtain the phase $\theta$, for which we can make use of the
fact that the unitary operation conserves the symmetry of $\mathcal{S}$,
$[U,\mathcal{S}]=0$. To do so, we need to introduce a reference state
$\ket{\psi_{\mathrm{Ref}}}$, which lies in a different sector of the initial
state, such that $\braket{\psi_{\mathrm{Ref}}}{U}{\psi_{0}}=0$. We assume that
$\braket{\psi_{\mathrm{Ref}}}{U}{\psi_{\mathrm{Ref}}}$ can be computed
classically. In addition, the following state can be prepared

$\ket{\phi_{0}}=U_{p}\ket{\psi_{0}}=\frac{1}{\sqrt{2}}(\ket{\psi_{\mathrm{Ref}}}+\ket{\psi_{0}}).$

As $r_{s}=|\braket{\psi_{0}}{U_{p}^{\dagger}UU_{p}}{\psi_{0}}|$ can be
measured on a computational basis, the phase $\theta$ can be computed.

The requirements of this protocol include the preparation of the superposition
state and efficient computation of
$\braket{\psi_{\mathrm{Ref}}}{U}{\psi_{\mathrm{Ref}}}$. Let us take the
fermionic Hamiltonian as an example (note that the Heisenberg model can be
regarded as a 1D Fermi-Hubbard model after the Jordan-Wigner transformation).
If the reference state is the vacuum state with the number of particles being
zero, then $\braket{\psi_{\mathrm{Ref}}}{U}{\psi_{\mathrm{Ref}}}$ can be
computed classically. The Hartree-Fock state takes the simple product-state
form $\ket{\psi_{0}}=\ket{1}^{\otimes N_{e}}\ket{0}^{\otimes
N-N_{e}}\in\\{\ket{N_{e}}\\}$, where the first $N_{e}$ qubits are prepared in
the $\ket{1}$ states and the rest of the qubits remain in the $\ket{0}$
states. To prepare the target superposition state, $\ket{\phi_{0}}$, the
Hadamard gate is applied to the first qubit, followed by a ladder of CNOT
gates applied up to the $N_{e}$th qubit, resulting in a total of $N_{e}-1$
CNOT gate operations.

To measure $\braket{\vec{\mathbf{j}}}{O}{\vec{\mathbf{i}}}$, we can use the
circuit in Fig. 1(c) to generate

$U_{p}^{\dagger}(\prod_{q=1}^{\nu}W_{i_{q}}S)^{\dagger}O\prod_{q=1}^{m}W_{i_{q}}SU_{p}\ket{\psi_{0}},$
(42)

and then measure on a computational basis.

In our case, we need to implement the Trotter formula $S$ and the Trotter
remainder $V$. For the implementation of $S$, we need to pair the terms to
ensure each individual term does not break the symmetry. Take the quantum
chemistry problems, $H=\hat{T}+\hat{V}$, for example. We group the term
$\hat{T}_{ij}=h_{ij}(\hat{a}_{i}^{\dagger}\hat{a}_{j}+\hat{a}_{j}^{\dagger}\hat{a}_{i})$
in the kinetic term $\hat{T}$, which conserve the particle number. Its
exponentiation $e^{-i\theta\hat{T}_{ij}}$ can be directly implemented, which
is similar to the implementation for the Heisenberg model case. For the
potential $V$, we can group
$g_{ijkl}(\hat{a}_{i}^{\dagger}\hat{a}_{j}^{\dagger}\hat{a}_{k}\hat{a}_{l}+\hat{a}_{l}^{\dagger}\hat{a}_{k}^{\dagger}\hat{a}_{j}\hat{a}_{i})$
and implement its exponentiation similarly.

The issue is that if we expand the Trotter remainder into Pauli bases, then
applying each individual Pauli term will break the symmetry. One may think of
expanding it into

$V(m)=US^{\dagger}=\sum_{k,l,i,j}\hat{T}_{i}^{k}\hat{V}_{j}^{l}m^{k+l}$ (43)

where $\hat{T}_{i}$ and $\hat{V}_{j}^{l}$ are components of the kinetic and
potential terms, respectively. Although each term conserves the symmetry, the
issue is that $\hat{T}_{i}$ or $\hat{V}_{j}$ is not unitary, and cannot be
implemented directly. To address this issue, we decompose the Hamiltonian into
SWAP and tensor products of Pauli-$Z$ operators rather than Pauli operators.
Let us give an example of the 1D Fermi-Hubbard model, which after the Jordan-
Wigner Transformation takes the form of

$H=J_{1}\sum_{i}X_{i}X_{i+1}+Y_{i}Y_{i+1}+J_{2}\sum_{i}Z_{i}Z_{i+1}+h_{z}\sum_{i}Z_{i}$
(44)

which is the XXZ Heisenberg Hamiltonian. We reformulate it into the following

$H=4J_{1}\sum_{i}\operatorname{SWAP}_{i,i+1}+(J_{2}-J_{1})\sum_{i}Z_{i}Z_{i+1}+h_{z}\sum_{i}Z_{i}.$
(45)

Note that the identity term in the Hamiltonian is a trivial term, so it is
always removed in the Hamiltonian in this work for simplicity. Now, each
individual term is unitary, and commutes with the particle number operator.

For electronic problems in Eq. 6 or more general cases, the Hamiltonian can be
reformulated with SWAP and Pauli $Z$ operators. The kinetic term,
$\hat{T}_{ij}=\frac{h_{ij}}{2}(X_{i}X_{j}\otimes_{k=i+1}^{j-1}Z_{k}+Y_{i}Y_{j}\otimes_{k=i+1}^{j-1}Z_{k})$
can be reformulated as

$\hat{T}_{ij}=h_{ij}(2\operatorname{SWAP}_{i,j}\otimes_{k=i+1}^{j-1}Z_{k}-\frac{1}{2}\otimes_{k=i}^{j}Z_{k}).$

The potential term $\hat{V}_{ijkl}$ can be similarly reformulated such that
each individual term commutes with the particle number operator.

We note that although we reformulate the Hamiltonian in the basis of SWAP and
Pauli Z operators, the implementation of the Trotter formulae remains the same
as that in the Pauli basis. Only the compensation term will be different
compared to the case with one ancillary qubit. As shown in subsubsection
III.4.2, the ancilla-free measurement strategy shows advantages when qubit
connectivity is restricted.

Hamiltonians | NN (ancilla-free) | NN (1-ancilla) | QETU dong2022ground  | arbitrary |
---|---|---|---|---|---
Heisenberg model | $\mathbf{d=\mathcal{O}(1)}$ | $d=\mathcal{O}(n)$ | $d=\mathcal{O}(n)$ | $d=\mathcal{O}(1)$ |
| $g=\mathcal{O}(n)$ | $g=\mathcal{O}(n)$ | $g=\mathcal{O}(n)$ | $g=\mathcal{O}(n)$ |
Electronic structure (Eq. 6) | $\mathbf{d=\mathcal{O}(n)}$ | ${d=\Omega(n^{2})}$ | $d=\Omega(n^{2})$ | $d=\mathcal{O}(n)$ |
| $g=\mathcal{O}(n^{2})$ | $g=\mathcal{O}(n^{3})$ | $g=\mathcal{O}(n^{3})$ | $g=\mathcal{O}(n\log n)$ |
|  |  $\mathbf{\mathbf{d=\mathcal{O}(n)}}$ ($[n/2]$ ancilla) |  |  |
Molecular Hamiltonians (Eq. 50) | $d=\mathcal{O}(n^{2})$ | $d=\Omega(n^{3})$ | $d=\Omega(n^{3})$ | $d=\mathcal{O}(n^{2})$ |
| $g=\mathcal{O}(n^{3})$ | $g=\mathcal{O}(n^{4})$ | $g=\mathcal{O}(n^{4})$ | $g=\mathcal{O}(n^{2}\log n)$ |
Table 4: Circuit depth $d$ and gate count $g$ in each segment when using our
scheme. We consider the restricted qubit structure with nearest-neighbour
architecture or arbitrary architecture (which has all-to-all connectivity). We
compare the ancilla-free or the 1-ancilla schemes. Note that when considering
the all-to-all connectivity, both schemes have the same gate and depth
complexity. We compare our method with the QETU. Note that the control-free
version of QETU has the same depth scaling as the controlled scheme.

###  III.4 Analysis of circuit depth and gate complexity

We first discuss the gate complexity for general Hamiltonians in subsubsection
III.4.1. In subsubsection III.4.2, we analyse the circuit depth and gate
complexity for various Hamiltonians of physical relevance when using either
our method or other early fault-tolerant quantum algorithms. A comparison with
a focus on circuit depth for each elementary block is shown in Table 4 which
supports the results in Table 3.

####  III.4.1 Gate complexity for general cases

The maximum resource appears in the case when $t_{j}=0$, in which case we need
to implement one controlled Trotter formula plus two controlled compensation
terms specified by ${V_{i}}$. The minimum resource appears in the case when
$t_{j}=t_{j}$, in which case we need to implement the Trotter formula plus two
controlled compensation terms specified by ${V_{i}}$. This is because a
controlled Trotter is more costly than a controlled compensation operation, as
analysed in Proposition 3. The resources in other cases with $t_{j}<t_{i}$
will be between Case I and Case II.

We define the following characters of the Hamiltonian,

$\displaystyle\operatorname{wt}(H)$
$\displaystyle:=\sum_{l=1}^{L}\mathrm{wt}\left(P_{l}\right)$ (46)
$\displaystyle\mathrm{wt}_{m}(H)$
$\displaystyle:=\max_{l}\mathrm{wt}\left(P_{l}\right)$

Here, wt$\left(P_{l}\right)$ indicates the weight of the Pauli operator
$P_{l}$, i.e., the number of $\\{X,Y,Z\\}$ terms in $P_{l}$. The following
result gives

######  Proposition 3 (Gate complexity analysis).

Suppose the composite LCU is realised by Algorithm 1 and the Hamiltonian
dynamics $U=e^{-iHt}$ is realised by Eq. 35 with the maximum segment number
$\nu$. The gate complexity of the circuit instance in eigenstate property
estimation is the following: CNOT gate number: $\nu(4\times
5^{k-1}\operatorname{wt}(H)+4\operatorname{wt_{m}}(H)+2\min(n,s_{c}\operatorname{wt_{m}}))-2L$,
single-qubit Pauli rotation gate number: $(2L+4)\nu$.

We give more details in the following. Overall, the circuit within each
segment consists of one controlled Trotter plus controlled compensation terms.
Each of the components within the Trotter formula takes the form of $e^{-iHm}$
with evolution time $m=t/\nu$. The cost by Trotter is given in Lemma 3.

Note that a controlled-Pauli-rotation gate, Ctrl-$\exp(-iP\theta/2)$ with
Pauli operator $P$, can be realised by $2(\operatorname{wt}(H)-1)$ CNOT gates,
a controlled single-qubit $Z$-axis Pauli rotation gate $R_{z}(\theta)$ with
rotation angle $\theta$, and some single-qubit Clifford gates. Furthermore, we
can decompose the controlled single-qubit Pauli rotation gate into two single-
qubit $Z$-axis rotation gates and two CNOT gates. Therefore, the cost for the
controlled Trotter operator Ctrl-$S_{2k}(m)$, in each segment, the number of
CNOT gates are $2\times 5^{k-1}(2\operatorname{wt}(H)-2L)$ and $2L$ controlled
single-qubit Pauli rotations.

The compensation term $V$ consists of a controlled multi-qubit Pauli rotation
operator and some controlled Pauli operators. The controlled multi-qubit Pauli
rotation operator can be realised by $2\operatorname{wt_{m}}(H)$ CNOT gates
and two single-qubit Pauli rotations. The Pauli gates are randomly drawn from
$\\{P_{l}\\}_{l}$ according to the corresponding probability distribution. In
the worst case, the gate sequence length is the truncated value $s_{c}$, and
thus the CNOT gate number is bounded by $s_{c}\operatorname{wt_{m}}(H)$. There
is a saturation when the number of the compensation Pauli operations reaches
$n$. To sum up, for the compensation term, the number of CNOT gates is upper
bounded by

$2\operatorname{wt_{m}}(H)+\min(n,s_{c}\operatorname{wt_{m}})\leq 3n$ (47)

and 2 single-qubit Pauli rotations. We note that since there is a saturation
of the gate count for the Trotter-error-compensation indicated by Eq. 47, we
take $s_{c}$ to be infinity in deriving the asymptotic scaling in gate
complexity in Theorem 1 and Theorem 2.

With the above results, it is easy to check that Case I is more costly than
Case II in Fig. 1. For observable dynamics described in Eq. 36, there are two
compensation terms, but simply one controlled Trotter Ctrl-$S_{2k}$, as
illustrated in Fig. 1(c1). Together with the cost for the Trotter formula, the
number of CNOT gates is upper bounded by

$\nu(4\times
5^{k}\operatorname{wt}(H)-2L+4\operatorname{wt_{m}}(H)+2\min(n,s_{c}\operatorname{wt_{m}})).$
(48)

The number of single-qubit Pauli rotations is $(2L+4)\nu$.

The gate complexity of the ancilla-free measurement strategy can be analysed
in a similar way. We emphasise that the advantage of the ancilla-free
measurement scheme is that it allows the compilation with nearest-neighbour
coupling only, which is discussed in subsubsection III.4.2 in detail.

Figure 3:  An example of how to use $[\frac{n}{2}]$ extra ancillary qubits to
realise controlled exponentiation of potential terms Ctrl-$e^{-ix\hat{V}}$
with a NN architecture. $[n/2]$ ancillary qubits are introduced to reduce the
depth from $\mathcal{O}(n^{3})$ to $\mathcal{O}(n)$. (a) Copy the data on the
qubit $A_{1}$ to $A_{2}$,…$A_{[n/2]}$, and the undo the copy operation. This
circuit is equivalent to using $A_{1}$ as the single-controlled qubit to
control all the other physical qubits $Q_{1}$ to $Q_{n}$. The advantage of
introducing the copy operation is that all the controlled rotation (shaded
blue box in (b)) can be realised by nearest-neighbour gates. (b) The circuit
block controlled $e^{-ix\hat{V}}$ consists of exponentiation of of nearest-
neighbour operations $e^{-iZ_{i}Z_{i+1}}$ followed by $[n/2]$ SWAP operations.
The circuit block in (b) is repeated $\mathcal{O}(n)$ times. An example is
shown in (c). The CNOT operations are omitted in the sub-figure. (c) 1D linear
architecture. $Q_{i}$ represents the $n$th physical qubit (encoding the $n$th
spin orbitals), and $A_{i}$ represents the $i$th ancillary qubit. The
ancillary qubits $A_{i}$ $i=1,2,..,[n/2]$ and actual qubits $Q_{i}$
$i=1,2,..n$ can be placed in the way shown in Step 1 in $\mathcal{O}(n)$
depth. The red arrow connecting $A_{i}$ and $Q_{j}$ is used to represent to
perform the exponentiation of $e^{-iZ_{j}Z_{k}}$ where $Q_{j}$ and $Q_{k}$ are
adjacent qubits. The black arrow connecting different qubits $Q_{j}$ and
$Q_{k}$ is used to represent to perform the corresponding CNOT operation as a
subroutine in the exponentiation of $e^{-iZ_{j}Z_{k}}$, followed by a SWAP
operation. The transformation from the shaded blue and orange boxes in Step 1
to Step 2 can realised by 2 SWAP gates (cyclic swap operation). The rest of
the transformation is realised in the same way. (d) 2D planer architecture.
The qubit connectivity is represented by the orange dashed line.

####  III.4.2 Circuit compilation for typical Hamiltonians and circuit depth
analysis

In this section, we show the circuit depth for typical Hamiltonians including
(1) Heisenberg models, (2) electronic structure in the plane wave dual basis
in Eq. 6, and (3) second-quantised molecular Hamiltonian which takes the
following form

${H}=\sum_{i,j=1}^{n}h_{ij}\hat{a}_{i}^{\dagger}\hat{a}_{j}+\frac{1}{2}\sum_{i,j,k,l,=1}^{n}V_{ijkl}\hat{a}_{i}^{\dagger}\hat{a}_{j}^{\dagger}\hat{a}_{k}\hat{a}_{l},$
(49)

where $n$ is the number of spin orbitals of the molecular system;
$\hat{a}_{i}^{\dagger}$ and $\hat{a}_{i}$ are the fermionic creation and
annihilation operators, respectively; $h_{ij}$ and $V_{ijkl}$ are the
corresponding coefficients for the one-body and two-body interactions,
respectively. For quantum chemistry Hamiltonians, there has been considerable
progress lee2021even ; motta2021low  in efficiently representing the
Hamiltonians with fewer terms and low weights, such as single factorisation
berry2019qubitization  and double factorisation von2021quantum . A common
strategy is to reformulate the two-body fermionic operators as a sum of
squared one-body operators by the Cholesky decomposition motta2021low ,
reformulated as

${H}=\hat{K}+\frac{1}{2}\sum_{\ell}^{\Gamma}\hat{L}_{\ell}^{2},$ (50)

where $\hat{K}$ and $\hat{L}_{\ell}$ are the one-body terms, the number of
terms is $\Gamma=\mathcal{O}(n)$ lee2021even , and the constant has been
removed.

We analyse the circuit depth for the task of eigenstate property estimation
with nearest-neighbour architecture. We will consider two scenarios: (1) the
qubit connectivity is restricted to a linear nearest-neighbour architecture;
(2) arbitrary connectivity is allowed. We will analyse the circuit depth and
gate complexity for typical Hamiltonians using either the ancilla-free scheme
or the one-ancilla scheme.

Within each segment, the circuit consists of two parts: the implementation of
the Trotter formula $S$ and the Trotter remainder $V$. As aforementioned, the
implementation of $S$ is the same for either the ancilla-free scheme or the
one-ancilla scheme. The only difference is the compensation term, which uses
the symmetry-conserved gates as the elementary gates, i.e., SWAP, Pauli-$Z$
gates and their exponentiations. However, we observe that the compensation
term requires at most $d=\mathcal{O}(n)$ and $g=\mathcal{O}(n)$ in the worse
case, which is fewer than the implementation of the Trotter formula in all
cases. In Table 4, we show the gate complexity (including both two-qubit
Clifford gates and non-Clifford gates) within each segment when taking the
qubit connectivity into account. Below, we elaborate on these results.

We first discuss the Heisenberg-type Hamiltonians given by Eq. 44. Each
component $e^{-im(J_{1}X_{i}X_{i+1}+J_{1}Y_{i}Y_{i+1}+J_{2}Z_{i}Z_{i+1})}$
with time duration $m$ can be realised by 3 CNOT gates and 3 single-qubit
Pauli rotation gates as proposed by smith2019simulating . When restricted to a
linear nearest-neighbour architecture, $d_{\textrm{Trotter}}=\mathcal{O}(1)$.
In the ancilla-free scenario, the compensation term is
$\operatorname{SWAP}_{i,i+1}$, Pauli-$Z$ operators and the exponentiations,
resulting in $d=\mathcal{O}(1)$ and $g=\mathcal{O}(n)$. In the one-ancilla
scenario, since we need to implement controlled two-qubit Pauli rotations,
when restricted to linear nearest-neighbour architecture, we need to swap the
ancillary qubit sequentially to qubit $1,...n$ and then perform the controlled
rotation. Then, we need to undo the swap to change back the ordering of the
qubits. The compensation term can be done when the ancillary qubit is adjacent
to the target controlled qubit. Compared to the ancilla-free scheme, we cannot
do it in parallel, resulting in $d=\mathcal{O}(n)$ and $g=\mathcal{O}(n)$.
When removing the restriction of the qubit connectivity, both scenarios have
$d=\mathcal{O}(1)$ and $g=\mathcal{O}(n)$.

For the electronic structure problem in Eq. 6, it can be grouped into
${H}=\hat{T}+\hat{V}$. Below, we discuss the cost for the ancilla-free scheme
and the one-ancilla scheme with a linear nearest-neighbour architecture. The
kinetic term is a non-interacting term which can be diagonalised as
$\hat{T}=\hat{C}(\sum_{i}\alpha_{i}\hat{n}_{i})\hat{C}^{\dagger}$. For the
implementation of $e^{-i\theta\hat{T}}$, a NN architecture requires
$d=\mathcal{O}(n)$ and $g=\mathcal{O}(n^{2})$. For the interacting potential
term $e^{-i\theta\hat{V}}$ which consists of $\mathcal{O}(n^{2})$ terms, a NN
architecture requires $d=\mathcal{O}(n)$ and $g=\mathcal{O}(n^{2})$.
Therefore, $d_{\textrm{Trotter}}=\mathcal{O}(n)$ and
$g_{\textrm{Trotter}}=\mathcal{O}(n^{2})$. For the compensation term, since
there is only one term in the form of
$\operatorname{SWAP}_{i,j}\otimes_{k=i+1}^{j-1}Z_{k}$ or
$\otimes_{k=i+1}^{j-1}Z_{k}$ or $\otimes_{k=i}^{j}Z_{k}$, which can be
implemented at most $d_{\textrm{Remainder}}=\mathcal{O}(n)$ and
$g_{\textrm{Remainder}}=\mathcal{O}(n)$. We note that for the 2D Fermi-Hubbard
model, which is a special case of Eq. 6, the compensation term can be
implemented with $d_{\textrm{Remainder}}=1$ given a planer nearest-neighbour
architecture. To sum up, a NN architecture requires $d=\mathcal{O}(n)$ and
$g=\mathcal{O}(n^{2})$.

Figure 4:  Resource estimation (the number of CNOT gates and T gates) for the
eigenstate property estimation for the $50$-site 1D Heisenberg Hamiltonian
with different target precision. The energy gap is fixed as $\Delta=0.5$. We
compare the gate count involved in different methods. The commutation
relations are taken into account. The second-order paired Taylor series method
is used in our result.

For the one-ancilla scheme, the first observation is that we do not need to
implement controlled $\hat{C}$. The rest of the diagonal part is similar to
the case of the Heisenberg model, which requires $d=\mathcal{O}(n)$ and
$g=\mathcal{O}(n)$. The implementation of $e^{-i\theta\hat{T}}$ requires
$d=\mathcal{O}(n)$ and $g=\mathcal{O}(n^{2})$. However, the implementation of
$e^{-i\theta\hat{V}}$ requires $d=\Omega(n^{2})$ and $g=\Omega(n^{2})$ in a NN
architecture. The challenge of realising the controlled $\hat{V}$ with one
control qubit is that at each time, the control qubit can only control one
rotation $R_{z}(\theta_{ij})$ given a pair of $(i,j)$. Therefore, the
controlled rotation cannot be realised in parallel. Naively, we could first
apply a cyclic swap operation and then realise each individual
$R_{z}(\theta_{ij})$, with the total depth $d=\mathcal{O}(n^{3})$. A lower
bound of the depth complexity of realising the controlled $\hat{V}$ is
$\Omega(n^{2})$. The compensation term requires at most $d=\mathcal{O}(n)$ and
$g=\mathcal{O}(n)$. Therefore, a NN architecture requires
$d=\mathcal{O}(n^{2})$ and $g=\mathcal{O}(n^{2})$.

The depth can be reduced to $\mathcal{O}(n)$ provided $[n/2]$ ancillary
qubits. The result is summarised in Proposition 4.

######  Proposition 4.

When restricted to linear nearest-neighbour qubit connectivity, a controlled
exponentiation Ctrl-$e^{i\theta{H}}$ with electronic problems’ Hamiltonian
${H}$ in Eq. 6 and an angle $\theta$ can be realised in depth $\mathcal{O}(n)$
with $[n/2]$ ancillary qubits. The circuit implementation is illustrated in
Fig. 3. When restricted to only one ancillary qubit, the circuit depth is
$\Omega(n^{2})$.

A graphic proof is shown in Fig. 3. Fig. 3(a,b) shows the overall structure of
the circuit with $\mathcal{O}(n)$ depth. Fig. 3(a) shows how to copy the
information of $A_{1}$ to the rest of ancilla and finally give it back to
$A_{1}$. By using the circuit illustrated in Fig. 3(c) for 1D architecture and
Fig. 3(d) for 2D architecture, the controlled-$H$ can be realised with only NN
operations.

When considering an arbitrary architecture, the gate count can be further
reduced from $\mathcal{O}(n^{2})$ to $\mathcal{O}(n\log n)$ using the results
in low2018hamiltonian . The exponentiation of $e^{-ix\hat{T}}$ can be
implemented by FFFT,
$\textrm{FFFT}e^{-i\sum_{i}t_{i}\hat{n}_{i}}\textrm{FFFT}^{\dagger}$, which
can be implemented with $g=\mathcal{O}(n\log n)$ and $d=\mathcal{O}(\log n)$.
The exponentiation of $e^{-ix\hat{V}}$ can be implemented with
$d=\mathcal{O}(n)$ and $g=\mathcal{O}(n\log n)$. In addition, we discuss the
case of molecular Hamiltonians in Eq. 50 when Cholesky decomposition is
applied. The dominant cost comes from the potential terms. We summarise the
results in Table 4.

Our method shows advantages in circuit depth when considering a linear NN
architecture. We emphasise that the circuit depth and gate complexity of our
ancilla-free scheme match the complexity of implementing Trotter formulae,
even though the eigenstate task is more complicated than the dynamics problem.
Table 3 displays the total gate and depth complexity concerning both CNOT and
T gates, which can be obtained by using the cost within each segment shown in
Table 4.

QETU method requires one ancillary qubit because of the phase iteration. For
certain Hamiltonians, if there exisists a single Pauli operator $K_{j}$ such
that it anticommutes with each component of $H$, then the evolution can be
implemented in a control-free way. For example, for Heisenberg models, we can
divide the Hamiltonian into three terms $H=H_{X}+H_{Y}+H_{Z}$ where $H_{X}$,
$H_{Y}$, $H_{Z}$ contains tensor products of Pauli $X$, $Y$, $Y$,
respectively. For $H_{X}$, we can choose
$K_{1}=\otimes_{i\in\mathrm{odd}}X_{i}\otimes_{i\in\mathrm{even}}Y_{i}$.
However, the controlled-K operations still need depth $d_{K}=\mathcal{O}(n)$,
as opposed to $d=\mathcal{O}(1)$ in our ancilla-free scheme. For electronic
structure problems, to reduce the Trotter error and implement the operations
in parallel, we can group the Hamiltonian into $\hat{T}$ and $\hat{V}$, in
which case it is difficult to find a Pauli operation $K$ such that it commutes
with the grouped term. If the Hamiltonian is not grouped, then the circuit
depth will be increased to $d=\mathcal{O}(n^{3})$ in a naive way of
implementing each individual term.

Figure 5:  Resource estimation (the number of CNOT gates and T gates) for the
eigenstate property estimation for the 2-local Hamiltonian,
$H=\sum_{i,j}X_{i}X_{j}+\sum_{i}Z_{i}$, in which case the commutator relation
between different Hamiltonian summands is ignored. Here $4$th-order random
Trotter formula is used, which is marked by an asterisk alongside Trotter in
the legend. We compare the gate count involved in different methods.

##  IV Resource estimations

In this section, we compare the resource costs associated with different
algorithms, focusing on the number of qubits, CNOT gate count, and non-
Clifford gate count across various models. We will consider the one-
dimensional Heisenberg model with NN coupling, 2-local spin model, and
molecules. The algorithms compared in this section include QPE, QSP, and QETU.
The QPE method relies on Hamiltonian simulation, which can be realised by the
Trottersation or qubitised quantum walk (QW), with the latter having a better
$\varepsilon$ scaling though with a larger overhead. We find that the fourth-
order (random) Trotter result is tested that performs the best among all the
Trotter methods, which is consistent with the results in childs2018toward .
For our RLCU method, we consider the Hamiltonian simulation using the $2k$th-
order paring with $k=0,1$ or $2$. To estimate the gate costs of the RLCU
algorithm, we set the normalisation factor $\mu=2$ to ensure that the sample
complexity of the RCLU algorithm is similar to other quantum algorithms.

Figure 6:  Resource estimation (the number of CNOT gates and T gates) for the
eigenstate property estimation for the 1D lattice Heisenberg Hamiltonian.
Commutation relations are considered. We compare the gate count involved in
different methods with a constant gap in (a,b). For QSP, the CNOT gate count’s
system size dependence is $\mathcal{O}(n^{1.97})$. Our method’s system size
dependence is $\mathcal{O}(n^{1.20})$. QETU’s system size dependence is
$\mathcal{O}(n^{1.32})$. (c) The number of CNOT gates for the eigenstate
property estimation for the 1D lattice Heisenberg Hamiltonian in Eq. 51. The
energy gap is fitted by an exponential decaying function.

The quantum circuit is synthesised to CNOT gates, single-qubit Clifford gates
and non-Clifford gates (including single-qubit Z-axis rotation Rz gates and T
gates). The requirements for the gate number in the representative algorithms
are estimated. Since the central objective is for the application in the early
FTQC or NISQ era, we mainly focus on the circuit depth in a single-run
experiment. Therefore, the amplitude amplification is not considered in the
algorithms, which can deterministically prepare the state closer to the true
ground state yet at the cost of a deeper circuit. Our algorithm requires
single-qubit Pauli rotation gate $R_{z}(\theta)$, which has a circuit
compilation overhead when it is compiled into T gates. The optimal ancilla-
free gate synthesis algorithm in Ref. ross2016optimal , which requires
$3\log_{2}(1/\varepsilon)+\mathcal{O}(\log\log(1/\varepsilon))~{}T$-gates to
approximate the $R_{z}(\theta)$ gate to a precision $\varepsilon$. It is worth
noting that $R_{z}(\theta)$ can be virtually implemented with real physical
devices. For superconducting devices, $R_{z}(\theta)$ indeed does not have to
be implemented physically, but rather it can be implemented by changing the
phase of the reference frame defined by the multi-level rotating frame. That
is, $R_{z}(\theta)$ is a virtual gate, and therefore, there is no physical
error in implementing $R_{z}(\theta)$.

Figure 7:  The number of CNOT and T gates for the eigenstate property
estimation for the P450 molecules with the A-type active space specified in
goings2022reliably . The simulation of fermionic dynamics with Trotterisation
is not optimised in this work. The total weight $\operatorname{wt}(H)=591$.

The elementary operations in the block encoding of $H$ is the Select and
Prepare operations. In the Appendix, we discuss the cost for these two
operations in detail, such that the cost for algorithms based on block
encoding of $H$ can be analysed. When we consider a real physical model, such
as quantum chemistry problems, the coefficient of the Hamiltonian is
constructed by calculating the integral and represents the feature of the
quantum system. Due to finite precision, there will be an amplitude encoding
error when we perform the $\operatorname{PREP}$ operation. To ensure that the
amplitude encoding error is less than a threshold, we require more qubits to
encode the coefficient. However, for the toy models, the absolute value of the
coefficient may not be essentially relevant for the actual physics. For
instance, we can manually set the interaction strength when we study the phase
transition. In this work, we include the amplitude encoding error in our
analysis when aiming for a realistic application. That is why we require more
qubits for the algorithm involving amplitude amplification.

For lattice Hamiltonians, the resource cost can be largely reduced by
considering the commutation relations of the Hamiltonian terms. First, we show
the gate number estimation for Heisenberg models with different target
precisions shown in Fig. 4. Our method shows advantages in both the CNOT gates
and the T gates. It is worth noting that although the logarithmic scaling does
not persist when taking account of the worst-case commutator bound, the actual
performance may be better. The precision dependence is similar to QSP. Fig. 6
clearly shows that our method has a better system size dependence than QSP,
both the gate count and the asymptotic scaling. In the numerics, we set the
energy gap $\Delta=0.01$ and the lower bound of the initial state overlap
$\eta=0.5$.

Next, we show the gate number estimation for spin models with an increasing
system size. We consider the XXZ Hamiltonian

$\displaystyle H=$
$\displaystyle-\sum_{i=1}^{n}(J_{x}X_{i}X_{i+1}+J_{y}Y_{i}Y_{i+1}+J_{z}Z_{i}Z_{i+1})$
(51) $\displaystyle+h_{x}\sum_{i=1}^{n}X_{i}+\sqrt{c^{2}-1}(Z_{1}-Z_{n})$

with the periodic boundary condition and the coupling strengths being
$J_{x}=J_{y}=1$ and $J_{z}=cJ_{x}$. When $h_{x}=0$, it has a gap
$\Delta(c)=4(c-1)$ when the system size is infinity. For $J_{z}=2$ and
$h_{x}=0$, it has an energy gap $\Delta=4$, and the degeneracy is $n+1$. When
$h_{x}$ increases, more excited states will emerge. For $n\leq 100$, the
energy gap is not very small, which can be fitted by an exponential decay
function, $\Delta=b\exp(-{n^{a}})$. We found that $a=1/2$ agrees quite well
with the actual gap. The energy gap could be obtained by fitting, and we have
$\Delta(20)=0.382$, $\Delta(50)=0.124$, $\Delta(100)=0.041$. The results
considering the energy gaps are shown in Fig. 6(c).

For general Hamiltonians without considering commutation relation, our method
may require more T gates than QSP because our method requires ${R_{z}}$ gates,
which have a large overhead when they are synthesised into T gates. To verify
this point, we consider the two-local Hamiltonian. The results are shown in
Fig. 5. It is worth noting that our method only requires rotation gates
$R_{z}(\theta)$ with the rotation angles being identical and small. Therefore,
further reductions may be made by exploiting the protocols in
campbell2016efficient ; duclos2015reducing  and kivlichan2020improved , and we
leave it to future work.

We give the resource cost for the cytochrome enzyme (P450) molecule with the A
active space goings2022reliably . We use the circuits in yordanov2020efficient
for compiling the single- and double-excitation operators, which saves gate
count compared to naive implementing each Trotterised operation of the
fermionic operator. However, we did not optimise the cost for Trotterisation
in every step, in which each term is nonlocal which is why the cost is quite
large. The methods for simulating fermionic Hamiltonians can be directly used
childs2021theory ; kivlichan2020improved ; su2021nearly . We leave the
resource estimations incorporating these fermion simulation methods to future
work.

##  V Outlook and discussion

In this work, we provide a full-stack quantum algorithm based on random
composite LCU for eigenstate property and eigenenergy estimation. Previous
works have established near-optimal query complexity by querying either the
block encoding of the Hamiltonian $H$ (QSP) or the time evolution $e^{-iHt}$
(QETU). Here, we show that the random sampling algorithm can achieve high
precision and good system size scaling when considering gate complexity. We
present an ancilla-free strategy by choosing the elementary operators to be
symmetry-conserved operators (SWAP and Pauli-$Z$ operators). Our method can be
implemented with linear NN architecture, under which constraints, our method
shows advantages over other existing schemes. This is particularly useful for
implementation on NISQ or early FTQC devices since, for most algorithms, due
to the restriction of connectivity, there could be a large overhead when
compiling it into the nearest-neighbour operation. We show a concrete gate
count analysis for various physical models on the CNOT gate and T gate
concerning the circuit synthesis. It shows advantages in both the asymptotic
scaling and the actual gate counts. Our result provides a clear illustration
of the resource requirements for physical systems of interest, such as the
Ising model, the Hubbard model, and molecular systems, with a focus on the
application of quantum computing in the NISQ and early FTQC era.

Similar to other spectral filter methods, our method cannot overcome the
fundamental limitations set up by the initial state and the energy gap. Ref.
lee2023evaluating  numerically investigated its scaling for a few examples,
possibly exponentially vanishing with increasing system size, aligning with
the complexity conjecture. However, on the other side, even if we are given a
sufficiently good initial state with $\eta$ close to $1$, it is BQP-hard to
prepare the ground state of 2-local Hamiltonians (known as the guided local
Hamiltonian problem) gharibian2022dequantizing ; gharibian2022improved .
Therefore, it does not rule out the possibility of quantum speedup over
classical computing. Similar to other random-sampling spectral filter methods
( zhang2022computing ; lin2021heisenberglimited ), the sample complexity with
respect to $\eta$ is $\mathcal{O}(\eta^{-2})$ for which the optimal scaling is
$\mathcal{O}(\eta^{-1/2})$ achieved by QSP or QETU with amplitude
amplification lin2020near , though amplitude amplification also introduces a
large circuit compilation cost. However, the random sampling method does bring
some non-trivial advantages. When considering qubit connectivity, the depth
complexity is nearly independent of the system size, which cannot be achieved
with protocols relying on coherent implementations, particularly QSP or QETU.
Our method is useful in the NISQ or early FTQC eras where the major issue is
the maximum gate (depth) complexity.

Finally, we remark that recent works consider preparing the ground state by
simulation of an open system dynamics described by Lindbladins ding2023single
or a thermal process chen2023local , which may partially overcome the
limitations set up by the initial state lee2023evaluating . Though they
consider simulating the open system dynamics, a fundamental building block of
the method in ding2023single  is querying the real-time evolution. The
resource cost can be analysed using the methods developed in this work. We
leave a detailed comparison between the spectral filter methods and the one
using open system dynamics for future work.

###### Acknowledgements.

J.S. would like to thank Xiao Yuan for early discussion on the project and
Yukun Zhang, Andrew Green for useful related discussions on the manuscript.
J.S. would like to express gratitude to Georg Schusteritsch for providing the
molecular Hamiltonians used in this work. J.S. and M.S.K. thank the UK EPSRC
through EP/Y004752/1, EP/W032643/1, and the KIST Open Research Program. T.G.
thanks UKRI MR/S031545/1, and EPSRC EP/X018180/1 and EP/W032635/1. J.S. would
like to thank Imperial Faculty Researcher Mobility Grant and Schmidt AI in
Science Fellowship.

## References

  * (1) Kishor Bharti, Alba Cervera-Lierta, Thi Ha Kyaw, Tobias Haug, Sumner Alperin-Lea, Abhinav Anand, Matthias Degroote, Hermanni Heimonen, Jakob S. Kottmann, Tim Menke, Wai-Keong Mok, Sukin Sim, Leong-Chuan Kwek, and Alán Aspuru-Guzik.  Noisy intermediate-scale quantum (nisq) algorithms, 2021.
  * (2) Sam McArdle, Suguru Endo, Alan Aspuru-Guzik, Simon C Benjamin, and Xiao Yuan.  Quantum computational chemistry.  Reviews of Modern Physics, 92(1):015003, 2020.
  * (3) Bela Bauer, Sergey Bravyi, Mario Motta, and Garnet Kin-Lic Chan.  Quantum algorithms for quantum chemistry and quantum materials science.  Chemical Reviews, 120(22):12685–12717, 2020.
  * (4) Alexander M Dalzell, Sam McArdle, Mario Berta, Przemyslaw Bienias, Chi-Fang Chen, András Gilyén, Connor T Hann, Michael J Kastoryano, Emil T Khabiboulline, Aleksander Kubica, et al.  Quantum algorithms: A survey of applications and end-to-end complexities.  arXiv preprint arXiv:2310.03011, 2023.
  * (5) Julia Kempe, Alexei Kitaev, and Oded Regev.  The complexity of the local hamiltonian problem.  Siam journal on computing, 35(5):1070–1097, 2006.
  * (6) Seunghoon Lee, Joonho Lee, Huanchen Zhai, Yu Tong, Alexander M Dalzell, Ashutosh Kumar, Phillip Helms, Johnnie Gray, Zhi-Hao Cui, Wenyuan Liu, et al.  Evaluating the evidence for exponential quantum advantage in ground-state quantum chemistry.  Nature Communications, 14(1):1952, 2023.
  * (7) A Yu Kitaev.  Quantum measurements and the abelian stabilizer problem.  arXiv preprint quant-ph/9511026, 1995.
  * (8) Brendon L Higgins, Dominic W Berry, Stephen D Bartlett, Howard M Wiseman, and Geoff J Pryde.  Entanglement-free heisenberg-limited phase estimation.  Nature, 450(7168):393–396, 2007.
  * (9) Emanuel Knill, Gerardo Ortiz, and Rolando D Somma.  Optimal quantum measurements of expectation values of observables.  Physical Review A, 75(1):012328, 2007.
  * (10) Patrick Rall.  Faster coherent quantum algorithms for phase, energy, and amplitude estimation.  Quantum, 5:566, 2021.
  * (11) Richard Meister and Simon C Benjamin.  Resource-frugal hamiltonian eigenstate preparation via repeated quantum phase estimation measurements.  arXiv preprint arXiv:2212.00846, 2022.
  * (12) John M Martyn, Zane M Rossi, Andrew K Tan, and Isaac L Chuang.  Grand unification of quantum algorithms.  PRX Quantum, 2(4):040203, 2021.
  * (13) Zhiyan Ding and Lin Lin.  Even shorter quantum circuit for phase estimation on early fault-tolerant quantum computers with applications to ground-state energy estimation.  PRX Quantum, 4(2):020331, 2023.
  * (14) Yimin Ge, Jordi Tura, and J. Ignacio Cirac.  Faster ground state preparation and high-precision ground energy estimation with fewer qubits.  Journal of Mathematical Physics, 60(2):022202, 2019.
  * (15) Trevor Keen, Eugene Dumitrescu, and Yan Wang.  Quantum algorithms for ground-state preparation and green’s function calculation.  arXiv preprint arXiv:2112.05731, 2021.
  * (16) Dong An, Jin-Peng Liu, and Lin Lin.  Linear combination of hamiltonian simulation for nonunitary dynamics with optimal state preparation cost.  Physical Review Letters, 131(15):150603, 2023.
  * (17) Yongdan Yang, Bing-Nan Lu, and Ying Li.  Accelerated quantum monte carlo with mitigated error on noisy quantum computer.  PRX Quantum, 2(4):040361, 2021.
  * (18) Kianna Wan, Mario Berta, and Earl T. Campbell.  A randomized quantum algorithm for statistical phase estimation, 2021\.
  * (19) Lin Lin and Yu Tong.  Heisenberg-limited ground state energy estimation for early fault-tolerant quantum computers, 2021.
  * (20) Pei Zeng, Jinzhao Sun, and Xiao Yuan.  Universal quantum algorithmic cooling on a quantum computer.  arXiv preprint arXiv:2109.15304, 2021.
  * (21) Sirui Lu, Mari Carmen Bañuls, and J. Ignacio Cirac.  Algorithms for quantum simulation at finite energies.  PRX Quantum, 2:020321, May 2021.
  * (22) Ruizhe Zhang, Guoming Wang, and Peter Johnson.  Computing ground state properties with early fault-tolerant quantum computers.  Quantum, 6:761, 2022.
  * (23) Mingxia Huo and Ying Li.  Shallow trotter circuits fulfil error-resilient quantum simulation of imaginary time, 2021.
  * (24) Guoming Wang, Daniel Stilck França, Ruizhe Zhang, Shuchen Zhu, and Peter D Johnson.  Quantum algorithm for ground state energy estimation using circuit depth with exponentially improved dependence on precision.  Quantum, 7:1167, 2023.
  * (25) Zhiyan Ding, Haoya Li, Lin Lin, HongKang Ni, Lexing Ying, and Ruizhe Zhang.  Quantum multiple eigenvalue gaussian filtered search: an efficient and versatile quantum phase estimation method.  arXiv preprint arXiv:2402.01013, 2024.
  * (26) Min-Quan He, Dan-Bo Zhang, and ZD Wang.  Quantum gaussian filter for exploring ground-state properties.  Physical Review A, 106(3):032420, 2022.
  * (27) Guoming Wang, Daniel Stilck França, Gumaro Rendon, and Peter D Johnson.  Faster ground state energy estimation on early fault-tolerant quantum computers via rejection sampling.  arXiv preprint arXiv:2304.09827, 2023.
  * (28) Samson Wang, Sam McArdle, and Mario Berta.  Qubit-efficient randomized quantum algorithms for linear algebra.  PRX Quantum, 5(2):020324, 2024.
  * (29) Lin Lin and Yu Tong.  Near-optimal ground state preparation.  Quantum, 4:372, 2020.
  * (30) Lin Lin and Yu Tong.  Optimal polynomial based quantum eigenstate filtering with application to solving quantum linear systems.  Quantum, 4:361, 2020.
  * (31) András Gilyén, Yuan Su, Guang Hao Low, and Nathan Wiebe.  Quantum singular value transformation and beyond: exponential improvements for quantum matrix arithmetics.  In Proceedings of the 51st Annual ACM SIGACT Symposium on Theory of Computing, pages 193–204, 2019.
  * (32) Toby S Cubitt.  Dissipative ground state preparation and the dissipative quantum eigensolver.  arXiv preprint arXiv:2303.11962, 2023.
  * (33) Zhiyan Ding, Lin Lin, et al.  Single-ancilla ground state preparation via lindbladians.  arXiv preprint arXiv:2308.15676, 2023.
  * (34) Chi-Fang Chen, Hsin-Yuan Huang, John Preskill, and Leo Zhou.  Local minima in quantum systems.  arXiv preprint arXiv:2309.16596, 2023.
  * (35) Jarrod R McClean, Jonathan Romero, Ryan Babbush, and Alán Aspuru-Guzik.  The theory of variational hybrid quantum-classical algorithms.  New Journal of Physics, 18(2):023023, feb 2016.
  * (36) Marco Cerezo, Andrew Arrasmith, Ryan Babbush, Simon C Benjamin, Suguru Endo, Keisuke Fujii, Jarrod R McClean, Kosuke Mitarai, Xiao Yuan, Lukasz Cincio, et al.  Variational quantum algorithms.  Nature Reviews Physics, 3(9):625–644, 2021.
  * (37) Yihao Liu, Min-Quan He, and ZD Wang.  Variational quantum eigensolvers with quantum gaussian filters for solving ground-state problems in quantum many-body systems.  arXiv preprint arXiv:2401.13459, 2024.
  * (38) William J Huggins, Bryan A O’Gorman, Nicholas C Rubin, David R Reichman, Ryan Babbush, and Joonho Lee.  Unbiasing fermionic quantum monte carlo with a quantum computer.  Nature, 603(7901):416–420, 2022.
  * (39) Jia-Jin Feng and Biao Wu.  Escaping local minima with quantum circuit coherent cooling.  Physical Review A, 109(3):032405, 2024.
  * (40) Lucas K Kovalsky, Fernando A Calderon-Vargas, Matthew D Grace, Alicia B Magann, James B Larsen, Andrew D Baczewski, and Mohan Sarovar.  Self-healing of trotter error in digital adiabatic state preparation.  Physical Review Letters, 131(6):060602, 2023.
  * (41) Mario Motta, Chong Sun, Adrian TK Tan, Matthew J O’Rourke, Erika Ye, Austin J Minnich, Fernando GSL Brandão, and Garnet Kin-Lic Chan.  Determining eigenstates and thermal states on a quantum computer using quantum imaginary time evolution.  Nature Physics, 16(2):205–210, 2020.
  * (42) Xiaosi Xu and Ying Li.  Quantum-assisted monte carlo algorithms for fermions.  Quantum, 7:1072, 2023.
  * (43) Yukun Zhang, Yifei Huang, Jinzhao Sun, Dingshun Lv, and Xiao Yuan.  Quantum computing quantum monte carlo.  arXiv preprint arXiv:2206.10431, 2022.
  * (44) Kasra Hejazi, Mario Motta, and Garnet Kin Chan.  Adiabatic quantum imaginary time evolution.  arXiv preprint arXiv:2308.03292, 2023.
  * (45) Amara Katabarwa, Katerina Gratsea, Athena Caesura, and Peter D Johnson.  Early fault-tolerant quantum computing.  arXiv preprint arXiv:2311.14814, 2023.
  * (46) Dhruv Devulapalli, Eddie Schoute, Aniruddha Bapat, Andrew M. Childs, and Alexey V. Gorshkov.  Quantum routing with teleportation, 2022.
  * (47) Note that Ref. zeng2021universal  extends to eigenstate property estimation.
  * (48) Yulong Dong, Lin Lin, and Yu Tong.  Ground state preparation and energy estimation on early fault-tolerant quantum computers via quantum eigenvalue transformation of unitary matrices.  arXiv preprint arXiv:2204.05955, 2022.
  * (49) Pei Zeng, Jinzhao Sun, Liang Jiang, and Qi Zhao.  Simple and high-precision hamiltonian simulation by compensating trotter error with linear combination of unitary operations.  arXiv preprint arXiv:2212.04566, 2022.
  * (50) Isaac H Kim, Ye-Hua Liu, Sam Pallister, William Pol, Sam Roberts, and Eunseok Lee.  Fault-tolerant resource estimate for quantum chemical simulations: Case study on li-ion battery electrolyte molecules.  Physical Review Research, 4(2):023019, 2022.
  * (51) Markus Reiher, Nathan Wiebe, Krysta M Svore, Dave Wecker, and Matthias Troyer.  Elucidating reaction mechanisms on quantum computers.  Proceedings of the national academy of sciences, 114(29):7555–7560, 2017.
  * (52) Joshua J Goings, Alec White, Joonho Lee, Christofer S Tautermann, Matthias Degroote, Craig Gidney, Toru Shiozaki, Ryan Babbush, and Nicholas C Rubin.  Reliably assessing the electronic structure of cytochrome p450 on today’s classical computers and tomorrow’s quantum computers.  Proceedings of the National Academy of Sciences, 119(38):e2203533119, 2022.
  * (53) Dave Wecker, Bela Bauer, Bryan K Clark, Matthew B Hastings, and Matthias Troyer.  Gate-count estimates for performing quantum chemistry on small quantum computers.  Physical Review A, 90(2):022305, 2014.
  * (54) Earl Campbell.  Random compiler for fast hamiltonian simulation.  Physical review letters, 123(7):070503, 2019.
  * (55) Joonho Lee, Dominic W Berry, Craig Gidney, William J Huggins, Jarrod R McClean, Nathan Wiebe, and Ryan Babbush.  Even more efficient quantum computations of chemistry through tensor hypercontraction.  PRX Quantum, 2(3):030305, 2021.
  * (56) Ryan Babbush, Craig Gidney, Dominic W Berry, Nathan Wiebe, Jarrod McClean, Alexandru Paler, Austin Fowler, and Hartmut Neven.  Encoding electronic spectra in quantum circuits with linear t complexity.  Physical Review X, 8(4):041015, 2018.
  * (57) Dominic W Berry, Craig Gidney, Mario Motta, Jarrod R McClean, and Ryan Babbush.  Qubitization of arbitrary basis quantum chemistry leveraging sparsity and low rank factorization.  Quantum, 3:208, 2019.
  * (58) Ian D Kivlichan, Craig Gidney, Dominic W Berry, Nathan Wiebe, Jarrod McClean, Wei Sun, Zhang Jiang, Nicholas Rubin, Austin Fowler, Alán Aspuru-Guzik, et al.  Improved fault-tolerant quantum simulation of condensed-phase correlated electrons via trotterization.  Quantum, 4:296, 2020.
  * (59) Earl T Campbell.  Early fault-tolerant simulations of the hubbard model.  Quantum Science and Technology, 7(1):015007, 2021.
  * (60) Andrew M Childs, Yuan Su, Minh C Tran, Nathan Wiebe, and Shuchen Zhu.  Theory of trotter error with commutator scaling.  Physical Review X, 11(1):011020, 2021.
  * (61) Andrew M Childs, Dmitri Maslov, Yunseong Nam, Neil J Ross, and Yuan Su.  Toward the first quantum simulation with quantum speedup.  Proceedings of the National Academy of Sciences, 115(38):9456–9461, 2018.
  * (62) Ryan Babbush, Nathan Wiebe, Jarrod McClean, James McClain, Hartmut Neven, and Garnet Kin-Lic Chan.  Low-depth quantum simulation of materials.  Physical Review X, 8(1):011044, 2018.
  * (63) Ian D Kivlichan, Jarrod McClean, Nathan Wiebe, Craig Gidney, Alán Aspuru-Guzik, Garnet Kin-Lic Chan, and Ryan Babbush.  Quantum simulation of electronic structure with linear depth and connectivity.  Physical review letters, 120(11):110501, 2018.
  * (64) Guang Hao Low and Nathan Wiebe.  Hamiltonian simulation in the interaction picture.  arXiv preprint arXiv:1805.00675, 2018.
  * (65) Thomas E O’Brien, Stefano Polla, Nicholas C Rubin, William J Huggins, Sam McArdle, Sergio Boixo, Jarrod R McClean, and Ryan Babbush.  Error mitigation via verified phase estimation.  PRX Quantum, 2(2):020317, 2021.
  * (66) Cristian L Cortes and Stephen K Gray.  Quantum krylov subspace algorithms for ground-and excited-state energy estimation.  Physical Review A, 105(2):022417, 2022.
  * (67) Mario Motta, Erika Ye, Jarrod R McClean, Zhendong Li, Austin J Minnich, Ryan Babbush, and Garnet Kin-Lic Chan.  Low rank representations for quantum simulation of electronic structure.  npj Quantum Information, 7(1):83, 2021.
  * (68) Vera von Burg, Guang Hao Low, Thomas Häner, Damian S Steiger, Markus Reiher, Martin Roetteler, and Matthias Troyer.  Quantum computing enhanced computational catalysis.  Physical Review Research, 3(3):033055, 2021.
  * (69) Adam Smith, MS Kim, Frank Pollmann, and Johannes Knolle.  Simulating quantum many-body dynamics on a current digital quantum computer.  npj Quantum Information, 5(1):106, 2019.
  * (70) Neil J Ross and Peter Selinger.  Optimal ancilla-free clifford+ t approximation of z-rotations.  Quantum Inf. Comput., 16(11&12):901–953, 2016.
  * (71) Earl T Campbell and Joe O’Gorman.  An efficient magic state approach to small angle rotations.  Quantum Science and Technology, 1(1):015007, 2016.
  * (72) Guillaume Duclos-Cianci and David Poulin.  Reducing the quantum-computing overhead with complex gate distillation.  Physical Review A, 91(4):042315, 2015.
  * (73) Yordan S Yordanov, David RM Arvidsson-Shukur, and Crispin HW Barnes.  Efficient quantum circuits for quantum computational chemistry.  Physical Review A, 102(6):062612, 2020.
  * (74) Yuan Su, Hsin-Yuan Huang, and Earl T Campbell.  Nearly tight trotterization of interacting electrons.  Quantum, 5:495, 2021.
  * (75) Sevag Gharibian and François Le Gall.  Dequantizing the quantum singular value transformation: hardness and applications to quantum chemistry and the quantum pcp conjecture.  In Proceedings of the 54th Annual ACM SIGACT Symposium on Theory of Computing, pages 19–32, 2022.
  * (76) Sevag Gharibian, Ryu Hayakawa, François Le Gall, and Tomoyuki Morimae.  Improved hardness results for the guided local hamiltonian problem.  arXiv preprint arXiv:2207.10250, 2022.
  * (77) Yukun Zhang, Yifei Huang, Jinzhao Sun, Dingshun Lv, and Xiao Yuan.  Quantum computing quantum monte carlo.  arXiv preprint arXiv:2206.10431, 2022.
  * (78) Alex Bocharov, Martin Roetteler, and Krysta M Svore.  Efficient synthesis of universal repeat-until-success quantum circuits.  Physical review letters, 114(8):080502, 2015.
  * (79) Guang Hao Low and Isaac L Chuang.  Hamiltonian simulation by qubitization.  Quantum, 3:163, 2019.
  * (80) Guang Hao Low and Isaac L Chuang.  Optimal hamiltonian simulation by quantum signal processing.  Physical review letters, 118(1):010501, 2017.
  * (81) Michele Mosca and Priyanka Mukhopadhyay.  A polynomial time and space heuristic algorithm for t-count.  Quantum Science and Technology, 7(1):015003, 2021.
  * (82) Dmitri Maslov.  Advantages of using relative-phase toffoli gates with an application to multiple control toffoli optimization.  Physical Review A, 93(2):022311, 2016.
  * (83) Andrew M Childs, Aaron Ostrander, and Yuan Su.  Faster quantum simulation by randomization.  Quantum, 3:182, 2019.
  * (84) Michael A Nielsen and Isaac Chuang.  Quantum computation and quantum information, 2002.

###### Contents

  1. I Introduction
  2. II Framework and main results
    1. II.1 Overview of the results
  3. III Methods
    1. III.1 Composite LCU form of the spectral filter
    2. III.2 Observable estimation
    3. III.3 Estimation by quantum circuits: general strategy and ancilla-free strategy
      1. III.3.1 General cases
      2. III.3.2 Ancilla-free composite LCU formula and the corresponding measurement scheme for symmetry conserved systems
    4. III.4 Analysis of circuit depth and gate complexity
      1. III.4.1 Gate complexity for general cases
      2. III.4.2 Circuit compilation for typical Hamiltonians and circuit depth analysis
  4. IV Resource estimations
  5. V Outlook and discussion
  6. A Construction of the spectral filter
    1. A.1 Composite LCU formulae for decomposing a nonunitary operator
    2. A.2 Spectral filter by composite LCU formulae
  7. B Eigenstate property estimation
    1. B.1 Estimator
    2. B.2 Error analysis
    3. B.3 Proof of Theorem 1 (Problem 1)
      1. B.3.1 Gate complexity analysis
    4. B.4 Discussions on molecular Hamiltonians
    5. B.5 Effect of energy error on eigenstate property expectation estimation
    6. B.6 Discretised version
  8. C Eigenenergy estimation (Problem 2)
  9. D Circuit compilation and gate cost for block-encoding-based methods
    1. D.1 Stage Setting
    2. D.2 Gate cost
      1. D.2.1 The index enumeration circuit
    3. D.3 The amplitude encoding, select gates, and reflection gates
    4. D.4 Gate cost for Trotter methods
  10. E Ground state property estimation by quantum signal processing proposed by Lin and Tong
    1. E.1 Overview
  11. F Ground state property estimation with phase estimation
    1. F.1 Complexity of phase estimation
    2. F.2 Hamiltonian simulation by Trotterisation
    3. F.3 Hamiltonian simulation by qubitised quantum walk

In Appendix A, we discuss the construction of the spectral filter with
composite LCU. We discuss the error in observation estimation when given an
LCU form. In Appendix B, we discuss the building blocks of the eigenstate
property estimation algorithm outlined in Algorithm 1. We will follow the
proof idea illustrated in Fig. 2 and Sec. III to prove the main theorems for
eigenstate property estimation in the main text. In Appendix C, we prove the
result for eigenenergy estimation in the main text. Appendix D analyses the
circuit compilation cost for block-encoding-based methods. Appendix E and
Appendix F discuss the cost when using QSP and QPE.

##  Appendix A Construction of the spectral filter