# Calculating the many-body density of states on a digital quantum computer

Alessandro Summer<EMAIL_ADDRESS>School of Physics, Trinity College Dublin,
Dublin 2, Ireland  Cecilia Chiaracane<EMAIL_ADDRESS>School of Physics,
Trinity College Dublin, Dublin 2, Ireland  Mark T. Mitchison
<EMAIL_ADDRESS>School of Physics, Trinity College Dublin, Dublin 2,
Ireland  John Goold<EMAIL_ADDRESS>School of Physics, Trinity College Dublin,
Dublin 2, Ireland

###### Abstract

Quantum statistical mechanics allows us to extract thermodynamic information
from a microscopic description of a many-body system. A key step is the
calculation of the density of states, from which the partition function and
all finite-temperature equilibrium thermodynamic quantities can be calculated.
In this work, we devise and implement a quantum algorithm to perform an
estimation of the density of states on a digital quantum computer which is
inspired by the kernel polynomial method. Classically, the kernel polynomial
method allows to sample spectral functions via a Chebyshev polynomial
expansion. Our algorithm computes moments of the expansion on quantum hardware
using a combination of random state preparation for stochastic trace
evaluation and a controlled unitary operator. We use our algorithm to estimate
the density of states of a non-integrable Hamiltonian on the Quantinuum H1-1
trapped ion chip for a controlled register of 18 qubits. This not only
represents a state-of-the-art calculation of thermal properties of a many-body
system on quantum hardware, but also exploits the controlled unitary evolution
of a many-qubit register on an unprecedented scale.

Tags

Quantum Simulation, Kernel Polynomial Method, DQC1, Density of States, Trapped
Ions

##  I Introduction

The idea of using one quantum system to efficiently simulate another one was
the vision of Feynman over 40 years ago [1]. This paradigm is known as quantum
simulation [2, 3, 4, 5] and is expected to be one of the first real
applications of the current generation of quantum computers [6]. In
particular, recent progress has been made in simulating the dynamics of
strongly correlated many-body systems on current devices [7, 8, 9, 10, 11],
albeit with systems which are still too small to compete with calculations on
classical super-computing architectures. The hope is that the achievable
system sizes will eventually become large enough to surpass what is
classically possible.

In terms of using quantum simulators to extract eigen energies of many-body
systems, early ideas include algorithms based on quantum Fourier transform
[12] such as quantum phase estimation [13, 14, 15] and adiabatic state
preparation [16]. The development of algorithms for the extraction of ground-
state energies is central for the promise of being able to perform quantum
chemistry and materials simulations on quantum computers [17, 18, 19, 20] and
ground state energy calculation is a target of many variational quantum
algorithms [21, 22]. Results for finite temperature and excited states are
more scarce. However recent proposals to measure finite temperature
expectation values on hardware include sampling [23, 24, 25] and imaginary
time evolution [26] and more recently algorithms which may have potential for
computing micro-canonical expectation values were proposed in [27].

The more general idea of using quantum computers to do statistical mechanics
is a topic which is gaining traction [27, 28]. In this work we focus on
developing an algorithm that gives a coarse-grained estimate of the density of
states (DOS) based on the classical kernel polynomial method (KPM) [29]. The
KPM provides a reconstruction of a spectral function by means of a Chebyshev
polynomial expansion, weighted by suitable kernels to damp the Gibbs
oscillations that occur due to finite series truncation. Chebyshev moments are
computed iteratively by applying functions of the Hamiltonian on some initial
state. This step is a challenge to implement on quantum hardware.

Block encoding of a Hamiltonian is deeply connected with the Chebyshev
polynomials [30], implementing the Hamiltonian as a quantum walk as exploited
in the context of the KPM in [31] and more generally to estimate physical
properties in [32, 33]. An alternative is to compute the Chebyshev moments
iteratively in a variational quantum algorithm [34] or otherwise overcoming
the problem of implementing the Chebyshev polynomials using suitably defined
Fourier ones [35, 36].

In this work we devise a hybrid algorithm which uses a combination of pseudo-
random state preparation, Hadamard test and Suzuki-Trotter (ST) decomposition
[37] to evaluate Chebyshev moments. These moments are then used in the
standard KPM expansion. We use an arc-cosine approximation of the Hamiltonian
to implement Chebyshev polynomials from standard ST decomposition and
implement our algorithm on the Quantinuum H1-1 trapped ion quantum simulator
[38]. We were able to approximate the DOS of a non-integrable spin chain for
up to 18 qubits using a single ancillary qubit. Our simulations represent one
of the first explorations of the use of near term quantum computers for
calculations in statistical mechanics.

In Section II we introduce the KPM method and discuss its use in the context
of statistical mechanics. We discuss how to compute the DOS and how a pseudo
random state can be used for stochastic trace estimation. In Section III we
explain the quantum algorithm for extracting Chebyshev polynomials and discuss
the subroutines for random state preparation and implementing the arc-cosine
approximation of the Hamiltonian. In Section IV we then introduce the model we
simulate on hardware and the corresponding gate decomposition used to
implement the controlled unitary. We display our results for estimations of
DOS computed using our hybrid algorithm for systems sizes of $12$ and $18$
qubits.

##  II Classical Kernel Polynomial Method for the density of states

###  II.1 Density of states

In this section we give an overview of the classical KPM and discuss how it is
used to calculate the DOS [29]. For a system of $L$ qubits with Hamiltonian
$\hat{H}$, the DOS is defined as

$g(E)=\frac{1}{2^{L}}\sum_{k=0}^{2^{L}-1}\delta(E-E_{k}),$ (1)

where we denote the energy eigenvalues by $E_{k}$ and the corresponding
eigenvectors by $\ket{k}$, i.e. $\hat{H}\ket{k}=E_{k}\ket{k}$. The DOS gives
access to all thermodynamic properties: in particular, the canonical partition
function can be evaluated as

$\mathcal{Z}(\beta)=\int Ee^{-\beta E}g(E)dE.$ (2)

Differentiation of $\mathcal{Z}$ yields any desired thermodynamic quantity:
for example, the energy

$E(\beta)=-\frac{\partial(\log\mathcal{Z}(\beta))}{\partial\beta},$ (3)

and the entropy $S(\beta)=\beta(E(\beta)-F(\beta))$ where
$F(\beta)=-\beta^{-1}\log(\mathcal{Z(\beta)})$ is the free energy.

In order to extract the DOS from a system of size $L$ one would typically need
exact diagonalization, which requires memory resources scaling as
$\mathcal{O}(2^{3L})$. In contrast, the KPM described in the following section
is able to approximate with memory scaling only as $\mathcal{O}(2^{L})$
(combined with stochastic evaluation of the trace, analysed in Section II.4).

###  II.2 Kernel polynomial method

Consider some function $f(x)$ defined on the interval $x\in[-1,1]$. The KPM
provides an optimal approximation of this function by a finite series of $M$
Chebyshev polynomials. Mathematically, it is defined as

$f_{\text{KPM}}(x)=\frac{1}{\pi\sqrt{1-x^{2}}}\sum_{m=0}^{M}\gamma_{m}^{M}\mu_{m}T_{m}(x),$
(4)

where $\gamma_{m}^{M}$ are the kernel coefficients used to damp Gibbs
oscillations, $T_{m}$ are the Chebyshev polynomials, and $\mu_{m}$ are the
Chebyshev moments. While a more detailed discussion can be found in Appendix
A, the polynomials are generally defined as

$T_{m}(x)=\cos(m\arccos(x))\quad\text{with}\ m\in\mathbb{N}_{0}.$ (5)

For example, $T_{0}(x)=1$, $T_{1}(x)=x$, and all higher Chebyshev polynomials
obey the recursion relation

$T_{m}(x)=2xT_{m-1}(x)-T_{m-2}(x).$ (6)

The KPM expansion Eq. (4) is thus reconstructed by computing the corresponding
Chebyshev moments

$\mu_{m}(x)=\int_{-1}^{1}f(x)T_{m}(x)\mathrm{d}x.$ (7)

The KPM can be adapted to estimate general spectral functions involving a
given quantum mechanical Hamiltonian. In this work, we will be interested
specifically in using it to get an estimate for the DOS. We note that recently
the KPM has been implemented using tensor network techniques [39], and was
used to study thermalisation [40] as well as to extract the DOS of simple
lattice gauge theories [41]. One clear advantage of the KPM is that it does
not suffer from the sign problem that is synonymous with Monte-Carlo
simulations.

###  II.3 Chebyshev moments of the DOS

Let us discuss the key steps in computing the KPM approximation of the DOS,
$g(E)$. Since the domain of Chebyshev polynomials is $[-1,1]$, $\hat{H}$ must
have a spectral norm $||\hat{H}||\leq 1$. If not, $\hat{H}$ can be normalised
as

$\hat{H}\mapsto\frac{\hat{H}-a}{b},$ (8)

with

$a\coloneqq\frac{E_{\text{max}}+E_{\text{min}}}{2}\quad\text{and }\quad
b\coloneqq\frac{E_{\text{max}}-E_{\text{min}}}{2-\varepsilon}$ (9)

where $E_{\text{max}}$ and $E_{\text{min}}$ are the largest and smallest
eigenvalues of $\hat{H}$. Additionally, a small cutoff $\varepsilon$ is
introduced to avoid stability issues that can arise close to the boundaries of
the spectrum. After this rescaling is performed, the expression of the moments
$\mu_{m}$ becomes

$\begin{split}\mu_{m}\coloneqq&\int_{-1}^{1}g(E)T_{m}(E)\mathrm{d}E\\\
=&\int_{-1}^{1}\frac{1}{2^{L}}\sum_{k=0}^{2^{L}-1}\delta(E-E_{k})T_{m}(E)\mathrm{dE}\\\
=&\frac{1}{2^{L}}\sum_{k}T_{m}(E_{k})=\frac{1}{2^{L}}\sum_{k}\bra{k}T_{m}(\hat{H})\ket{k}\\\
=&\frac{1}{2^{L}}\text{Tr}[T_{m}(\hat{H})].\end{split}$ (10)

###  II.4 Stochastic trace evaluation

The first task in the extraction of Chebyshev moments of the KPM is the
efficient estimation of the trace of an operator $\hat{X}$ acting on $L$
qubits, as in Eq. (10). For example computational complexity of determining an
element on the main diagonal, $X_{ii}\coloneqq\bra{i}\hat{X}\ket{i}$, is
$\mathcal{O}(2^{L})$, and therefore determining all elements on the main
diagonal requires $\mathcal{O}(2^{2L})$ operations [42]. A common alternative
method to this computational procedure is stochastic trace estimation. The
main idea is to estimate $\hat{X}$ on a set of randomly chosen states. Let
$\\{\ket{r}\\}$ be a set of $R$ random states on $L$ qubits:

$\ket{r}=\sum_{i=0}^{2^{L}-1}c_{ri}\ket{i}\quad\text{with }c_{ri}\in\
\mathbb{C},$ (11)

so that the stochastic estimate will be

$\Theta=\frac{1}{R}\sum_{r}\bra{r}\hat{X}\ket{r}.$ (12)

If $\forall r$ the statistical average is $0$, i.e.

$\mathopen{\hbox{\set@color${\langle}$}\kern-1.94444pt\leavevmode\hbox{\set@color${\langle}$}}c_{ri}\mathclose{\hbox{\set@color${\rangle}$}\kern-1.94444pt\leavevmode\hbox{\set@color${\rangle}$}}=0,$
(13)

(where the
$\mathopen{\hbox{\set@color${\langle}$}\kern-1.94444pt\leavevmode\hbox{\set@color${\langle}$}}\cdot\mathclose{\hbox{\set@color${\rangle}$}\kern-1.94444pt\leavevmode\hbox{\set@color${\rangle}$}}$
average is over $i$), and if

$\begin{split}\mathopen{\hbox{\set@color${\langle}$}\kern-1.94444pt\leavevmode\hbox{\set@color${\langle}$}}c_{ri}c_{r^{\prime}j}\mathclose{\hbox{\set@color${\rangle}$}\kern-1.94444pt\leavevmode\hbox{\set@color${\rangle}$}}&=0,\\\
\mathopen{\hbox{\set@color${\langle}$}\kern-1.94444pt\leavevmode\hbox{\set@color${\langle}$}}c_{ri}^{*}c_{r^{\prime}j}\mathclose{\hbox{\set@color${\rangle}$}\kern-1.94444pt\leavevmode\hbox{\set@color${\rangle}$}}&=\delta_{rr^{\prime}}\delta_{ij},\end{split}$
(14)

then the variance in the estimate of $\text{Tr}[\hat{X}]$ would be

$\big{(}\delta\Theta\big{)}^{2}=\dfrac{1}{R}\left(\text{Tr}[\hat{X}^{2}]+(\mathopen{\hbox{\set@color${\langle}$}\kern-1.94444pt\leavevmode\hbox{\set@color${\langle}$}}\left\lvert
c_{ri}\right\rvert^{4}\mathclose{\hbox{\set@color${\rangle}$}\kern-1.94444pt\leavevmode\hbox{\set@color${\rangle}$}}-2)\sum_{j=0}^{2^{L}-1}X_{jj}^{2}\right)$
(15)

(for more details see the overview in [29]). Therefore, the relative error
will scale as $\mathcal{O}\big{(}1/\sqrt{R2^{L}}\big{)}$, which means that
fewer random states will be needed to achieve a fixed level of precision as
the size of the system increases. Note that if the coefficients of the random
states are distributed as a Gaussian, then the variance of $\Theta$ only
depends on $\text{Tr}[\hat{X}^{2}]$.

The stochastic trace is used in conjunction with the recursion relation Eq.
(6) to estimate the Chebyshev moments. Starting with a random state
$\ket{r}=\ket{r_{0}}$, one can define

$\ket{r_{1}}=\hat{H}\ket{r_{0}},$ (16)

and recursively generate a series of $M$ vectors

$\ket{r_{m}}=2\hat{H}\ket{r_{m-1}}-\ket{r_{m-2}}=T_{m}(\hat{H})\ket{r}.$ (17)

The moments can then be obtained from the overlap
$\braket{r_{m}}{r}=\braket{r}{T_{m}(\hat{H})}{r}$.

One potential issue with computing the moments in this way is that errors make
the precision of $m$-th moment dependent on errors from the evaluation of the
previous ones. As mentioned earlier, the computation of the $\mu_{m}$ is the
major bottleneck of the classical KPM method, limiting the size of the systems
it can be applied to. In the next section we describe our algorithm to compute
the $m$-th Chebyshev moment on a quantum computer.

##  III Kernel polynomial method hybrid quantum algorithm

Figure 1: The proposed circuit to compute the Chebyshev moments on the quantum
computer is shown. A pseudo-random state is generated by means of a random
circuit on a register of $L$ qubits and an ancillary qubit is prepared in the
$\ket{+}$ state with a Hadamard operation. A unitary operation
$e^{im\hat{H}_{K}}$ is then performed on the register and controlled on the
ancillary qubit (yellow rectangular). This is pursued through a controlled
unitary and a rotation $Z$, with $\beta$ the coefficient of the identity in
$\hat{H}_{K}$. The subsequent gates and measurements on the controlled qubit
yield $\text{Tr}[\cos(m\hat{H}_{K})]$ and $\text{Tr}[\sin(m\hat{H}_{K})]$,
with $\hat{H}_{K}$ connected to the arc-cosine expansion defined by Eq. (38)
which can be then connected to the Chebyshev polynomials of the Hamiltonian of
order $m$.

Fig. 1 shows the circuit we use to compute Chebyshev moments on quantum
hardware. Our proposal is reminiscent of the DQC1 protocol [43], which is a
sub-universal computational paradigm that requires $L+1$ qubits, an ancillary
one (that will act as a control) and $L$ qubits initialised in a maximally
mixed state. The circuit consists of a Hadamard test where a controlled
unitary is applied to $L$ qubits and then the control qubit is measured. We
simulate the maximally mixed state by a randomisation procedure. We now
discuss the two main parts of the circuit in detail. Section III.1 analyses
the random circuit we use to mimic the identity register and Section III.2 the
application of the controlled unitary for the extraction of Chebyshev moments.

###  III.1 Stochastic trace evaluation via state randomisation

Figure 2: Shown in this figure is the 6-layer randomised circuit used for
trace evaluation in the $L=12$ (the blue box on L-1 qubit register in Fig. 1).
($a$) A representation that illustrates the connectivity of the circuit, the
grey dots represent the qubits, and the colored lines are the 2-qubit gates
where they are numbered according to their application. Note that this
particular structure is chosen to exploit the all to all connectivity of the
Quantinuum H1-1 device [38]. ($b$) Gate decomposition of the random circuit.
The colored lines represent $\hat{ZZ}(\pi/2)$ gates, and the dark gray squares
are the single-qubit rotations, randomly selected from $\hat{X}(\pi/2)$,
$\hat{Y}(\pi/4)$ and $\hat{Z}(\pi/2)$. The latter two will be removed by the
compilation since they commute with $\hat{ZZ}(\pi/2)$ and anti-commute with
the other operators. Therefore, their effect is to add two more types of
single-qubit rotations.

The generation of random states on quantum computers has garnered considerable
recent attention recently due to their important role in benchmarking [44].
However, producing uniformly random states that sample the Haar distribution
is inefficient in the sense that the number of gates required scales
exponentially with the register size [45, 46]. T-designs are a type of circuit
that replicate moments of the Haar distribution up to order T or lower [47].
They offer some improvement in terms of gate efficiency, with the number of
gates required scaling polynomially with the number of qubits. Nonetheless,
T-designs are still costly and may be overly random for specific purposes.
Here we take an alternative approach based on pseudorandom state generation
[46]. This involves the generation of states that do not uniformly sample the
Haar distribution 111In some particular cases, it has been proven that the
pseudorandom states after a certain depth become T-designs [71], but still
possess the desired properties such as Eqs. (13) and (14). This method has
been shown to generate states that are sufficiently random in an efficient
way.

Using pseudorandom states to stochastically evaluate the trace has been
proposed in various papers [35, 49, 50, 51, 52], and recently used on quantum
hardware to extract high-temperature transport exponents [11]. In order to
generate pseudo-random states, one can use a circuit composed of alternating
layers of 2-qubit gates and layers with random single qubit rotations, as
suggested by Ref. [46] and adopted in Refs. [44, 49, 11] with small
variations. This random composition has become a widely accepted method for
generating this type of random state. In Ref. [49] the random state is aimed
at generation on a quantum computer where the qubits are connected in a ring
geometry. Layers of 2-qubit gates connecting even-odd and odd-even qubits are
alternated and in between them, there are layers of single-qubit gates
randomly chosen among $\\{\hat{X}(\pi/2),\hat{Y}(\pi/2),\hat{Z}(\pi/4)\\}$ so
that the same rotation is not applied to the same qubit sequentially.

In our case, we focus on a variation of this procedure that takes advantage of
the all-to-all connectivity and of the specific set of elementary gates that
can be implemented directly on the Quantinuum H1-1 trapped-ion-based quantum
computer (see Fig. 2). The gate set of the device Quantinuum H1-1 device
includes

$\begin{split}\hat{ZZ}(\theta)&=e^{-i\theta/2\hat{Z}\otimes\hat{Z}},\\\
\hat{Z}(\theta)&=e^{-i\theta/2\hat{Z}},\\\
\hat{U}_{1q}(\theta,\phi)&=e^{-i\theta/2(\cos(\phi)\hat{X}+\sin(\phi)\hat{Y})}.\\\
\end{split}$ (18)

The device that we will use also supports parallelisation (using Quantum
charge-coupled device (QCCD) architecture with five parallel gate zones [53]).
Inspired by existing techniques to create shallow randomisers, we change the
single qubit rotations to be chosen from
$\\{\hat{X}(\pi/2),\hat{Y}(\pi/4),\hat{Z}(\pi/2)\\}$. Note that the former two
rotations can be applied as a single $\hat{U}_{1q}(\theta,\phi)$ gate, while
the latter can be implemented virtually 222Anyway we can noticed that
$\hat{Z}$ commutes or anti-commutes with all the gates so it could be
implemented virtually independently from the device. The 2-qubit gate will now
be a $\hat{ZZ}(\pi/2)$ with a different connectivity. Each $\hat{ZZ}(\pi/2)$
will connect the qubits $2i$ and $(2i+p)\mod L$ with $i=0,\dots,L/2$ and $p$
an odd number called a jump. The jumps are chosen so that each qubit is
narrowly connected to the other. For instance, as shown in Fig. 2, the jumps
of the first half of the layers can be chosen with
$p_{\ell}=-(-1)^{\ell}(2s\ell+1)$ with $s\in\mathbb{N}_{0}$ and $\ell$ the
index of the layer, and with $p_{\ell}=-p_{-\ell}$ for the second half of the
layers. In particular, Fig. 2 shows the random compiler with $L=12$ and $s=1$.
Notice that $s=0$ reproduces the same pattern of 2-qubit gates as in [49].

Following Ref. [49], we quantify the randomising effect of our circuit by
checking how well the half-system entanglement entropy converges to the Page
value [55], which is the entanglement entropy for a typical Haar-random state.
The von Neumann entanglement entropy of a state $\rho_{r}=\lvert
r\rangle\\!\langle r\lvert$ on a space divided in subspaces $\mathcal{H}_{a}$
and $\mathcal{H}_{b}$ is:

$S(\rho_{r(a)})=\text{Tr}[\rho_{r(a)}\ln\rho_{r(a)}]$ (19)

with $\rho_{r(a)}=\text{Tr}_{(b)}[\lvert r\rangle\\!\langle r\lvert]$. As
shown in [55], this value converges to the Page value
$\log(\dim\mathcal{H}_{a})-\dim\mathcal{H}_{a}/(2\dim\mathcal{H}_{b})$ for
random pure states. For $\dim\mathcal{H}_{a}=\dim\mathcal{H}_{b}$ the it
becomes $\log(2^{L/2})-1/2.$

Figure 3: The convergence of different half-chain entanglement entropies to
the expected value, for $L=10,12,14,16$ and $18$ averaged over $20$ different
random states. The shadows are drawn by the values of the individual random
states. ($a$) The half-chain von Neumann entropy $S$ is shown as function of
random circuit depth and shows the saturation to the Page value [55] (dark
grey lines) for different system size for the circuit used in this work
(Par.). ($b$) Comparison of the relative error of the von Neumann half-chain
entropy $\varepsilon_{s}$ at $L=18$ of the three different random circuits:
the one proposed in Richter et al. [49] (Ric.), the one used in this work
(Par.) (see Fig. 2) and one where the 2-qubit gates are applied sequentially
(Seq.). All three circuits lead to fast convergence to a maximally bi-partite
entangled state. We note that the relative error in Seq. seems to be most
favorable but we have used Par. on the hardware here because less layers are
required to get a good estimate of the trace.

In Fig. 3, we present a classical simulation comparing the convergence rates
to the Page value of three different approaches: Ric., which is the approach
suggested by Richter et al. [49]; Par., the approach used in this work; and
Seq., a variation of Par. where the 2-qubit gates are applied as a chain. In
this latter method, the layers of 2-qubit gates are composed of a sequence of
gates, each of these has support overlapping over half of the previous and
half of the succeeding gate. Therefore, we will compare these three methods in
terms of the number of 2-qubit gates used instead of the number of layers.
Although this last constraint limits the implementation on Quantinuum’s
hardware it is helpful in analysing the three different random compilers. In
fact, while Par. is not the fastest at converging to the Page value, it
performs best at estimating the trace, as shown in Fig. 4. In contrast, the
states produced by Seq. converge to the Page value fastest, but they are not
as effective for trace estimation as those built with Par. This suggests that
the entanglement entropy alone is not sufficient to determine a suitable state
for stochastic trace estimation. In fact, we have observed that, as the
entropy converges to the Page value, the fourth moment of the coefficients
$c_{ri}$ converges to 2 (as required by Eq. (15)). This suggests that the
states become Gaussian-distributed at the same rate as they approach the Page
value.

Figure 4:  Convergence of relative error $\varepsilon_{t}$ of the stochastic
trace to the trace of the Hamiltonian defined by Eq. (IV.1) with $L=18$ for
the three random circuits used to generate Fig. 3. The results are averaged
over $20$ different random states from which the shadows are obtained.

###  III.2 Implementing the Chebyshev polynomials

Simulating the unitary evolution operator or other functions of the
Hamiltonian are central tasks in the field of quantum simulation. In recent
years there has been significant progress, with various approximations, like
qDrift [56, 57, 58, 59], LCU [60], Taylor expansion [61, 62] and qubitization
[63, 64].The latter relies on the sequential application of two oracles:
select ($\hat{S}$) and prepare ($\hat{P}$), defined as follows. Consider a
normalised Hamiltonian that can be decomposed into a sum of unitary operators
as

$\hat{H}=\sum_{p=1}^{P}{\omega_{p}}\hat{H}_{p}\qquad\text{with
}\hat{H}_{p}^{2}=\hat{I}\text{ and }\omega_{p}\in\mathbb{R}.$ (20)

The select operator is defined as

$\hat{S}\coloneqq\sum_{p=1}^{P}\lvert p\rangle\\!\langle
p\lvert_{(a)}\otimes\big{(}\hat{H}_{p}\big{)}_{(s)},$ (21)

where $(a)$ denotes an ancillary Hilbert space comprising $a$ qubits, where
$a\geq\lceil\log_{2}P\rceil$ [64]. The prepare operator acts as

$\hat{P}:\ket{0}_{(a)}\mapsto\sum_{p=1}^{P}\sqrt{\omega_{p}}\ket{p}_{(a)}\eqqcolon\ket{P}_{(a)}.$
(22)

These operators can be combined to obtain a block encoding of the Hamiltonian
as

$\bra{P}_{(a)}\hat{S}\ket{P}_{(a)}=\hat{H}.$ (23)

The idea behind qubitization is to exploit the block encoding of $\hat{H}$ in
a quantum walk to generate any function of $\hat{H}$. The walk operator is
defined as

$\begin{split}\hat{W}\coloneqq&\underbrace{\big{(}2\lvert P\rangle\\!\langle
P\lvert-\hat{I}\big{)}_{(a)}}_{\hat{R}_{(a)}}\hat{S}\\\
=&\sum_{p,q}\big{(}2\sqrt{\omega_{p}\omega_{q}}\lvert q\rangle\\!\langle
p\lvert-\lvert p\rangle\\!\langle
p\lvert\big{)}_{(a)}\otimes\big{(}\hat{H}_{p}\big{)}_{(s)}\end{split}$ (24)

where the operator $\hat{R}_{(a)}$ acts as a reflection about the $\ket{P}$
state. The walk operator can be decomposed as

$\begin{split}\hat{W}=&\bigoplus_{k}\begin{pmatrix}E_{k}&\sqrt{I-\left\lvert
E_{k}\right\rvert^{2}}\\\ -\sqrt{I-\left\lvert E_{k}\right\rvert^{2}}&E_{k}\\\
\end{pmatrix}\\\ =&\bigoplus_{k}e^{i\hat{Y}_{(k)}\arccos{E_{k}}}\end{split}$
(25)

where $E_{k}$ are the eigenvalues (with respective eigenstates $\ket{k}$) of
$\hat{H}$ and each term of the direct sum acts on the subspace
$\mathcal{H}_{k}$ generated by $\ket{\phi_{k}}\coloneqq\ket{P}\ket{k}$ and its
orthogonal state
$\ket{\phi_{k}^{\perp}}\propto(\hat{S}-E_{k}\hat{I})\ket{\phi_{k}}$. Likewise
the $\hat{Y}_{(k)}$ operator will act as a Pauli $\hat{Y}$ operator on this
subspace. Eq. (25) is also useful to understand how $\hat{W}$ is isomorphic to

$e^{i\hat{Y}\otimes\arccos\hat{H}}=\begin{pmatrix}\hat{H}&\sqrt{\hat{I}-\hat{H}^{2}}\\\
-\sqrt{\hat{I}-\hat{H}^{2}}&\hat{H}\end{pmatrix}$ (26)

that is the minimal block encoding of $\hat{H}$. A fundamental feature of
$\hat{W}$, on which the efficiency of qubitization is based, is that repeating
it $m$ times and projecting it on $\ket{P}$ generates the Chebyshev
polynomials:

$\bra{P}_{(a)}\hat{W}^{m}\ket{P}_{(a)}=\hat{T}_{m}(\hat{H}).$ (27)

Let us now consider the smallest decomposition of $\hat{H}$ in unitaries,
where Eq. (25) becomes easier to interpet. Given any Hamiltonian $\hat{H}$,
there is a unitary operator $\hat{U}_{H}$ such that

$\hat{H}=\frac{1}{2}(\hat{U}_{H}+\hat{U}_{H}^{\dagger})$ (28)

where

$\hat{U}_{H}\coloneqq\sum_{k}e^{i\arccos{E_{k}}}\lvert k\rangle\\!\langle
k\lvert=e^{i\arccos{\hat{H}}}$ (29)

It can be observed that $\hat{U}_{H}^{2}\neq\hat{I}$ as would seem to be
required by Eq. (20); however, we notice that this condition can be relaxed by
introducing a new walk operator
$\hat{V}\coloneqq\hat{R}_{(a)}\hat{S}^{\dagger}$ and alternating it with
$\hat{W}$.

As shown in detail in Appendix C, an alternative decomposition to Eq. (28) is

$\hat{H}=\frac{1}{2i}\big{(}e^{i\arcsin{\hat{H}}}-e^{-i\arcsin{\hat{H}}}\big{)}.$
(30)

For this decomposition of $\hat{H}$, the prepare operator simplifies to

$\hat{P}:\ket{0}_{(a)}\mapsto\frac{1}{\sqrt{2}}\ket{0}_{(a)}+\frac{1}{\sqrt{2}}\ket{1}_{(a)}=\ket{+}_{(a)}$
(31)

while the select operator becomes

$\hat{S}=e^{i\hat{Z}_{(a)}\otimes\arccos{\hat{H}}_{(s)}}=\begin{pmatrix}\hat{U}_{H}&0\\\
0&\hat{U}_{H}^{\dagger}\end{pmatrix}.$ (32)

The reflection operator becomes

$\hat{R}_{(a)}=\big{(}2\lvert+\rangle\\!\langle+\lvert-\hat{I}\big{)}_{(a)}=\hat{X}_{(a)}$
(33)

and the walk operators can be constructed as

$\hat{W}=\begin{pmatrix}0&\hat{U}_{H}\\\
\hat{U}_{H}^{\dagger}&0\end{pmatrix},\quad\hat{V}=\begin{pmatrix}0&\hat{U}_{H}^{\dagger}\\\
\hat{U}_{H}&0\end{pmatrix}.$ (34)

The isomorphism between $\hat{W}$ and Eq. (26) is now clear since it results
in mapping the ancilla from $\hat{Z}$ to $\hat{Y}$. The quantum walk now
generates

$\bra{+}_{(a)}\underbrace{\dots\hat{V}\hat{W}\hat{V}\hat{W}}_{m\text{ walk
operators}}\ket{+}_{(a)}=\hat{T}_{m}(\hat{H}).$ (35)

Finally, we notice that

$\hat{X}_{(a)}\hat{S}^{\dagger}\hat{X}_{(a)}=\hat{S}.$ (36)

Therefore, the iteration of the select operator already generates the
desirable walk:

$\bra{+}_{(a)}\hat{S}^{m}\ket{+}_{(a)}=\hat{T}_{m}(\hat{H}).$ (37)

We now want to utilise the decomposition in order to exploit it for the
evaluation of the DOS via our KPM inspired algorithm. The KPM generally works
better at the centre of the Hamiltonian spectrum, where exponentially many
states reside. Close to the spectral edges it can become unstable,
particularly at smaller system size, and hence less reliable [29]. Exploiting
the fact that exponentially many eigenstates are at the centre of the
spectrum, one can expand the $\arccos(\hat{H})$ around $E_{k}=0$:

$\begin{split}\arccos(\hat{H})&=\frac{\pi}{2}-\sum_{k=0}^{\infty}\frac{(2k)!}{2^{2k}(k!)^{2}(2k+1)}\hat{H}^{2k+1}\\\
&=\lim_{K\rightarrow\infty}\Bigg{[}\frac{\pi}{2}-\underbrace{\sum_{k=0}^{K}c_{k}\hat{H}^{2k+1}}_{\hat{H}_{K}}\Bigg{]}.\end{split}$
(38)

It follows that the select operator (now depending $K$) becomes

$\hat{S}_{K}=e^{i\hat{Z}\otimes(\pi/2-\hat{H}_{K})}=e^{i\pi\hat{Z}/2}e^{-i\hat{Z}\otimes\hat{H}_{K}},$
(39)

and

$(\hat{S}_{K})^{m}=e^{im\pi\hat{Z}/2}e^{-im\hat{Z}\otimes\hat{H}_{K}}.$ (40)

The select operator now can be implemented with just two controlled operators
as

$\hat{S}_{K}=\lvert 0\rangle\\!\langle 0\lvert\otimes\hat{U}_{H_{K}}+\lvert
1\rangle\\!\langle 1\lvert\otimes\hat{U}_{H_{K}}^{\dagger}.$ (41)

In what follows we choose to implement the $\hat{U}_{H_{K}}$ using the
standard ST decomposition since this decomposition is already known to require
adequate resources on the hardware platform [65, 66]. A key limitation of the
ST formulae in general is that the only operator it can implement is the
evolution operator. However, the combination of the arc-cosine approximation,
ST decomposition and qubitization allow us to approximate polynomials of
$\hat{H}$ which is sufficient for our goal. We leave the explicit qubitization
of $\hat{U}_{H}$ for future work.

From Eq. (40) we see that the Chebyshev moments can be written as

$\begin{split}&T_{2m}\big{(}\hat{H}\big{)}\simeq(-1)^{m}\cos\big{(}m\hat{H}_{K}\big{)},\\\
&T_{2m+1}\big{(}\hat{H}\big{)}\simeq(-1)^{m}\sin\big{(}m\hat{H}_{K}\big{)}.\end{split}$
(42)

Hence, a Hadamard test on top of the random state will correctly implement
this operation, as depicted in Fig. 1. Note that, as shown in Fig. 1, a final
$\hat{Z}(-\beta)$ rotation on the control qubit is needed to implement the
component of $\hat{H}_{K}$ that is proportional to the identity. For $K=0$
this rotation is simply

$\hat{Z}(-\beta)=e^{i(a/2b)\hat{Z}},$ (43)

i.e. $\beta=-a/b$, where $a$ and $b$ are defined in Eq. (8). Note that since
the random state is used to simulate the maximally mixed state, the circuit of
Fig. 1 is a variant of a DQC1.

##  IV Results on the Quantinuum H1-1 device

###  IV.1 Example Model to be simulated

As a test model to implement our algorithm on quantum hardware we choose the
non-integrable spin-$\frac{1}{2}$ XYZ Heisenberg chain with a staggered
interaction along the $Z$ direction:

$\displaystyle\hat{H}=$
$\displaystyle\sum_{i=0}^{L-1}J_{x}\hat{X}_{(i)}\hat{X}_{(i+1)}+J_{y}\hat{Y}_{(i)}\hat{Y}_{(i+1)}+J_{z}\hat{Z}_{(i)}\hat{Z}_{(i+1)}$
$\displaystyle+\sum_{i=0}^{L-1}(-1)^{i}\Lambda\hat{Z}_{(i)}\hat{Z}_{(i+1)}.$
(44)

By choosing $J_{z}=\Lambda$, we can add a small advantage from the perspective
of gate count, without adding any symmetry to the Hamiltonian. Indeed, the
exponentiation of $\hat{H}$ (Eq. (IV.1)) can be split into two non-commuting
terms: interactions between even-odd spins and between odd-even spins. The
latter, when $\Lambda=J_{z}$, have non-zero couplings only along the X and Y
direction. Then, even-odd terms require a gate composition of Fig. 5.($b$),
while odd-even terms require the even shallower Fig. 5.($a$).

Figure 5: Gate decompositions used to implement the controlled unitary (i.e.
the orange gate of Fig. 1) for the Hamiltonian of Eq. (IV.1). The controlled
ST decomposition of this Hamiltonian requires a controlled
$\exp\big{(}i\alpha\hat{X}\otimes\hat{X}+i\beta\hat{Y}\otimes\hat{Y}\big{)}$
for the odd-even spin pairs (a) and a controlled
$\exp\big{(}i\alpha\hat{X}\otimes\hat{X}+i\beta\hat{Y}\otimes\hat{Y}+i\gamma\hat{Z}\otimes\hat{Z}\big{)}$
for the even-odd spin pairs (b), with $\alpha,\beta,\gamma\in\mathbb{C}$. The
hatched area in (a) represents the part that for the first ST step will
naturally mix with the last layer of the random circuit, without requiring
additional gates; the hatched area of (b) in the last ST step can be
neglected. The (c) and (d) respectively represent the conversion of the
C-$\hat{Z}(2\alpha)$ and C-$\hat{X}(2\alpha)$ in the native gate set currently
supported by the System Model H1-1. Note that this particular gate
decomposition was chosen in order to require as few gates as possible on the
control qubit, to enhance parallelisation and reduce the probability of errors
occurring.

###  IV.2 Hardware results

Figure 6: Results for the Hamiltonian of Eq. (IV.1) at $L=12$ with
$J_{x}=1,J_{y}=1/3,J_{z}=\Lambda=1/2$. In ($a$) the estimates for the first 25
Chebyshev moments. The comparison is between Chebyshev moments obtained in
different ways: analytically (Cheb), applying the arc-cosine approximation at
the first order (Arccos), adding the ST approximation with one single step
(ST), and the results from the circuit of Fig. 1 from Quantinuum, the emulator
(H1-1e) and the quantum computer (H1-1). We used 4 different random circuits
with $j=5$ and 1000 shots. Each cross represents the result from a circuit
with a different random circuit. Figure ($b$) is composed by two parts. On
top: the DOS obtained with the five different estimates for the Chebyshev
moments (using a kernel with $M=25$) and the DOS from the exact
diagonalisation (ED). On bottom we reported the comparison of the DOS obtained
with different approximation of the Chebyshev moments w.r.t. the analytical
ones. Figure 7: Similarly to Fig. 6, the results for the Hamiltonian of Eq.
(IV.1) at $L=18$ with $J_{x}=1,J_{y}=1/3,J_{z}=\Lambda=1/2$. In ($a$) the
estimates for the first 50 Chebyshev moments. Here, for the Quantinuum
executions, we used 10 different random circuits with $j=4$ and 1000 shots.
Each cross represents the result from a circuit with a different random
circuit. Figure ($b$) is composed by two parts. On top: the DOS obtained with
the five different estimates for the Chebyshev moments and the DOS from the
exact diagonalisation (ED). Since we computed $\mu_{2}$ only from the H1-1
System, we also included what the DOS estimate from the $\mu_{2}$ from the
real hardware combined with the estimates from the emulator would look like
(H1-1*). On bottom we reported the comparison of the DOS obtained with
different approximation of the Chebyshev moments w.r.t. the analytical ones.

We are now in a position to test our quantum algorithm on physical hardware
using the model described above. To evaluate the effectiveness of the method
at each stage of approximation, we conducted quantum and classical simulations
of the chain defined by Eq. (IV.1) at $L=12$ and $L=18$. In Fig. 6, we compare
the Chebyshev moments obtained by various methods: by analytical calculation
(using eigenvalues obtained from ED), by approximating the $\arccos$ function,
by including the ST approximation, and by simulating the circuit of Fig. 1
with and without noise. To generate the random state, we implemented a series
of gates on the $L$-qubit register via the Par. method, as described in
Section III.1. The Chebyshev moments were computed using the methods described
in Section III.2, taking the arccosine expansion to order $K=0$, i.e.
$\arccos(\hat{H})\approx\pi/2-\hat{H}$, and using a single ST step.
Remarkably, we find that even with these parameters it is possible to achieve
an appreciable precision. We found that hardware errors overcome any
attainable improvements from using arccosine expansions with $K>0$ or more
than a single ST step. More detail on the error analysis is given in Appendix
D.

We found that for $L=12$, the KPM expansion reaches a sufficient degree of
convergence to ED values at $M\sim 25$. Due to resource limitations on the
Quantinuum device we used only four different random states, running 1000
shots for each one of these. We simulated up to $m=7$ on the quantum hardware,
and approximate $\mu_{m}\approx 0$ for $m>7$, since these moments are too
small to be distinguished from zero with our resources (see Appendix D). Fig.
6.($b$) shows that the moments extracted from the quantum computer are almost
indistinguishable from the exact values within the bulk of the spectrum, and
that the KPM with $M=25$ is able to accurately reconstruct the DOS using these
moments.

We next attempted the same calculation on a register of $L=18$, which uses
$19$ out of $20$ qubits currently available on the H1-1 System [38]. Here the
KPM requires $M\sim 50$ to achieve an accurate approximation, while the
moments that can be simulated accurately range up to $m=11$, beyond which
their values become too small to be distinguished from zero given the number
of shots we have used. Simulating a larger system requires even more
resources, and due to the higher costs associated with these larger circuits,
we ultimately only simulated $\mu_{2}$ on the real hardware for the $L=18$
case. Our results demonstrate that the emulator (H1-1e) consistently produces
accurate outcomes and we compute the rest of the moments with it.
Additionally, we note that the fidelity of the $\hat{ZZ}$ gate appears to
improve with smaller angles. Our circuit heavily employs these gates (see Fig.
5), and any discrepancies in the emulator results can be attributed to an
underestimation of these gates’ fidelities. We report the results from this
last system in Fig. 7, where in the plot of the DOS we add a projection
obtained by combining the results from hardware and emulator together.

##  V Conclusions

In summary, we have performed the first estimation of the density of states of
an non-integrable many-body quantum system on a digital quantum simulator. We
have designed and implemented a quantum algorithm that exploits a combination
of the Hadamard test, Suzuki-Trotter decomposition and random state
preparation to extract Chebyshev moments. Proof-of-principle hardware
simulations were performed on registers of $L=12$ and $L=18$ qubits on the
Quantinuum H1-1 ion trap quantum computer, obtaining a good approximation to
the DOS for a non-integrable Hamiltonian in the bulk of the spectrum
(corresponding to high microcanonical temperatures). We explored in detail the
crucial subroutines of stochastic trace evaluation and controlled evolution
with arccosine approximation. We believe that our quantum hardware results
represent the current state-of-the-art, in terms of both the generation of
pseudo-random states and the implementation of controlled unitary operations
on a many-qubit register. We emphasise that the accuracy of our hardware
results has been limited primarily by financial constraints, and not by
fundamental resource scalings nor even by noise on the H1-1 device.

For the DOS we found it was sufficient to take the arccosine expansion to very
low order ($K=0$), which is ultimately due to the concentration of energy
levels at the centre of the spectrum in large systems. We note that our KPM-
inspired approach can easily be tailored to compute finite-temperature
expectation values in the diagonal and micro-canonical ensembles, in addition
to other spectral functions such as the Lehmann representation of multi-time
correlation functions. In these cases, it may be necessary to consider higher-
order expansions (i.e. $K>0$) to account for features away from the centre of
the spectrum. Our methods could also be combined with other quantum algorithms
tailored to compute ground-state and low-lying excited state properties [17,
18, 19, 20, 21, 22], in order to estimate thermodynamic properties across the
full range of temperature scales.

Our estimation of the DOS on current quantum hardware represents an important
step forward towards quantum statistical mechanics calculations on quantum
computers. As the devices improve, we expect that this algorithm and
subroutines can be used to extract useful approximations to thermodynamic
properties in regimes not accessible to state-of-the-art classical numerical
techniques for strongly correlated systems.

##  VI Acknowledgements

JG would like to thank Sabrina Maniscalco for inviting him to take part in an
unconference in Lapland where some of the first discussions related to this
work took place. We thank Microsoft Ireland, in particular Kieran McCorry, for
providing generous funding to run this project and and for providing access to
the Quantinuum machine through Microsoft Azure Quantum. JG is supported by a
SFI-Royal Society University Research Fellowship and acknowledges funding from
European Research Council Starting Grant ODYSSEY (Grant Agreement No. 758403).
MTM is supported by a Royal Society-Science Foundation Ireland University
Research Fellowship (URF\R1\221571), and acknowledges funding from the
European Commission via the Horizon Europe project ASPECTS (Grant Agreement
No. 101080167).

##  Appendix A Chebyshev polynomials

Given a $f$ function, $\\{g_{M}\\}$, a family of orthogonal functions s.t.

$f(x)\approx g_{M}(x)=\sum_{m=0}^{M}a_{m}\Phi_{m}(x)$ (45)

is said to be a good approximation if it approximates $f$ in at least the
square norm:

$\left\lVert
f(x)-g_{M}(x)\right\rVert_{2}=\left[\int\left(f(x)-g_{M}(x)\right)^{2}dx\right]^{1/2}.$
(46)

A family of functions that are frequently used for this scope is the given by
Fourier series. In this case, the polynomials are:

$g_{M}(x)=\frac{1}{2}a_{0}+\sum_{m=1}^{M}\left(a_{m}\cos(mx)+b_{m}\sin(mx)\right)$
(47)

where the coefficients are

$\begin{split}a_{m}&=\frac{1}{\pi}\int f(x)\cos(mx)dx\\\
b_{m}&=\frac{1}{\pi}\int f(x)\sin(mx)dx\end{split}$ (48)

The Fourier decomposition works well with signals, namely with processes that
happen to be periodic and extended in time (its convergence domain is an
infinite strip, symmetric around the real axes [67]). Similarly, we can
evaluate $f$ (and the coefficients $a_{m}$ and $b_{m}$) just in the
$[-\pi,\pi]$ interval. The Fourier series is then exponentially convergent for
periodic functions with derivatives bounded in $[-\pi,\pi]$, namely its
coefficients $a_{m}$ and $b_{m}$ decrease exponentially in $m$:

$a_{m},b_{m}\sim\mathcal{O}[\exp(-qn^{r})]$ (49)

with $n\gg 1$ and $q$ a constant for some $r>0$. If $f$ is even(odd) all the
sine(cosine) coefficients will cancel out and the expansion becomes the
Fourier sine(cosine) series.

Many physical phenomena are not periodic but bounded (occur in a limited
space). So, to make a function behave periodically one can apply a change of
variable

$x=\cos(\theta)$

with then $\theta\in[-\pi,\pi]$ $\rightarrow$ $x\in[-1,1]$. Noticing that
$f(\cos(\theta))=f(\cos(-\theta))$ the expansion reduces to the Fourier cosine
series:

$g_{M}(\cos(\theta))=\sum_{m=0}^{M}a_{m}\cos(m\cos(\theta))$ (50)

where $\cos(m\cos(\theta))$ are the Chebychev polynomials
$T_{m}(\cos(\theta))$. Reformulating this last equation in $x$ we have

$g_{M}(x)=\sum_{m=0}^{M}a_{m}T_{m}(x)$ (51)

where now

$T_{m}(x)=\cos(m\arccos(x)).$ (52)

Due to this strong connection with the Fourier series, the Chebyshev series
inherits all its properties (among which the exponential convergence) with now
the advantage to be working also with non-periodic functions.

We can notice that $\cos(m\arccos(x))$ is actually a polynomial in $x$ by
looking at the identity

$\cos(2\phi)=\cos^{2}(\phi)-1$ (53)

which leads to the following iterative property of the Chebychev polynomials:

$\begin{split}T_{0}(x)&=1,\\\ T_{1}(x)&=x,\\\
T_{m+1}(x)&=2xT_{m}(x)-T_{m-1}(x).\end{split}$ (54)

##  Appendix B Higher orders of the arc-cosine expansion

Figure 8: Comparison of the first 50 analytical Chebyshev moments with
different approximations: the first order arc-cosine approximation (Arccos1),
the second, which involves the $3^{\text{rd}}$ power of the Hamiltonian
(Arccos3), the normal ST approximation at the first order (ST) and, its
version with the updated parameters (updated ST). After $m=15$ the difference
between these last two approximations becomes negligible.

Due to the concentration of eigenstates around zero, in a lot of physical
Hamiltonians, the first-order expansion is already enough to shape the DOS
with a good approximation. However, in case a better precision is required one
can truncate the $\arccos$ approximation to the second order (i.e. $K=1$). At
this order of approximation, a $k$-local Hamiltonian would become $3k$-local.
However, the terms that are actually $3k$-local contribute to a lesser degree.
Moreover, a large part of $\hat{H}^{3}$ will lie on the same operators of
$\hat{H}$. For instance, in the Hamiltonian of Eq. (IV.1) the only
coefficients in $\hat{H}^{3}$ that increase with the system size are the ones
of Pauli operators already present in $\hat{H}$. Indeed, one can simply
implement the same combination of Pauli operators of $\hat{H}$ by updating the
parameters to include the contributions from $\hat{H}^{3}/6$. These updated
parameters will be, considering the case $J_{z}+\Lambda=J_{x}$, for the
$J_{x}$ of the even sites

$J_{x}\bigg{(}1+\frac{1}{6}\Big{(}(9L/2-2)J_{x}^{2}+(3L-4)J_{y}^{2}+3\beta^{2}+6J_{y}\beta\Big{)}\bigg{)},$
(55)

for the $J_{x}$ of the odd sites

$J_{x}\bigg{(}1+\frac{1}{6}\Big{(}(9L/2-6)J_{x}^{2}+(3L-4)J_{y}^{2}+3\beta^{2}\Big{)}\bigg{)}.$
(56)

Similarly for $J_{y}$ of even sites:

$J_{y}+\dfrac{1}{6}\Big{(}(3L-2)J_{y}^{3}+(9L/2-4)J_{y}J_{x}^{2}+3J_{y}\beta^{2}+6J_{x}^{2}\beta\Big{)}$
(57)

and for the odd ones:

$J_{y}\bigg{(}1+\dfrac{1}{6}\Big{(}(3L-2)J_{y}^{2}+(9L/2-8)J_{x}^{2}+3\beta^{2}\Big{)}\bigg{)}.$
(58)

Then, for $z$-coupling our Hamiltonian only has even terms with a coefficient
updated as:

$J_{x}\bigg{(}1+\dfrac{1}{6}\Big{(}((9L/2-6)J_{x}^{2}+(3L-4)J_{y}^{2}+3\beta^{2}+6J_{y}\beta)\Big{)}\bigg{)}$
(59)

and finally, the coefficient of the identity becomes:

$\beta+\dfrac{1}{6}\Big{(}-\beta^{3}-3LJ_{y}^{2}-9/2L\beta
J_{x}^{2}-3LJ_{x}^{2}J_{y}\Big{)}.$ (60)

The comparison of how the ST decomposition is affected by using these
parameters instead of the initial one is reported in Fig. 8.

##  Appendix C Powers of the Hamiltonian

In the appendix, we see how a similar approach to the arc-cosine approximation
can be adapted to implement the powers of the Hamiltonian. In [68] Seki and
Yunoki devise a method to implement powers of a Hamiltonian as the finite time
derivative of the evolution operator, based on the idea that

$\hat{H}^{m}=i^{m}\frac{\mathrm{d}^{m}\hat{U}(t)}{\mathrm{d}t^{m}}\bigg{\lvert}_{t=0}.$
(61)

At a finite time $t$, the authors suggest taking the central derivative and
for the first power of the Hamiltonian the equation becomes

$\hat{H}\simeq\frac{i}{2t}\left(e^{-it\hat{H}}-e^{it\hat{H}}\right).$ (62)

This is equivalent to approximating the Hamiltonian as

$\hat{H}\simeq\frac{e^{it\hat{H}}-e^{-it\hat{H}}}{2it}=\frac{\sin(t\hat{H})}{t}.$
(63)

A closer look reveals that the first order approximation of the sine
corresponds to the first order approximation of the arc-sine in

$\begin{split}\hat{H}&=\frac{e^{i\arcsin{t\hat{H}}}-e^{-i\arcsin{t\hat{H}}}}{2it}\\\
&\simeq\frac{1}{2it}\left(e^{i(t\hat{H})_{K}}-e^{-i(t\hat{H})_{K}}\right).\end{split}$
(64)

Namely an alternative two-unitary decomposition to the one of Eq. (28) with
the arc-sine instead of the arc-cosine in Eq. (29). This offers a path to
increasing the precision of the method in a more controlled way. Furthermore,
the coefficients of the expansion can be tailored so that the powers of
$\hat{H}$ in the expansion will elide together:

$c_{k}=-\dfrac{(-1)^{k}}{(2k+1)!}-\frac{(-1)^{k^{\prime}}(c_{k^{\prime}})^{2k^{\prime}+1}}{(2k^{\prime}+1)!}\vspace{.1cm}$
(65)

if $\exists\ k^{\prime}\in\mathbb{N}$ s.t. $2k^{\prime}+1=\sqrt{2k+1}$, and

$c_{k}=-\frac{(-1)^{k}}{(2k+1)!}-\frac{(-1)^{n^{\prime}}c_{k^{\prime}}^{2n^{\prime}+1}}{(2n^{\prime}+1)!}-\frac{(-1)^{k^{\prime}}c_{n^{\prime}}^{2k^{\prime}+1}}{(2k^{\prime}+1)!}$
(66)

if $\exists\ k^{\prime},n^{\prime}\in\mathbb{N},k^{\prime}\neq n^{\prime}$
s.t. $2k+1=(2k^{\prime}+1)(2n^{\prime}+1)$, and otherwise

$c_{k}=-\frac{(-1)^{k}}{(2k+1)!}.$ (67)

In this way the error will scale as $\mathcal{O}((t\hat{H})^{2K+1})$.

As highlighted, the probability of successfully selecting the combination of
unitaries will depend on the magnitude of the trace of the product of $t$ and
$\hat{H}$. Therefore, if the Hamiltonian has a large trace, as

$\hat{H}=J\sum_{i}\hat{I}+\hat{X}_{(i)}\hat{X}_{(i+1)}+\hat{Y}_{(i)}\hat{Y}_{(i+1)}+\hat{Z}_{(i)}\hat{Z}_{(i+1)}$
(68)

it will be possible to take small time steps and still have a non-negligible
trace. Otherwise, if the Hamiltonian is normalised it will necessary to take
$t\lesssim 1$.

##  Appendix D Error analysis

The proposed method seeks to achieve efficient execution on NISQ devices by
prioritizing the creation of shallow, noise-resilient circuits over high
precision. This is done while adhering to the constraint of limited time on
the QPUs 333From this last, it is necessary to prioritize the requirement of
finding good (so that a few are enough) and shallow random circuits rather
than relying on MPS, made up of single layers of single rotation gates. In
fact, the cost of the circuit is dictated by the equation [38]
$\text{HQC}=5+\frac{\text{shots}}{5000}(N_{1q}+10N_{2q}+5N_{m})$ (69) where
$N_{1q},N_{2q}$ are the number of single and double qubit gates, respectively,
while $N_{m}$ is the number of qubits measured. We have to take into account
several sources of error in comparison to the analytical Chebyshev moment $m$.
The first is due to the arc-cosine approximation of Section III.2. In this
case, the error depends on the distribution of the eigenvalues of the
Hamiltonian. For odd polynomials, the total expansion will be

$\begin{split}&\
\hat{T}_{2m+1}(\hat{H})=(-1)^{m}\sin((2m+1)\hat{H}_{K\rightarrow\infty})=\\\
&=(-1)^{m}\sum_{n,k}\frac{(-1)^{n}}{(2n+1)!}\left[\frac{(2m+1)(2k)!}{2^{2k}(k!)^{2}(2k+1)}\hat{H}^{2k+1}\right]^{2n+1}\end{split}$
(70)

and similarly for the even terms

$\begin{split}&\
\hat{T}_{2m}(\hat{H})=(-1)^{m}\cos((2m)\hat{H}_{K\rightarrow\infty})=\\\
&=(-1)^{m}\sum_{n,k}\frac{(-1)^{n}}{(2n)!}\left[\frac{(2m)(2k)!}{2^{2k}(k!)^{2}(2k+1)}\hat{H}^{2k+1}\right]^{2n}.\end{split}$
(71)

The first contributes to the error are

$\varepsilon_{o}\simeq\frac{2m+1}{6}\text{Tr}[\hat{H}^{2K+3}]\ \text{ and }\
\varepsilon_{e}\simeq\frac{m^{2}}{18}\text{Tr}[\hat{H}^{4K+6}]$ (72)

for even and odd moments respectively. To roughly characterize it, we need to
introduce speculation about the spectrum of the system. In many physical
systems, the DOS can be modelled as a Gaussian distribution [70], idea behind
which seems also to be inspired the GIT method [32]. If the variance of the
DOS is $\Delta_{E}$ and $\bar{E}$ the average, the Gaussian curve that
describes the spectrum of $\hat{H}$ is

$g(E)\simeq\frac{1}{\sqrt{2\pi}\Delta_{E}}\exp\left(-\frac{(E-\bar{E})^{2}}{2\Delta_{E}^{2}}\right).$
(73)

For $K=0$, the contributes of Eq. (72) become

$\displaystyle\varepsilon_{o}=$ $\displaystyle\frac{2m+1}{6}\int
E^{3}g(E)\mathrm{d}E$ $\displaystyle\simeq$
$\displaystyle\frac{2m+1}{6}\bar{E}(\bar{E}^{2}+3\Delta_{E}^{2}),$ (74)
$\displaystyle\varepsilon_{e}=$ $\displaystyle\frac{m^{2}}{18}\int
E^{6}g(E)\mathrm{d}E\simeq$ $\displaystyle\simeq$
$\displaystyle\frac{m^{2}}{18}(\bar{E}^{6}+15\Delta_{E}^{2}\bar{E}^{4}+45\Delta_{E}^{4}\bar{E}^{2}+15\Delta_{E}^{6}).$
(75)

Note that following the normalisation of Eq. (8), both $\bar{E}$ and
$\Delta_{E}$ will be less than 1. Choosing $\Delta_{E}\geq\bar{E}$ (as it is
in our case) justifies why in our results we did not see better performances
at even rather than odd moments.

The second to the ST decomposition, to decompose the operator in simple gates,
if the ST order is the first and the amount of steps is $1$ the error will go
as $\mathcal{O}(5^{t}m^{1+\frac{1}{t}}/\varepsilon^{\frac{1}{t}})$ where $t$
is the ST order [37].
The third follows from the stochastic trace estimation of Section II.4, where
the error is $\mathcal{O}\big{(}1/\sqrt{R2^{L}}\big{)}$.
The fourth is due to the circuit and the finiteness of the number of shots:

$\sqrt{\frac{1-\mu_{m}}{R\
\text{shots}}}\sim\mathcal{O}\left(\frac{1}{\sqrt{R\ \text{shots}}}\right).$
(76)

## References

  * Feynman [1982] R. P. Feynman, International Journal of Theoretical Physics 21 (1982).
  * Lloyd [1996] S. Lloyd, Science 273, 1073 (1996).
  * Georgescu  _et al._ [2014] I. M. Georgescu, S. Ashhab, and F. Nori, Rev. Mod. Phys. 86, 153 (2014).
  * Tacchino  _et al._ [2020] F. Tacchino, A. Chiesa, S. Carretta, and D. Gerace, Advanced Quantum Technologies 3, 1900052 (2020).
  * Daley  _et al._ [2022] A. J. Daley, I. Bloch, C. Kokail, S. Flannigan, N. Pearson, M. Troyer, and P. Zoller, Nature 607, 667 (2022).
  * Preskill [2018] J. Preskill, Quantum 2, 79 (2018).
  * Zhukov  _et al._ [2018] A. Zhukov, S. Remizov, W. Pogosov, and Y. E. Lozovik, Quantum Information Processing 17, 1 (2018).
  * Cervera-Lierta [2018] A. Cervera-Lierta, Quantum 2, 114 (2018).
  * Francis  _et al._ [2020] A. Francis, J. Freericks, and A. Kemper, Physical Review B 101, 014411 (2020).
  * Smith  _et al._ [2019] A. Smith, M. Kim, F. Pollmann, and J. Knolle, npj Quantum Information 5, 1 (2019).
  * Keenan  _et al._ [2022] N. Keenan, N. Robertson, T. Murphy, S. Zhuk, and J. Goold, arXiv preprint arXiv:2208.12243 (2022).
  * Abrams and Lloyd [1999a] D. S. Abrams and S. Lloyd, Phys. Rev. Lett. 83, 5162 (1999a).
  * Kitaev [1997] A. Y. Kitaev, Russian Mathematical Surveys 52, 1191 (1997).
  * Cleve  _et al._ [1998] R. Cleve, A. Ekert, C. Macchiavello, and M. Mosca, Proceedings of the Royal Society of London. Series A: Mathematical, Physical and Engineering Sciences 454, 339 (1998).
  * Abrams and Lloyd [1999b] D. S. Abrams and S. Lloyd, Phys. Rev. Lett. 83, 5162 (1999b).
  * Aspuru-Guzik  _et al._ [2005] A. Aspuru-Guzik, A. D. Dutoi, P. J. Love, and M. Head-Gordon, Science 309, 1704 (2005).
  * Kassal  _et al._ [2011] I. Kassal, J. D. Whitfield, A. Perdomo-Ortiz, M.-H. Yung, and A. Aspuru-Guzik, Annual Review of Physical Chemistry 62, 185 (2011), pMID: 21166541.
  * Hastings  _et al._ [2015] M. B. Hastings, D. Wecker, B. Bauer, and M. Troyer, Quantum Information & Computation 15, 1 (2015).
  * Cao  _et al._ [2019] Y. Cao, J. Romero, J. P. Olson, M. Degroote, P. D. Johnson, M. Kieferová, I. D. Kivlichan, T. Menke, B. Peropadre, N. P. Sawaya, _et al._ , Chemical reviews 119, 10856 (2019).
  * de Leon  _et al._ [2021] N. P. de Leon, K. M. Itoh, D. Kim, K. K. Mehta, T. E. Northup, H. Paik, B. Palmer, N. Samarth, S. Sangtawesin, and D. Steuerman, Science 372 (2021).
  * Cerezo  _et al._ [2021] M. Cerezo, A. Arrasmith, R. Babbush, S. C. Benjamin, S. Endo, K. Fujii, J. R. McClean, K. Mitarai, X. Yuan, L. Cincio, and P. J. Coles, Nature Reviews Physics 3, 625 (2021).
  * Bharti  _et al._ [2022] K. Bharti, A. Cervera-Lierta, T. H. Kyaw, T. Haug, S. Alperin-Lea, A. Anand, M. Degroote, H. Heimonen, J. S. Kottmann, T. Menke, W.-K. Mok, S. Sim, L.-C. Kwek, and A. Aspuru-Guzik, Rev. Mod. Phys. 94, 015004 (2022).
  * Temme  _et al._ [2011] K. Temme, T. J. Osborne, K. G. Vollbrecht, D. Poulin, and F. Verstraete, Nature 471, 87 (2011).
  * Chowdhury  _et al._ [2017] A. N. Chowdhury, , and R. D. Somma, Quantum Information and Computation 17, 41 (2017).
  * Cohn  _et al._ [2020] J. Cohn, F. Yang, K. Najafi, B. Jones, and J. K. Freericks, Phys. Rev. A 102, 022622 (2020).
  * Motta  _et al._ [2019] M. Motta, C. Sun, A. T. K. Tan, M. J. O’Rourke, E. Ye, A. J. Minnich, F. G. S. L. Brandão, and G. K.-L. Chan, Nature Physics 16, 205 (2019).
  * Lu  _et al._ [2021] S. Lu, M. C. Bañuls, and J. I. Cirac, PRX Quantum 2, 020321 (2021).
  * Schuckert  _et al._ [2022] A. Schuckert, A. Bohrdt, E. Crane, and M. Knap, arXiv preprint arXiv:2206.01756 (2022).
  * Weiße  _et al._ [2006] A. Weiße, G. Wellein, A. Alvermann, and H. Fehske, Rev. Mod. Phys. 78, 275 (2006).
  * Childs  _et al._ [2017] A. M. Childs, R. Kothari, and R. D. Somma, SIAM Journal on Computing 46, 1920 (2017).
  * Rall [2020] P. Rall, Phys. Rev. A 102, 022408 (2020).
  * Roggero [2020] A. Roggero, Phys. Rev. A 102, 022409 (2020).
  * Rajput  _et al._ [2022] A. Rajput, A. Roggero, and N. Wiebe, Quantum 6, 780 (2022).
  * Jensen  _et al._ [2022] P. W. K. Jensen, P. D. Johnson, and A. A. Kunitsa, arXiv preprint arXiv:2206.09881 (2022).
  * Wang  _et al._ [2022] H. Wang, J. Nan, X. Qiu, and X. Li, arXiv preprint arXiv:2202.01170 (2022).
  * Hartse and Roggero [2022] J. Hartse and A. Roggero, arXiv preprint arXiv:2211.00790 (2022).
  * Suzuki [1991] M. Suzuki, Journal of Mathematical Physics 32, 400 (1991).
  * Quantinuum [2023] Quantinuum, “H1-1 system,”  (December 2022 - January 2023).
  * Wolf  _et al._ [2014] F. A. Wolf, I. P. McCulloch, O. Parcollet, and U. Schollwöck, Phys. Rev. B 90, 115124 (2014).
  * Yang  _et al._ [2020] Y. Yang, S. Iblisdir, J. I. Cirac, and M. C. Bañuls, Phys. Rev. Lett. 124, 100602 (2020).
  * Papaefstathiou  _et al._ [2021] I. Papaefstathiou, D. Robaina, J. I. Cirac, and M. C. Bañuls, Phys. Rev. D 104, 014514 (2021).
  * Jin  _et al._ [2021] F. Jin, D. Willsch, M. Willsch, H. Lagemann, K. Michielsen, and H. D. Raedt, Journal of the Physical Society of Japan 90, 012001 (2021).
  * Knill and Laflamme [1998] E. Knill and R. Laflamme, Physical Review Letters 81, 5672 (1998).
  * Boixo  _et al._ [2018] S. Boixo, S. V. Isakov, V. N. Smelyanskiy, R. Babbush, N. Ding, Z. Jiang, M. J. Bremner, J. M. Martinis, and H. Neven, Nature Physics 14, 595 (2018).
  * Nielsen and Chuang [2002] M. A. Nielsen and I. Chuang, American Association of Physics Teachers (2002).
  * Emerson  _et al._ [2003] J. Emerson, Y. S. Weinstein, M. Saraceno, S. Lloyd, and D. G. Cory, Science 302, 2098 (2003).
  * Ambainis and Emerson [2007] A. Ambainis and J. Emerson, Twenty-Second Annual IEEE Conference on Computational Complexity (CCC’07) , 129 (2007).
  * Note [1] In some particular cases, it has been proven that the pseudorandom states after a certain depth become T-designs [71].
  * Richter and Pal [2021] J. Richter and A. Pal, Phys. Rev. Lett. 126, 230501 (2021).
  * Seki and Yunoki [2022] K. Seki and S. Yunoki, Phys. Rev. B 106, 155111 (2022).
  * Goto  _et al._ [2021] S. Goto, R. Kaneko, and I. Danshita, Phys. Rev. B 104, 045133 (2021).
  * Goto  _et al._ [2023] S. Goto, R. Kaneko, and I. Danshita, Phys. Rev. B 107, 024307 (2023).
  * Pino  _et al._ [2021] J. M. Pino, J. M. Dreiling, C. Figgatt, J. P. Gaebler, S. A. Moses, M. S. Allman, C. H. Baldwin, M. Foss-Feig, D. Hayes, K. Mayer, C. Ryan-Anderson, and B. Neyenhuis, Nature 592, 209 (2021).
  * Note [2] Anyway we can noticed that $\hat{Z}$ commutes or anti-commutes with all the gates so it could be implemented virtually independently from the device.
  * Page [1993] D. N. Page, Phys. Rev. Lett. 71, 1291 (1993).
  * Campbell [2019] E. Campbell, Phys. Rev. Lett. 123, 070503 (2019).
  * Chen  _et al._ [2021] C.-F. Chen, H.-Y. Huang, R. Kueng, and J. A. Tropp, PRX Quantum 2, 040305 (2021).
  * Faehrmann  _et al._ [2022] P. K. Faehrmann, M. Steudtner, R. Kueng, M. Kieferova, and J. Eisert, Quantum 6, 806 (2022).
  * Berry  _et al._ [2020] D. W. Berry, A. M. Childs, Y. Su, X. Wang, and N. Wiebe, Quantum 4, 254 (2020).
  * Childs and Wiebe [2012] A. M. Childs and N. Wiebe, Quantum Information and Computation 12, 901 (2012).
  * Berry  _et al._ [2015a] D. W. Berry, A. M. Childs, R. Cleve, R. Kothari, and R. D. Somma, Phys. Rev. Lett. 114, 090502 (2015a).
  * Berry  _et al._ [2015b] D. W. Berry, A. M. Childs, and R. Kothari, in  _2015 IEEE 56th Annual Symposium on Foundations of Computer Science_ (2015) pp. 792–809.
  * Low and Chuang [2017] G. H. Low and I. L. Chuang, Phys. Rev. Lett. 118, 010501 (2017).
  * Low and Chuang [2019] G. H. Low and I. L. Chuang, Quantum 3, 163 (2019).
  * Childs  _et al._ [2018] A. M. Childs, D. Maslov, Y. Nam, N. J. Ross, and Y. Su, Proceedings of the National Academy of Sciences 115, 9456 (2018).
  * Childs  _et al._ [2021] A. M. Childs, Y. Su, M. C. Tran, N. Wiebe, and S. Zhu, Phys. Rev. X 11, 011020 (2021).
  * Boyd [1989] J. P. Boyd,  _Chebyshev & Fourier Spectral Methods_, Lecture Notes in Engineering (Springer, 1989).
  * Seki and Yunoki [2021] K. Seki and S. Yunoki, PRX Quantum 2, 010333 (2021).
  * Note [3] From this last, it is necessary to prioritize the requirement of finding good (so that a few are enough) and shallow random circuits rather than relying on MPS, made up of single layers of single rotation gates. In fact, the cost of the circuit is dictated by the equation [38] $\text{HQC}=5+\frac{\text{shots}}{5000}(N_{1q}+10N_{2q}+5N_{m})$ (77) where $N_{1q},N_{2q}$ are the number of single and double qubit gates, respectively, while $N_{m}$ is the number of qubits measured.
  * Hartmann  _et al._ [2005] M. Hartmann, G. Mahler, and O. Hess, Journal of Statistical Physics 119, 1139 (2005).
  * Harrow and Mehraban [2018] A. Harrow and S. Mehraban, arXiv preprint arXiv:1809.06957 (2018).