# Quantum Ranging Enhanced TDoA Localization

###### Abstract

Localization is critical to numerous applications. The performance of
classical localization protocols is limited by the specific form of distance
information and suffer from considerable ranging errors. This paper foresees a
new opportunity by utilizing the exceptional property of entangled quantum
states to measure a linear combination of target-anchor distances.
Specifically, we consider localization with quantum-based TDoA measurements.
Classical TDoA ranging takes the difference of two separate measurements.
Instead, quantum ranging allows TDoA estimation within a single measurement,
thereby reducing the ranging errors. Numerical simulations demonstrate that
the new quantum-based localization significantly outperforms conventional
algorithms based on classical ranging, with over 50% gains on average.

Index Terms—  quantum ranging, localization, TDoA.

##  1 Introduction

Numerous techniques have been proposed to achieve accurate localization [1],
including range-based methods and range-free methods. Under the classical
range-based localization framework, sensors determine their locations via
nodes with known positions, generally referred to as anchors [2]. Existing
localization protocols rely on ranging techniques including Time of Arrival
(ToA) [3], Time Difference of Arrival (TDoA) [4], and Received Signal Strength
(RSS) [5], among which the information on distances between the target
location and the anchor is encoded into the signal in various ways. Various
modalities have been exploited, such as radio signals (e.g., Wi-Fi, RFID, UWB,
etc), acoustic signals, visible light, magnetic signals, etc.

Due to the limit of signal transmission form under the classical framework,
conventional ranging merely enables measuring single sensor-anchor distances
(ToA, RSS, etc) or the difference of two such distances (TDoA), which
corresponds to a circle and a hyperbolic curve on a plane, respectively. In
this paper, we provide the first glance to introduce the theory of quantum
mechanics [6] and quantum metrology [7, 8, 9] into the field of localization,
which we call quantum-enhanced ranging. Designing a simple entangled probing
quantum state allows a concurrent distance measurement to multiple anchors, in
the form of any $\pm 1$ combination of individual sensor-anchor distances,
within one-shot ranging. This allows a new localization paradigm with reduced
ranging errors.

Localization using quantum ranging, however, raises a non-convex objective
minimization problem. We solve the problem by relaxing it into a convex conic
optimization [10]. Simulations are conducted by simulating erroneous ranging
in the TDoA under the quantum ranging framework, and the comparison to
classical TDoA and CRLB validates the advantage of our proposed localization
technique. The numerical results exhibit that with the aid of quantum ranging,
our method significantly outperforms conventional ranging methods, breaking
the limit of classical systems.

The rest of the paper is organized as follows: §2 introduces quantum ranging,
presents the proposed algorithm, and derives a closed-form expression for
CRLB. §3 presents the evaluation results, and §4 concludes the paper.

##  2 System Model

###  2.1 Quantum Ranging

We propose a new ranging protocol based on the specific property of the
entangled quantum state that incorporates multiple sensor-anchor distances
into the total distance within one-shot ranging. The process of obtaining
distance information is enabled by the time-variation property of an isolated
quantum state. To see this, we first consider the Schrödinger equation [6] on
a quantum state $\ket{\psi(t)}$ that depicts the evolution of entangled
photons state with invariant Hamiltonian $H$ when emitted into the space

$i\hbar\frac{\partial}{\partial t}\ket{\psi(t)}=H\ket{\psi(t)}$ (1)

Suppose the initial state $\ket{\psi(0)}$ is prepared to be an eigenstate of
$H$ corresponds to eigenvalue $E$ under energy expression, Eq.(1) is solved by

$\ket{\psi(t)}=\exp\left(-i\frac{E}{\hbar}t\right)\ket{\psi(0)}$ (2)

When employing a two-qubit addition probe state
$\ket{\psi(0)}=\frac{1}{\sqrt{2}}\left(\ket{0}+\ket{1}\right)$, with state
energy $E_{0},E_{1}$ at states $\ket{0},\ket{1}$ respectively, the evolved
state [11] can be written as

$\ket{\psi(t)}=\frac{1}{\sqrt{2}}e^{-i\frac{E_{0}}{\hbar}t}\left(\ket{0}+e^{-i\theta}\ket{1}\right):~{}\theta=\frac{E_{1}-E_{0}}{\hbar}t$
(3)

Utilizing the methodologies in quantum phase estimation, the information on
the photon’s time of flight (TOF) $t$ can be extracted. Typically, a binary
quantum measurement scenario $\\{\Pi_{0},\Pi_{1}\\}$ is utilized [12, 13],
where $\Pi_{0}=\ket{\psi(0)}\bra{\psi(0)}$ is the projection operator onto the
space of the initial probe state, and $\Pi_{1}=\textbf{id}-\Pi_{0}$ where id
is the identity operator on the Hilbert space $\mathscr{H}$. The probabilities
of obtaining two resulting states are given by

$p_{0}=\left<\psi(t)|\Pi_{0}|\psi(t)\right>=\cos^{2}\frac{\theta}{2},~{}~{}p_{1}=1-p_{0}=\sin^{2}\frac{\theta}{2}$
(4)

Suppose in the $d$-dimensional space, the photon travels from the sensor
location $\bm{x}$ and gets reflected by the anchor $\bm{a}$ where
$\bm{x},\bm{a}\in\mathbb{R}^{d}$. Once get recollected by the sensor, its
state is parametrized by the total phase shift $\theta$, from which the total
distance can be obtained with the additional light speed $c$.

To generalize, we can realize measuring $\pm 1$-combination of the
aforementioned sensor-anchor distances. More specifically, suppose a total of
$n$ anchors $\\{\bm{a}_{1},\dots,\bm{a}_{n}\\}$ are available, for the
$k^{\text{th}}$ ranging, an index set $I_{k}\subseteq\\{1,\dots,n\\}$ picks
several anchors from all available ones. For anchor $i\in I_{k}$, we use
$\omega_{i,k}\in\\{-1,1\\}$ to indicate the sign assignment. To realize such a
ranging scenario, we prepare an entangled probe state correspondingly, which
is explicitly expressed as

$\ket{\psi_{\text{init}}}=\frac{1}{\sqrt{2}}\sum_{j=0,1}\bigotimes_{i\in
I_{k}}\ket{\pi_{i,k}^{j}},~{}~{}\pi^{j}_{i,k}=\omega_{i,k}+j~{}\text{mod}~{}2$
(5)

After the independent time-variation of each photon, we arrive at the
resulting state

$\displaystyle\ket{\psi_{\text{ret}}}$
$\displaystyle=\frac{1}{\sqrt{2}}e^{-i\gamma}\left(\bigotimes_{i\in
I_{k}}\ket{\pi_{i,k}^{0}}+e^{-i\theta}\bigotimes_{i\in
I_{k}}\ket{\pi_{i,k}^{1}}\right)$ (6) $\displaystyle\gamma$
$\displaystyle=\frac{E_{0}}{\hbar}\sum_{\begin{subarray}{c}i\in I_{k}\\\
\omega_{i,k}=1\end{subarray}}\frac{2\|\bm{x}-\bm{a}_{i}\|}{c}+\frac{E_{1}}{\hbar}\sum_{\begin{subarray}{c}i\in
I_{k}\\\ \omega_{i,k}=-1\end{subarray}}\frac{2\|\bm{x}-\bm{a}_{i}\|}{c}$
$\displaystyle\theta$ $\displaystyle=\frac{2(E_{1}-E_{0})}{\hbar
c}\left(\sum_{i\in I_{k}}\omega_{i,k}\|\bm{x}-\bm{a}_{i}\|\right)$

By applying an analogous binary measurement scenario, the estimation on
$\theta$ can be done based on the phase-probability relation as shown in
Eq.(4). The eventual distance combination $\sum_{i\in
I_{k}}\omega_{i,k}\|\bm{x}-\bm{a}_{i}\|$ can thus be obtained by computing
$\hbar c\theta/2(E_{1}-E_{0})$.

###  2.2 Localization through Relaxation

The geometry of our proposed quantum ranging localization can be stated as
follows: On a $\mathbb{R}^{d}$ plane there are a set $A$ of $n$ anchors
$A=\\{\bm{a}_{1},\dots,\bm{a}_{n}\\}$ are available, a sensor with unknown
location $\bm{x}$ attempts to determine its position through $m$ times of
quantum ranging with results $\bm{d}=(d_{1},\dots,d_{m})$. In the presence of
phase estimation error, each distance can be modeled as
$d_{k}=\hat{d_{k}}+\epsilon_{k}$ where $\epsilon_{k}:1\leq k\leq m$ are
independent and identically distributed (i.i.d.) zero-mean Gaussian noise.
Under these assumptions, the localization problem can be formulated as the
following optimization problem:

$\min_{\bm{x}}\sum_{k=1}^{m}\left(\sum_{i\in
I_{k}}\omega_{i,k}\|\bm{x}-\bm{a}_{i}\|-d_{k}\right)^{2}$ (7)

To realize the TDoA ranging, we restrict $|I_{k}|=2$ for all $1\leq k\leq m$,
and the two anchors included by $I_{k}$ are assigned to different signs. Let
$I_{k}=\\{\alpha_{k,1},\alpha_{k,2}\\}$ where $\alpha_{k,v},v=1,2$ indicates
the anchor index included. The Eq.7 is non-convex, but can be transformed into
a conic convex optimization problem. For simplicity, auxiliary variable
$\bm{y}\in\mathbb{R}^{n}$ is introduced to linearize the objective function.
Let incidence matrix $P\in\mathbb{R}^{m\times n}$, the problem can be written
as follows:

$\begin{array}[]{cc}\min\limits_{\bm{x},\bm{y}}&\|P\bm{y}-\bm{d}\|^{2}\\\
\text{s.t.}&y_{i}=\|\bm{x}-\bm{a}_{i}\|\\\
&P=\begin{pmatrix}\bm{e}(\alpha_{1,1},\alpha_{1,2})&\cdots&\bm{e}(\alpha_{m,1},\alpha_{m,2})\end{pmatrix}^{T}\\\
&\forall~{}~{}1\leq i\leq n,\quad 1\leq k\leq m\end{array}$ (8)

where the notation $\bm{e}(i,j)$ stands for the vector in $\mathbb{R}^{n}$
that takes value $1$ at its $i^{\text{th}}$ entry, value $-1$ at its
$j^{\text{th}}$ entry, and remaining entries being $0$. In order to deal with
the non-convex term in the constraint, $\bm{y}$ can be raised to a higher
dimensional space by introducing $\bm{Y}\in\mathbb{S}^{n}$ and
$\gamma\in\mathbb{R}^{+}$. By combining the relaxation techniques proposed
independently by Tseng [14] and Biswas et al. [15] in the sensor network
localization problem, we can reformulate Eq.8 as follows:

$\begin{array}[]{cc}\min\limits_{\bm{x},\bm{y},\bm{Y},\gamma}&P^{T}P\bullet\bm{Y}-2\bm{d}^{T}P\bm{y}\\\
\text{s.t.}&Y_{i,i}=\gamma-2\bm{a}_{i}^{T}\bm{x}+\|\bm{a}_{i}\|^{2},\quad\forall
1\leq i\leq n\\\
&Y_{i,j}\geq\left|\gamma-(\bm{a}_{i}+\bm{a}_{j})^{T}\bm{x}+\bm{a}_{i}^{T}\bm{a}_{j}\right|,~{}\forall
1\leq i<j\leq n\\\ &\gamma\geq\|\bm{x}\|^{2},\quad
y_{i}\geq\|\bm{x}-\bm{a}_{i}\|,\quad\forall 1\leq i\leq n\\\
&\begin{bmatrix}\bm{Y}&\bm{y}\\\ \bm{y}^{T}&1\end{bmatrix}\succeq
0,\quad\bm{Y}\in\mathbb{S}^{n}\end{array}$ (9)

where the operation $A\bullet B$ stands for $\text{Tr}(A^{T}B)$ for symmetric
matrices $A,B$, and $\mathbb{S}^{n}$ is the space of $n\times n$ real
symmetric matrices. Eq.(9) is a fusion of semidefinite programming (SDP) and
second-order cone programming (SOCP), and can be equivalently stated as a pure
SDP problem. It contains decision variables and constraints both up to scale
$O(n^{2})$. An interior-point method [16] can be applied to solve it
efficiently.

###  2.3 Cramér-Rao Lower Bound on Localization Error

It is known that the noise is distance dependent, thus it is convenient to
assume that $\epsilon_{k}\propto{d}_{k}$. We suppose the scalar
$\eta:=\epsilon_{k}/{d}_{k}\in[0,1]$ is set to be identical for all
measurements $k\in\\{1,\dots,m\\}$. By the assumption that
$\epsilon_{k}\sim\mathcal{N}(0,\eta^{2}{d}_{k}^{2})$ as previously suggested,
the probability density function (PDF) of obtaining distance $\bm{d}$
conditioned on real location $\bm{x}$ can be expressed as:

$\begin{array}[]{c}p(\bm{d}|\bm{x})=\\\
\prod\limits_{k=1}^{m}\frac{1}{\sqrt{2\pi}\eta{d}_{k}}\exp\left\\{-\frac{(\|\bm{x}-\bm{a}_{\alpha_{k,1}}\|-\|\bm{x}-\bm{a}_{\alpha_{k,2}}\|-d_{k})^{2}}{2\eta^{2}{d}_{k}^{2}}\right\\}\end{array}$
(10)

An estimator $\hat{\bm{x}}$ on $\bm{x}$ can be obtained by applying the
maximum likelihood estimation (MLE) method, ie,
$\hat{\bm{x}}:=\mathop{\arg\max}_{\bm{x}}p(\bm{d}|\bm{x})$. Maximizing the
likelihood function is equivalent to solving Eq.7 with a weighted version of
the objective function, i.e., the objective function is written as
$P^{T}W^{2}P\bullet\bm{Y}-2\bm{d}^{T}WP\bm{y}$ where
$W=\text{diag}\left\\{1/d_{1},\dots,1/d_{m}\right\\}$.

In the relaxed problem, all parameters are linearized, and thus the estimator
$\hat{\bm{x}}$ induced is an unbiased one. The Cramér-Rao Lower Bound (CRLB)
allows us to derive a lower bound for the localization technique performance.
The Fisher Information Matrix (FIM) $J\in\mathbb{R}^{d\times d}$ is defined as
follows [17]:

$\displaystyle J$ $\displaystyle=-\mathbb{E}\left[\frac{\partial^{2}\log
p(\bm{d}|\bm{x})}{\partial\bm{x}\partial\bm{x}^{T}}\right]$ (11)
$\displaystyle=\sum_{k=1}^{m}\frac{1}{\eta^{2}d_{k}^{2}}\left(\frac{\bm{x}-\bm{a}_{\alpha_{k,1}}}{\|\bm{x}-\bm{a}_{\alpha_{k,1}}\|}-\frac{\bm{x}-\bm{a}_{\alpha_{k,2}}}{\|\bm{x}-\bm{a}_{\alpha_{k,2}}\|}\right)$
$\displaystyle\times\left(\frac{\bm{x}-\bm{a}_{\alpha_{k,1}}}{\|\bm{x}-\bm{a}_{\alpha_{k,1}}\|}-\frac{\bm{x}-\bm{a}_{\alpha_{k,2}}}{\|\bm{x}-\bm{a}_{\alpha_{k,2}}\|}\right)^{T}$

An intuitive metric for localization performance can be given as follows:
Suppose in $r$ times of simulations under the same testbed, the real sensor
position is $\bm{x}^{(l)}$ while the solution returned by arbitrary algorithm
is $\hat{\bm{x}}^{(l)}$. We define the mean error (ME) to specify the solution
deviation:

$\text{ME}:=\frac{1}{r}\sum_{l=1}^{r}\left\|\bm{x}^{(l)}-\hat{\bm{x}}^{(l)}\right\|$
(12)

A non-trivial lower-bound for the expectation of ME in repeated experiments,
by Jenson’s inequality and the CRLB, can be stated as follows:

$\displaystyle\mathbb{E}\left[\text{ME}\right]$
$\displaystyle=\frac{1}{r}\sum_{l=1}^{r}\mathbb{E}\left[\left\|\bm{x}^{(l)}-\hat{\bm{x}}^{(l)}\right\|\right]$
(13)
$\displaystyle\overset{\text{Jenson}}{\geq}\frac{1}{r}\sum_{l=1}^{r}\sqrt{\mathbb{E}\left[(\bm{x}^{(l)}-\hat{\bm{x}}^{(l)})^{T}(\bm{x}^{(l)}-\hat{\bm{x}}^{(l)})\right]}$
$\displaystyle\overset{\text{CRLB}}{\geq}\frac{1}{r}\sum_{l=1}^{r}\sqrt{\text{Tr}\left(J_{l}^{-1}\right)}$

Where $J_{l}$ is the FIM corresponds to the $l^{\text{th}}$ experiment,
determined by the ground truth $\bm{x}^{(l)}$ and the noise factor $\eta$.

##  3 Evaluation

###  3.1 Simulation Settings

In this section, we perform simulation results based on the aforementioned
relaxed problem in the 3-dimensional space. The whole testbed is set on a
classical computer, where the impact of quantum phase estimation error is
simulated by adjusting the noise coefficient $\eta$ throughout the
experiments. Both sensors and anchors are restricted in a cubic region, where
for all sensor position $\bm{x}$ and anchor deployment $\bm{a}$, it holds that
$\bm{x}\in[-D_{x},D_{x}]^{3}$ and $\bm{a}\in[-D_{a},D_{a}]^{3}$ for some fixed
scalars $D_{x},D_{a}\in\mathbb{R}^{+}$. To include both the near-field
($\bm{x}\in\text{conv}(\bm{a}_{1},\dots,\bm{a}_{n})$) and far-field
($\bm{x}\not\in\text{conv}(\bm{a}_{1},\dots,\bm{a}_{n})$) cases, where
$\text{conv}(\cdot)$ stands for the convex hull of a set of points on
$\mathbb{R}^{3}$, we set $D_{x}>D_{a}$. At each noise level, $r=2000$
simulations are conducted.

Moreover, as is mentioned in the previous context, we restrict to explore the
improvement of quantum-assisted localization under an analogous signing to the
TDoA, i.e., distances $\|\bm{x}-\bm{a}_{i}\|-\|\bm{x}-\bm{a}_{j}\|$ for $i\neq
j$ are measured at each time. Details of parameter values are listed in Tab.1.

Table 1: Parameters of testbed geometry parameter | value
---|---
$D_{x}$ |  $2.0$ (m)
$D_{a}$ |  $1.0$ (m)
$n$ | $16$
$m$ | $8$
$\begin{array}[]{c}\text{anchor}\\\ \text{placement}\end{array}$ | $\left\\{\begin{array}[]{c}(1,1,1),(-1,1,1),(1,-1,1)\\\ (-1,-1,1),(1,1,-1),(-1,1,-1)\\\ (1,-1,-1),(-1,-1,-1),(0,0,1)\\\ (0,0,-1),(0,1,0),(0,-1,0)\\\ (1,0,0),(-1,0,0)\\\ (0,0,0.5),(0,0,-0.5)\par\end{array}\right\\}$
$\begin{array}[]{c}\text{ranging}\\\ \text{scenario}\\\ (\text{$y_{i}:=\|\bm{x}-\bm{a}_{i}\|$})\end{array}$ | $\left\\{\begin{array}[]{c}y_{1}-y_{2},y_{3}-y_{4},y_{5}-y_{6}\\\ y_{7}-y_{8},y_{9}-y_{10},y_{11}-y_{12}\\\ y_{13}-y_{14},y_{15}-y_{16}\end{array}\right\\}$
$\eta$ |  $0\%-4\%$, step size $1\%$

In the framework of TDoA with quantum ranging, where the distance error
originates merely from the estimation of phase shift, the error is single-fold
and proportionate to the total distance
$d_{ij}=\|\bm{x}-\bm{a}_{i}\|-\|\bm{x}-\bm{a}_{j}\|$. In the classical TDoA
ranging, the two distances $\|\bm{x}-\bm{a}_{i}\|$ and $\|\bm{x}-\bm{a}_{j}\|$
are separately distributed, and the resulting error of $d_{ij}$ is the
superposition of term-wise errors.

By the above assumption, we model the noise integrated into distance ranging
in following two ways, which subsequently correspond to two TDoA setups:

  * •

Quantum-assisted TDoA. Under noise level $\eta$, the noisy distance can be
modelled as
$d_{k}=(\|\bm{x}-\bm{a}_{\alpha_{k,1}}\|-\|\bm{x}-\bm{a}_{\alpha_{k,2}}\|)(1+\eta\varepsilon_{k})$,
where $\varepsilon_{k}\sim\mathcal{N}(0,1)$.

  * •

Classical TDoA. Each sensor-anchor distance is integrated with noise, i.e.,
the noisy distance is expressed as
$d_{k}=\|\bm{x}-\bm{a}_{\alpha_{k,1}}\|(1+\eta\varepsilon_{k,1})-\|\bm{x}-\bm{a}_{\alpha_{k,2}}\|(1+\eta\varepsilon_{k,2})$,
where $\varepsilon_{k,v}\sim\mathcal{N}(0,1)$ for all $k$ and $v\in\\{1,2\\}$.
It will serve as a baseline for the performance of our proposed quantum-
ranging localization.

It is worth noticing that since the difference between two distances is
measured, for any $\bm{y}^{*}$ falling into the feasible region,
$\bm{y}:=\bm{y}^{*}+\tau\bm{1}$ will yield another feasible $\bm{y}$ with
respect to the objective function, where $\bm{1}$ is the vector with all
entries being $1$ in the $\mathbb{R}^{n}$ space. To avoid large deviations of
solution returned by the solver from the optimal one, we add a penalty term to
the primitive objective function in Eq.(9). With penalty coefficient
$\delta\in\mathbb{R}^{+}$, the problem now is

$\begin{array}[]{cc}\min\limits_{\bm{x},\bm{y},\bm{Y},\gamma}&P^{T}P\bullet\bm{Y}-2\bm{d}^{T}P\bm{y}+\delta\cdot\text{Tr}(\bm{Y})\\\
\text{s.t.}&Y_{i,i}=\gamma-2\bm{a}_{i}^{T}\bm{x}+\|\bm{a}_{i}\|^{2},\quad\forall
1\leq i\leq n\\\
&Y_{i,j}\geq\left|\gamma-(\bm{a}_{i}+\bm{a}_{j})^{T}\bm{x}+\bm{a}_{i}^{T}\bm{a}_{j}\right|,\forall
1\leq i<j\leq n\\\ &\gamma\geq\|\bm{x}\|^{2},\quad
y_{i}\geq\|\bm{x}-\bm{a}_{i}\|,\quad\forall 1\leq i\leq n\\\
&\begin{bmatrix}\bm{Y}&\bm{y}\\\ \bm{y}^{T}&1\end{bmatrix}\succeq
0,\quad\bm{Y}\in\mathbb{S}^{n}\end{array}$ (14)

The CRLB as a function of the ground truth of sensor location and noise factor
for each single randomly generated test case is also computed to serve as a
benchmark, against which the efficiency of our proposed algorithm is tested.

###  3.2 Numerical Results

Our code written in Python generates the simulation settings and solves the
mixed SOCP/SDP problem and baselines. It calls MOSEK through cvxpy (Version
1.3) when solving the conic programming problems (proposed algorithm and pure
SDP), and runs on a 3.2 GHz AMD Ryzen 7 5800H processor and 16 GB RAM. The
aforementioned penalty coefficient $\delta$ is set to be $6\times 10^{-7}$
throughout the experiment.

Fig.1 shows the performance measured by ME of the result solved by our
proposed algorithm Eq.(9) with penalty under both quantum and classical
settings. It can be observed that apart from the noiseless case, i.e.,
$\eta=0$, TDoA with quantum ranging consistently outperforms that with
traditional ranging. With the increase in noise level, the performance of
quantum-assisted TDoA is resilient and the localization error yields a slow
trend of increase, which is analogous to the trend of CRLB.

Fig. 1: ME under various settings with respect to noise

The underlying optimization problem yields decision variables and constraints
both up to scale $O(n^{2})$. To the best of our knowledge [18], a solver would
return the optimal solution within $O(n^{5.246})$ iterations. Tab.2 shows the
respective runtime of the localization algorithms under different ranging
settings. It can be observed that our proposed technique trades off
computational efficiency for localization accuracy. Rather than the algorithm
utilizing SOCP/SDP we suggested, simpler algorithms of higher efficiency can
be designed for specific ranging scenarios, which provides a direction for
future research.

Algorithm | Average time cost (sec)
---|---
quantum-assisted ranging | $0.4348$
classical ranging | $0.4334$
Table 2: Time cost of algorithms throughout the simulation

##  4 Conclusion

In this paper, we describe the first TDoA localization technique under the
framework of quantum-enhanced distance ranging. Conic relaxation methods are
deployed to solve the raised non-convex optimization problem. The analysis of
CRLB is conducted to provide a theoretical lower bound for the localization
model. The numerical simulations demonstrate that the proposed algorithm
significantly outperforms and breaks the limit of TDoA localization algorithms
with classical ranging. Future work shall include designing optimal
measurement scenarios when illumination in the environment is taken into
account, finding generic ranging scenarios with specific anchor placement
geometry, and developing more efficient localization algorithms.

##  5 Acknowledgments

This work is supported in part by National Natural Science Foundation of China
(NSFC) under grant No. 62222216, Guangdong Basic and Applied Basic Research
Foundation (Project No. 2022A1515010340) and by the Hong Kong Research Grant
Council (RGC) through the Early Career Scheme (ECS) grant 27310822, the
General Research Fund (GRF) grants No. 17303923.

## References

  * [1] Chenshu Wu, Zheng Yang, and Yunhao Liu,  Wireless Indoor Localization: A Crowdsourcing Approach,  Springer, 2018.
  * [2] Guoqiang Mao, Barış Fidan, and Brian D.O. Anderson,  “Wireless sensor network localization techniques,”  Computer Networks, vol. 51, no. 10, pp. 2529–2553, 2007.
  * [3] Saleh Al-Jazzar, Mounir Ghogho, and Desmond McLernon,  “A joint toa/aoa constrained minimization method for locating wireless devices in non-line-of-sight environment,”  IEEE Transactions on Vehicular Technology, vol. 58, no. 1, pp. 468–472, 2009.
  * [4] Regina Kaune, Julian Hörst, and Wolfgang Koch,  “Accuracy analysis for tdoa localization in sensor networks,”  in 14th International Conference on Information Fusion, 2011, pp. 1–8.
  * [5] R.W. Ouyang, A.K.-S. Wong, and Chin-Tau Lea,  “Received Signal Strength-Based Wireless Localization via Semidefinite Programming: Noncooperative and Cooperative Schemes,”  IEEE Transactions on Vehicular Technology, vol. 59, no. 3, pp. 1307–1318, Mar. 2010.
  * [6] David J. Griffiths and Darrell F. Schroeter,  Introduction to Quantum Mechanics,  Cambridge University Press, 3 edition, Aug. 2018.
  * [7] Vittorio Giovannetti, Seth Lloyd, and Lorenzo Maccone,  “Quantum Metrology,”  Physical Review Letters, vol. 96, no. 1, pp. 010401, Jan. 2006.
  * [8] Vittorio Giovannetti, Seth Lloyd, and Lorenzo Maccone,  “Advances in quantum metrology,”  Nature Photonics, vol. 5, no. 4, pp. 222–229, Mar. 2011.
  * [9] Vittorio Giovannetti, Seth Lloyd, and Lorenzo Maccone,  “Quantum-enhanced measurements: Beating the standard quantum limit,”  Science, vol. 306, no. 5700, pp. 1330–1336, 2004.
  * [10] Stephen Boyd and Lieven Vandenberghe,  Convex Optimization,  Cambridge University Press, 1 edition, Mar. 2004.
  * [11] J. Anandan and Y. Aharonov,  “Geometry of quantum evolution,”  Physical Review Letters, vol. 65, no. 14, pp. 1697–1700, Oct. 1990.
  * [12] Luca Pezze’ and Augusto Smerzi,  “Quantum theory of phase estimation,”  2014,  Publisher: arXiv Version Number: 1.
  * [13] John Watrous,  The Theory of Quantum Information,  Cambridge University Press, 1 edition, Apr. 2018.
  * [14] Paul Tseng,  “Second‐Order Cone Programming Relaxation of Sensor Network Localization,”  SIAM Journal on Optimization, vol. 18, no. 1, pp. 156–185, Jan. 2007.
  * [15] P. Biswas, T.-C. Liang, K.-C. Toh, Y. Ye, and T.-C. Wang,  “Semidefinite Programming Approaches for Sensor Network Localization With Noisy Distance Measurements,”  IEEE Transactions on Automation Science and Engineering, vol. 3, no. 4, pp. 360–371, Oct. 2006.
  * [16] Christoph Helmberg, Franz Rendl, Robert J. Vanderbei, and Henry Wolkowicz,  “An Interior-Point Method for Semidefinite Programming,”  SIAM Journal on Optimization, vol. 6, no. 2, pp. 342–361, May 1996.
  * [17] Steven M. Kay,  Fundamentals of statistical signal processing,  Prentice Hall signal processing series. Prentice-Hall PTR, Englewood Cliffs, N.J, 1993.
  * [18] Haotian Jiang, Tarun Kathuria, Yin Tat Lee, Swati Padmanabhan, and Zhao Song,  “A faster interior point method for semidefinite programming,”  in 2020 IEEE 61st Annual Symposium on Foundations of Computer Science (FOCS), 2020, pp. 910–918.


# Continuous Hamiltonian dynamics on noisy digital quantum computers without
Trotter error

Etienne Granet and Henrik Dreyer  Quantinuum, Leopoldstrasse 180, 80804
Munich, Germany

###### Abstract

We introduce an algorithm to compute Hamiltonian dynamics on digital quantum
computers that requires only a finite circuit depth to reach an arbitrary
precision, i.e. achieves zero Trotter error with finite depth. This finite
number of gates comes at the cost of an attenuation of the measured
expectation value by a known amplitude, requiring more shots per circuit. The
algorithm generalizes to time-dependent Hamiltonians, appearing for example in
adiabatic state preparation. This makes it particularly suitable for present-
day relatively noisy hardware that supports only circuits with moderate depth.

_Introduction._ — Hamiltonian dynamics is one of the most promising
application of current and near-term quantum computers Feynman (2018); Lloyd
(1996). It is believed to be a task that can be performed exponentially faster
on quantum computers and that finds multiple applications, either directly or
indirectly as subroutines of other quantum algorithms Abrams and Lloyd (1999);
Harrow  _et al._ (2009); Farhi  _et al._ (2000); Yang  _et al._ (2023); Lu
_et al._ (2021); Schuckert  _et al._ (2022). Given a Hamiltonian $H$ and a
simulation time $t$, there exist several ways to implement the time evolution
operator $U(t)=e^{itH}$ on a digital quantum computer where only one or two-
qubit gates can be used. The simplest approach uses product formulas such as
Trotterization to decompose the dynamics into a sequence of elementary gates
$U(t)\approx e^{itH_{1}}...e^{itH_{n}}$ Lloyd (1996); Suzuki (1990, 1991).
These methods have been generalized and improved in many ways Berry  _et al._
(2007); Childs and Su (2019); Childs  _et al._ (2021); Low  _et al._ (2023);
Zeng  _et al._ (2022); Cho  _et al._ (2022); Tomesh  _et al._ (2021); Aharonov
and Ta-Shma (2003); Childs and Wiebe (2012); Leadbeater  _et al._ (2023);
Fitzpatrick  _et al._ (2021); Mc Keever and Lubasch (2023). In particular, the
introduction of randomness Campbell (2017); Childs  _et al._ (2019); Faehrmann
_et al._ (2022); Wan  _et al._ (2022); Wang  _et al._ (2023); Gong  _et al._
(2023) such as in the qDRIFT algorithm Campbell (2019); Ouyang  _et al._
(2020); Chen  _et al._ (2021); Nakaji  _et al._ (2023); Kiss  _et al._ (2023);
Hagan and Wiebe (2022); Huggins  _et al._ (2021); Berry  _et al._ (2020);
Pocrnic  _et al._ (2023) have particularly good theoretical performance for
non-sparse Hamiltonians $H$. These algorithms have in common that they only
_approximate_ the continuous time evolution with finite-depth circuits, and
display discretization errors (sometimes called _Trotter errors_) that vanish
only in the limit of infinitely many gates. This is particularly problematic
for adiabatic state preparation (since Trotter errors generically lead to an
effective heating of the system) Aspuru-Guzik  _et al._ (2005); Farhi  _et
al._ (2000); Albash and Lidar (2018) and in applications like quantum
chemistry Su  _et al._ (2021); Reiher  _et al._ (2017); Bauer  _et al._
(2020); McArdle  _et al._ (2020); Motta and Rice (2022); Kim  _et al._ (2022)
where extreme precision is desired. More sophisticated algorithms have better
theoretical scalings, such as linear combination of unitaries Childs and Wiebe
(2012); Berry  _et al._ (2015, 2014), quantum signal processing Low and Chuang
(2017, 2019); Low and Wiebe (2018); Low  _et al._ (2016); Kikuchi  _et al._
(2023) or quantum walks Berry and Childs (2009); Childs (2010), but require a
large resource overhead. As a consequene, at least in the near term, exact
continuous time dynamics is considered to be the prerogative of non-gate-
based, analog quantum simulators Buluta and Nori (2009); Blatt and Roos
(2012); Aspuru-Guzik and Walther (2012); Bloch  _et al._ (2012); Georgescu
_et al._ (2014).

In this Letter we introduce an algorithm for the computation of Hamiltonian
dynamics of observables $\langle\mathcal{M}(t)\rangle=\langle
0|e^{iHt}\mathcal{M}e^{-iHt}|0\rangle$ on digital quantum computers that has
no Trotter error even with a finite number of gates, for any hermitian
Hamiltonian $H$ and unitary $\mathcal{M}$. The expectation value
$\langle\mathcal{M}(t)\rangle$ is obtained as the average over random circuits
drawn from a judiciously chosen distribution, multiplied by a known
amplification factor. The angle $\tau$ of the gates entering the circuit is a
free parameter of our algorithm, that only modifies the amplification factor
and not the precision of the result. It allows one to tune the number of gates
in the circuit at the cost of an increased number of shots on a quantum
computer, which makes it particularly suitable for present-day hardware that
supports only moderate circuit depths. The algorithm straightforwardly
generalizes to time-dependent Hamiltonians, achieving again zero Trotter error
with finite circuit depth. We demonstrate the algorithm with numerical
simulations on a 2D Ising model and for the electronic structure of the
stretched water molecule.

$U_{X}$$U_{X}$$U_{ZZ}$$U_{ZZ}$$U_{ZZ}$$U_{ZZ}$$U_{X}$$U_{X}$$U_{ZZ}$$U_{X}$$U_{X}^{-1}$$U_{X}$$U_{X}$$U_{ZZ}^{-1}$$U_{X}$$U_{ZZ}$$\mathcal{M}$$U_{ZZ}^{-1}$$U_{ZZ}^{-1}$$U_{X}^{-1}$$U_{X}^{-1}$$U_{X}^{-1}$$U_{ZZ}^{-1}$$U_{ZZ}^{-1}$$U_{ZZ}^{-1}$$U_{X}^{-1}$$U_{X}^{-1}$
Figure 1: Example of a tetris configuration for a transverse-field Ising
Hamiltonian, with $U_{O}=e^{i\tau O}$ acting on qubits at the bottom of the
drawing. The average over such configurations converges to the exact time-
evolved expectation value $\langle 0|e^{iHt}\mathcal{M}e^{-iHt}|0\rangle$, up
to a known multiplicative constant.

_Simulating small-angle gates with large-angle gates._ —Let us first introduce
the basic mechanism that underpins our algorithm. For a real number $0\leq
p\leq 1$, an angle $\tau$ and $O$ any operator such that $O^{2}=I$, we have
the equality

$1-p+pe^{i\tau O}=\lambda e^{i\tau^{\prime}O}\,,$ (1)

where

$\tan\tau^{\prime}=\frac{p\sin\tau}{1-p+p\cos\tau}\,,\qquad\lambda=\frac{p\sin\tau}{\sin\tau^{\prime}}\,.$
(2)

This relation can be used to realize an effective angle $0<\tau^{\prime}<\tau$
by applying $e^{i\tau O}$ with probability

$p=\frac{\tan\tau^{\prime}}{\sin\tau+(1-\cos\tau)\tan\tau^{\prime}}\,.$ (3)

To see this, let us consider a circuit containing $G$ gates
$e^{i\tau^{\prime}O}$ that produces the wave function $|\psi\rangle$. For a
subset $S$ of these gates, we denote $|\psi_{S}\rangle$ the wave function
obtained with the same circuit but with deleting the gates
$e^{i\tau^{\prime}O}$ contained in $S$ and replacing the remaining gates
$e^{i\tau^{\prime}O}$ not in $S$ by $e^{i\tau O}$. Then, for any observable
$\mathcal{M}$, by repeatedly applying (1), we have

$\sum_{S,S^{\prime}}(1-p)^{|S|+|S^{\prime}|}p^{2G-|S|-|S^{\prime}|}\langle\psi_{S^{\prime}}|\mathcal{M}|\psi_{S}\rangle=\lambda^{2G}\langle\psi|\mathcal{M}|\psi\rangle\,,$
(4)

where the sum runs over all the possible subsets $S,S^{\prime}$ of gates to
delete in the original circuit. The expectation value
$\langle\psi|\mathcal{M}|\psi\rangle$ that originally involves a circuit with
angle $\tau^{\prime}>0$ can thus be exactly computed with a circuit involving
only larger angles $\tau>\tau^{\prime}$ but in which some gates are randomly
removed. This however comes at the cost of an attenuation factor
$\lambda^{2G}<1$ which increases the total number of shots required to
evaluate $\langle\psi_{S^{\prime}}|\mathcal{M}|\psi_{S}\rangle$ for different
$S,S^{\prime}$ to reach a given precision on
$\langle\psi|\mathcal{M}|\psi\rangle$. The same reasoning applies to gates
with a negative angle $\tau^{\prime}<0$, with then $\tau<\tau^{\prime}$.

_Continuous time limit._ —We now consider a Hamiltonian $H$ that we write as a
linear combination of $N$ operators $O_{n}$ that satisfy $O_{n}^{2}=I$, i.e.

$H=\sum_{n=1}^{N}c_{n}O_{n}\,,$ (5)

which is always possible since a Hamiltonian can be written as a sum of
products of Pauli matrices. We would like to compute the expectation value of
an observable $\mathcal{M}$ within the wave function
$|\psi(t)\rangle=e^{iHt}|\psi(0)\rangle$ obtained by continuous time
evolution. According to the Trotter-Suzuki formula, we can write

$e^{itH}=\underset{\tau^{\prime}\to
0^{+}}{\lim}(e^{i\tau^{\prime}c_{1}O_{1}}...e^{i\tau^{\prime}c_{N}O_{N}})^{t/\tau^{\prime}}\,.$
(6)

For each $n=1,...,N$, we choose $0<\tau_{n}\leq\pi/2$ and implement
$e^{i\tau^{\prime}c_{n}O_{n}}$ using the previously explained protocol with
angle $\tau_{n}$. In the sequence of gates appearing on the right-hand side of
(6), instead of each gate $e^{i\tau^{\prime}c_{n}O_{n}}$ we thus apply
$e^{i\tau_{n}{\rm sgn}\,(c_{n})O_{n}}$ with probability $p_{n}$ given by

$p_{n}=\frac{\tau^{\prime}|c_{n}|}{\sin\tau_{n}}+\mathcal{O}((\tau^{\prime})^{2})\,.$
(7)

The corresponding attenuation factor $\lambda_{n}$ is

$\lambda_{n}=1-\tau^{\prime}|c_{n}|\tan(\tau_{n}/2)+\mathcal{O}((\tau^{\prime})^{2})\,.$
(8)

In the limit $\tau^{\prime}\to 0$, the probability of picking a gate $p_{n}\to
0$ becomes a rare event and we converge to a Poisson process. Namely, for each
gate $e^{i\tau_{n}{\rm sgn}\,(c_{n})O_{n}}$ we obtain a sequence of gate times
$0<t_{n}^{(1)}<...<t_{n}^{(m_{n})}<t$ drawn from a Poisson process with rate
$|c_{n}|/\sin\tau_{n}$. For a given realization $T$ of all these times for
different $n$’s, we denote $|\psi_{T}\rangle$ the wave function obtained by
applying the gates $e^{i\tau_{n}{\rm sgn}\,(c_{n})O_{n}}$ ordered by time of
occurrence $t_{n}^{(i)}$, irrespective of $n$. We call each of these
configurations a _tetris_ , as each gate “falls” randomly on the initial state
with rate $|c_{n}|/\sin\tau_{n}$, which is reminiscent of the eponymous game,
see Fig 1. Then, denoting $\mathbb{E}$ the expectation value with respect to
two independent tetrises $T,T^{\prime}$, we have

$\langle\psi(t)|\mathcal{M}|\psi(t)\rangle=\frac{\mathbb{E}[\langle\psi_{T^{\prime}}|\mathcal{M}|\psi_{T}\rangle]}{\lambda_{\rm
att}}\,,$ (9)

with the total attenuation

$\lambda_{\rm
att}=\exp\left(-2t\sum_{n=1}^{N}|c_{n}|\tan(\tau_{n}/2)\right)\,.$ (10)

This is the fundamental relation that defines our algorithm. Remarkably, the
exact expectation value at time $t$ under continuous time dynamics can be
obtained _without_ Trotter error, by only applying gates with application time
$\tau_{n}>0$ fixed. Moreover, the average number of gates in each circuit is
$2t\sum_{n=1}^{N}\frac{|c_{n}|}{\sin\tau_{n}}$, which is finite and only
controlled by $\tau_{n}$ and not by the precision wanted on the result.

_Statement of the algorithm._ — We state here the algorithm deduced from the
previous explanations. It takes as input a Hamiltonian $H$ decomposed as in
(5), a unitary observable $\mathcal{M}$, a time $t$ and an initial state
$|\psi(0)\rangle$, and outputs $\langle\psi(t)|\mathcal{M}|\psi(t)\rangle$. It
takes $0<\tau_{n}\leq\pi/2$ as $N$ parameters.

  1. 1.

Prepare the system in the initial state.

  2. 2.

For each $n\in\\{1,...,N\\}$, draw an integer $m_{n}$ from a Poisson
distribution with parameter $\frac{|c_{n}|t}{\sin|\tau_{n}|}$. Then draw
$m_{n}$ real numbers $t_{n}^{(i)}$, where $i\in\\{1,...,m_{n}\\}$, uniformly
at random between $0$ and $t$. These $M\equiv\sum_{n=1}^{N}m_{n}$ numbers are
collected into a set $T$.

  3. 3.

Deduce then the sequence $k_{1},...,k_{M}\in\\{1,...,N\\}$ such that the index
$n$ corresponding to the $j$-th smallest element of $T$ is $k_{j}$. For
$m=1,...,M$ in this order, apply the gate $e^{i\tau_{k_{m}}{\rm
sgn}\,(c_{k_{m}})O_{k_{m}}}$ on the system.

  4. 4.

Apply the observable $\mathcal{M}$ on the system.

  5. 5.

Repeat steps (2) and (3) with $\tau_{n}$ replaced by $-\tau_{n}$.

  6. 6.

Measure the overlap with the initial state and divide the result by
$\lambda_{\rm att}$ given in (10). On a quantum computer, this step requires
to perform a Hadamard test.

  7. 7.

Repeat steps (1) to (6) a number of times and average the results.

The algorithm shares similarities with qDRIFT Campbell (2019), as they are
both random compilation algorithms. However the algorithm we present in this
Letter does not have any Trotter error and any angle $\tau_{n}$ can be used.
Our algorithm also requires to perform a backward propagation, i.e. to
implement $e^{itH}$ and $e^{-itH}$ explicitly in the circuit. The drawing
process of the random gates also differ: the tetrises have to be drawn with
the Poisson process described earlier, whereas using the same drawing process
as in qDRIFT would lead to a systematic bias at finite $t,N,\tau$.

_Noisy gates and optimal angles._ — In our algorithm, the gate angles
$\tau_{n}$ are parameters that can be chosen freely and do not deteriorate the
precision. We will now show that there exist optimal $\tau_{n}$’s that only
depend on the level of noise of the quantum computer on which the algorithm is
run. Let us assume that the application of each gate $e^{i\tau_{n}{\rm
sgn}\,(c_{n})O_{n}}$ comes with a signal attenuation $e^{-r_{n}}$ due to
depolarizing noise Yang  _et al._ (2023). This scenario is realistic in state-
of-the-art devices and typical two-qubit gate fidelities are $e^{-r}\sim
99.8\%$, in which case $r\sim 2\mathrm{e}{-3}$ Moses and Others (2023).

There are on average $\frac{t|c_{n}|}{\sin\tau_{n}}$ such gates per tetris.
The damping due to hardware imperfection is thus

$q_{\rm
att}=\exp\left(-2t\sum_{n=1}^{N}\frac{r_{n}|c_{n}|}{\sin\tau_{n}}\right)\,.$
(11)

The total attenuation of the signal is $\lambda_{\rm att}q_{\rm att}$.
Assuming $r_{n}$ small, the optimal times $\tau_{n}^{*}$ that minimize the
total attenuation are

$\tau_{n}^{*}=\sqrt{2r_{n}}\,.$ (12)

Then the minimal total attenuation is

$\lambda_{\rm att}q_{\rm
att}=\exp\left(-2t\sum_{n=1}^{N}\sqrt{2r_{n}}|c_{n}|\right)\,.$ (13)

We note the scaling in $t\sqrt{r_{n}}$ for continuous time dynamics up to time
$t$, whereas the exponent for discrete time dynamics with $t$ steps would
scale as $tr_{n}$.

$0$$1$$2$$-0.5$$0$$0.5$$1$$t$$\langle Z\rangle$

$10^{3}$$10^{4}$$10^{5}$$0.13$$0.14$$0.15$number of samples

Figure 2: Expectation value of $Z$ in the 2D Ising model at $h=3$ in size
$3\times 4$ computed with (9). Left panel: as a function of time $t$, with
gate angles $\tau=0.04$, with $100$ tetrises (teal) and $1000$ tetrises
(orange). For this value of $\tau$, there are around $3000\times t$ gates
$e^{i\tau O}$ per circuit. Right panel: as a function of the number of
samples, for $t=1$ and $\tau=0.04$ (purple) and $\tau=0.08$ (olive). The exact
value is indicated by the black continuous curves in both panels. In the right
panel, the dashed black line indicates the value obtained with a Trotter
decomposition with a time step $\tau=0.04$.

_Scaling comparison with Trotter._ — We now compare the efficiency of our
algorithm with Trotter-Suzuki. Specifically, we wish to estimate the number of
shots $M$ required in our case and in Trotter to obtain a given precision
$\epsilon$ on the result. For simplicity, we assume a depolarizing noise
model, causing a constant signal attenuation $e^{-r}$ for all gates. In our
algorithm, the error after $M$ shots is $\frac{1}{\sqrt{M\lambda_{\rm
att}q_{\rm att}}}$. At the optimal gate time $\tau^{*}=\sqrt{2r}$ we obtain
thus the following required number of shots $M_{\text{Tetris}}$ for a
precision $\epsilon$

$M_{\text{Tetris}}\sim\frac{e^{4t\sqrt{2r}\sum_{n=1}^{N}|c_{n}|}}{\epsilon^{2}}\,.$
(14)

One notes that the precision $\epsilon$ appears only as a multiplicative
factor $1/\epsilon^{2}$, and not in the exponential. In the Trotter algorithm,
any choice of Trotter step $\tau$ comes with a certain minimal error, that
scales as $Nt^{2}C\tau$ for the first order Trotter formula, with some
coefficient $C$. To reach precision $\epsilon$, we thus necessarily have to
take $\tau<\epsilon/(Nt^{2}C)$, hence run a circuit with more than
$N^{2}t^{3}C/\epsilon$ gates. Then the imperfections in the hardware incur a
signal attenuation of order $e^{-N^{2}t^{3}Cr/\epsilon}$. The number of shots
thus has to be at least

$M_{\text{Trotter}}\sim\frac{e^{2N^{2}t^{3}Cr/\epsilon}}{\epsilon^{2}}\,.$
(15)

Higher-order Trotter formulas achieve better scalings with fractional powers
of $\epsilon$ appearing in the exponential, but are generally too costly to be
implemented on present-day hardware. Comparing with the exponent in (14), we
thus see that under the assumption of a depolarizing noise model our algorithm
requires exponentially fewer shots for a given precision $\epsilon$ when

$\epsilon<\frac{t^{2}N^{2}C\sqrt{r}}{2\sqrt{2}\sum_{n=1}^{N}|c_{n}|}\,.$ (16)

Alternatively, without assuming a depolarizing noise model, we can say that
the noise sets an upper limit $\sim 1/r$ for the circuit depth on the
hardware. This sets a lower limit for the precision $\epsilon\sim Nt^{3}Cr$
that can be achieved with a Trotter decomposition, whereas in our algorithm
arbitrary precision can be reached by increasing the number of shots.

_Time-dependent Hamiltonian._ — The algorithm can be generalized
straightforwardly to dynamics with a time-dependent Hamiltonian $H(t)$.
Decomposing

$H(t)=\sum_{n=1}^{N}c_{n}(t)O_{n}\,,$ (17)

with time-dependent coefficients $c_{n}(t)$, and denoting
$|\psi(t)\rangle=\mathcal{T}\exp\left(i\int_{0}^{t}\text{d}sH(s)\right)|\psi(0)\rangle$,
the expectation value of operators $\mathcal{M}$ can again be written as an
average over tetrises (9). The tetrises $|\psi_{T}\rangle$ are produced by
drawing times $0<t_{n}^{(1)}<...<t_{n}^{(m_{n})}<t$ from a Poisson process
with time-dependent rate $|c_{n}(s)|/\sin\tau_{n}$, for each $n=1,...,N$, and
applying on the initial wave function each gate $e^{i\tau_{n}{\rm
sgn}\,(c_{n}(s))O_{n}}$ ordered by time of occurrence $s=t_{n}^{(i)}$. The
attenuation factor is then obtained by replacing $t|c_{n}|$ by
$\int_{0}^{t}\text{d}s|c_{n}(s)|$.

In practice, we proceed as follows. We introduce the function

$z_{n}(u)=\int_{0}^{u}\text{d}s|c_{n}(s)|\,,$ (18)

as well as $z_{n}^{-1}(u)$ its reciprocal, i.e. such that
$z_{n}^{-1}(z_{n}(u))=u$. We replace step (2) of the time-independent case by

  1. 2.

For each $n\in\\{1,...,N\\}$, draw an integer $m_{n}$ from a Poisson
distribution with parameter $\frac{z_{n}(t)}{\sin|\tau_{n}|}$. Then draw
$m_{n}$ real numbers $\tilde{t}_{n}^{(i)}$, where $i\in\\{1,...,m_{n}\\}$,
uniformly at random between $0$ and $z_{n}(t)$, and set
$t_{n}^{(i)}=z_{n}^{-1}(\tilde{t}_{n}^{(i)})$. These
$M\equiv\sum_{n=1}^{N}m_{n}$ numbers $t_{n}^{(i)}$ are collected into a set
$T$.

Then in step (3), the gate applied is $e^{i\tau_{k_{m}}{\rm
sgn}\,(c_{k_{m}}(s))O_{k_{m}}}$ where $s\in\\{t_{n}^{(i)}\\}$ is the time at
which the gate is applied. For the backward propagation one has to reverse the
order of the gates. The total attenuation $\lambda_{\rm att}$ of step (6) is

$\lambda_{\rm
att}=\exp\left(-2\sum_{n=1}^{N}z_{n}(t)\tan(\tau_{n}/2)\right)\,.$ (19)

This algorithm produces the _exact_ time evolution of the wave function under
the time-dependent Hamiltonian $H(t)$, without errors arising from
discretizing the continuous function $c_{n}(s)$ or Trotter errors.

$0$$1$$2$$-2.8$$-2.6$$-2.4$$-2.2$$-2$$T_{f}$energy per site

$3$$3.5$$4$$\cdot 10^{-2}$$0$$0.5$$1$$1.5$$2$$t$$R(t)$

Figure 3: Left: _Adiabatic state preparation._ Energy per site at time
$t=T_{f}$ as a function of $T_{f}$, for the Ising model in size $3\times 4$
with the time-dependent magnetic field $h(t)$ quoted in the text, starting
from the $Z=1$ product state: noiseless simulations with $10^{5}$ tetrises and
$\tau=0.04$ (teal, the shade indicating one standard deviation), and exact
value (black). Right: _Chemistry Hamiltonian._ Ratio $R(t)$ in (23) as a
function of $t$, exact value (black), noiseless simulations (teal) and noisy
simulations (orange) with a two-body gate fidelity $0.9982$. In all cases the
angles are given by (12). Each run of the algorithm uses $2000$ tetrises and
$20$ shots per tetris, and each tetris has in average $144$ two-body gates.

_Background Hamiltonian._ — The algorithm has the following simple
improvement. Let us consider a subset of gates
$\mathcal{O}\subset\\{O_{n}\\}_{n=1,..,N}$ that all commute among themselves,
i.e. $[O,O^{\prime}]=0$ for all $O,O^{\prime}\in\mathcal{O}$, and denote
$H_{\rm background}=\sum_{O_{p}\in\mathcal{O}}c_{p}O_{p}$. If for these gates
we take the limit of zero angle $\tau_{p}\to 0$ in the algorithm, the Poisson
process almost always picks gates in $\mathcal{O}$ with very short angle. But
since they commute this can be implemented with a finite product of $e^{i\tau
O}$ with $\tau$ obtained by summing the different angles. Hence the algorithm
is modified as follows. In step (2) we only draw times corresponding to the
gates _not_ in $\mathcal{O}$, and use them to produce an ordered sequence of
gates as in step (3). However, between each application of these gates
$e^{i\tau_{n}O_{n}}$ and $e^{i\tau_{m}O_{m}}$ at times $s_{n}<s_{m}$, we apply
$e^{i\delta sH_{\rm background}}$ with $\delta s=s_{m}-s_{n}$ the difference
between the two instants at which the gates are applied. This reduces the
attenuation factor $\lambda_{\rm att}$, which is now computed _only_ including
the gates _not_ in $\mathcal{O}$. The Hamiltonian $H_{\rm background}$ thus
plays the role of a “background” Hamiltonian that is always applied on the
system, but on top of which gates $e^{i\tau_{n}O_{n}}$ with
$O_{n}\notin\mathcal{O}$ are applied.

_Numerical tests: square lattice Ising model._ — We present numerical tests of
our algorithm on the 2D square lattice Ising model in a transverse field $h$,
with Hamiltonian

$H=-\sum_{\langle i,j\rangle}Z_{i}Z_{j}-h\sum_{j}X_{j}\,,$ (20)

with periodic boundary conditions. In Figure 2 we show $\langle Z(t)\rangle$,
with the initial state being the $Z=+1$ product state, as computed with
formula (9). We show the convergence of (9) towards the exact value in the
limit of large number of tetrises, even if the gate time $\tau$ is finite. The
convergence is faster when the angle $\tau$ is small, but the corresponding
quantum circuits involve then more gates.

We then consider the time-dependent case of our algorithm with the adiabatic
preparation of the ground state at $h_{f}=2.5$, starting from $h=0$. For a
final time $T_{f}$ we consider the time-dependent magnetic field
$h(t)=h_{f}\sin(\tfrac{\pi}{2}\sin(\tfrac{\pi t}{2T_{f}})^{2})^{2}$ Wei  _et
al._ (2023). In the left panel of Figure 3 we show the final energy per site
obtained at $t=T_{f}$, as a function of $T_{f}$, together with the error bars
for a fixed number of shots. We see the agreement with the exact values, as
well as the broadening of the error bars as $T_{f}$ grows at fixed $\tau$, as
imposed by (10).

_Noisy numerical tests: quantum chemistry._ — We now consider a noisy
application of our algorithm on a non-sparse Hamiltonian. Hamiltonians
describing molecular electronic structure problems are typically decomposed in
a basis set of orbitals, where they read

$H=\sum_{ij=1}^{N}h_{ij}c^{\dagger}_{i}c_{j}+\sum_{ijkl=1}^{N}h_{ijkl}c^{\dagger}_{i}c^{\dagger}_{j}c_{k}c_{l}\,,$
(21)

with $c$’s canonical fermionic operators and $N$ the number of orbitals
considered. Written in terms of Pauli gates through a Jordan-Wigner
transformation, these molecular Hamiltonians involve thus $\mathcal{O}(N^{4})$
terms with each $\mathcal{O}(N)$ Pauli gates, incurring a prohibitively large
gate count for any Trotter implementation of the dynamics.

We consider a H2O water molecule in a STO-6G basis with $14$ orbitals, and use
our algorithm to study the Loschmidt echo

$\mathcal{L}(t)=\langle HF|e^{itH}|HF\rangle\,,$ (22)

where $|HF\rangle$ denotes the Hartree-Fock ground state. We take a geometry
with an angle H-O-H of $105^{o}$ and with an elongated distance H-O of $2.2$
Å. This geometry is chosen so as to have a more significant difference between
the exact ground state and $|HF\rangle$. Our algorithm for expectation values
is straightforwardly modified to compute $\mathcal{L}(t)$ by skipping steps
(4) and (5) and using $\sqrt{\lambda_{\rm att}}$ in step (6), as we have
$\mathcal{L}(t)=\mathbb{E}[\langle HF|\psi_{T}\rangle]/\sqrt{\lambda_{\rm
att}}$.

On a noisy hardware, the attenuation due to $q_{\rm att}$ has to be mitigated.
One strength of our algorithm is that one can reduce the depth of the
circuits, and so the noise, by increasing the angles $\tau$, without affecting
the precision. Alternatively, one may compute the ratio

$R(t)=\frac{\Im\mathcal{L}(t)}{\Re\mathcal{L}(t)}\,.$ (23)

This ratio should be noise-free if the noise on the hardware can be well
approximated by a depolarizing channel, while in principle containing
information about the energy of all states that have non-zero overlap with
$\ket{HF}$, including the ground state, as they appear in the oscillation
frequencies of the numerator and denominator. In the right panel of Figure 3
we show a noisy simulation of our algorithm that includes a depolarizing
channel after each $2$-qubit gate, and observe excellent agreement.

A few comments on a hardware implementation of our algorithm are in order.
Implementing (9) with an actual digital quantum circuit requires to perform a
Hadamard test, as one has to average the amplitude
$\langle\psi_{T^{\prime}}|\mathcal{M}|\psi_{T}\rangle$ over the tetrises
$T,T^{\prime}$, and not the absolute value square. One also has to choose a
certain number of shots $M$ to do per tetris, i.e. how to distribute the
resources between number of tetrises and accuracy of each tetris. Neglecting
compilation time, the optimal number of shots can be seen to be $M=1$.
Separately, the numerical data reported in this Letter underestimates the
efficiency of the algorithm, since no background Hamiltonian techniques were
used. A more quantitative study of its efficiency would be interesting to
conduct.

_Interpolation between quantum time evolution and path integral
representation._ — Some comments about the upper and lower limits of the angle
$\tau$ are in order. When $\tau\to 0$, the different gates $e^{i\tau O_{n}}$
commute at order $\tau$ and there are increasingly more gates drawn from the
Poisson process. With probability $1$ in this limit, a single tetris
$|\psi_{T}\rangle$ is equal to the exact time evolution $e^{iHt}|0\rangle$,
and the attenuation factor (10) is equal to $1$. When $\tau=\pi/2$, we have
$e^{i\tau O_{n}}=iO_{n}$, and so when $O_{n}$ is a string of Pauli operators,
the gates can be implemented with only one-qubit gates. Each tetris
$|\psi_{T}\rangle$ remains thus a product state in the $Z$ basis and describes
a classical trajectory in discrete time. Thus (9) becomes akin to a path
integral representation, in the sense that a quantum expectation value is
expressed as an averaging over exponentially many classical trajectories.

These two limits provide an interesting intepretation of our algorithm as an
_interpolation_ between a single quantum trajectory for $\tau=0$, and
exponentially many classical trajectories for $\tau=\pi/2$. Intermediate
values of $\tau$ near $\pi/2$ allow one to decrease the entanglement present
in each tetris $|\psi_{T}\rangle$ compared to the exact time-evolved state
$e^{iHt}|0\rangle$, but at the cost of increasing the number of samples
required through the attenuation factor (10). This could find applications in
the classical simulation of quantum dynamics, where the major obstacle is the
entanglement in $e^{iHt}|0\rangle$ that grows linearly with time. One could
instead implement (9) classically by choosing $\tau>0$ such that each tetris
$|\psi_{T}\rangle$ has little enough entanglement to be implemented with
tensor network techniques. This would come at the cost of summing over
multiple samples because of the attenuation factor (10), but this sampling can
be efficiently parallellized.

_Discussion._ — We presented a quantum algorithm to compute continuous
Hamiltonian dynamics that requires a circuit depth that is independent of the
desired precision $\epsilon$, contrary to previously known algorithms that
require an infinite depth when $\epsilon\to 0$. We show that by applying
randomly gates on the initial state according to a well-chosen distribution,
and multiplying the result by a known amplification factor, one achieves zero
Trotter error while having finite-depth circuits. This amplification factor
however increases the number of shots required to reach a given precision.

The algorithm and its mechanisms can be generalized in a number of ways. For
example, decomposing a circuit into Clifford gates and $T$-gates, one can
implement each $T$-gate by applying a $e^{iZ\pi/4}$-gate or no gate both with
probability $1/2$, which is a $S$-gate times a global phase. Because of the
attenuation factor $\lambda=\cos(\pi/8)$ for each replaced $T$-gate, this
enables one to sample classically from a circuit with $t$ many $T$-gates in a
time $e^{2|\log(\cos\tfrac{\pi}{8})|t}$ multiplied by the simulation time of
the resulting Clifford circuits. The exponent matches the currently best known
exponent Bravyi and Gosset (2016).

We thank David Zsolt Manrique for helpful comments on the draft. E.G.
acknowledges support by the Bavarian Ministry of Economic Affairs, Regional
Development and Energy (StMWi) under project Bench-QC (DIK0425/01).

## References

  * Feynman (2018) R. P. Feynman, in  _Feynman and computation_ (CRC Press, 2018) pp. 133–153.
  * Lloyd (1996) S. Lloyd, Science 273, 1073 (1996).
  * Abrams and Lloyd (1999) D. S. Abrams and S. Lloyd, Physical Review Letters 83, 5162 (1999).
  * Harrow  _et al._ (2009) A. W. Harrow, A. Hassidim, and S. Lloyd, Physical review letters 103, 150502 (2009).
  * Farhi  _et al._ (2000) E. Farhi, J. Goldstone, S. Gutmann, and M. Sipser, arXiv preprint quant-ph/0001106 (2000), 10.48550/arXiv.quant-ph/0001106.
  * Yang  _et al._ (2023) Y. Yang, A. Christianen, S. Coll-Vinent, V. Smelyanskiy, M. C. Bañuls, T. E. O’Brien, D. S. Wild, and J. I. Cirac, arXiv preprint arXiv:2303.08461 (2023), https://doi.org/10.48550/arXiv.2303.08461.
  * Lu  _et al._ (2021) S. Lu, M. C. Bañuls, and J. I. Cirac, PRX Quantum 2, 020321 (2021), https://doi.org/10.1103/PRXQuantum.2.020321.
  * Schuckert  _et al._ (2022) A. Schuckert, A. Bohrdt, E. Crane, and M. Knap, Phys. Rev. B 107, L140410 (2022), https://doi.org/10.1103/PhysRevB.107.L140410.
  * Suzuki (1990) M. Suzuki, Physics Letters A 146, 319 (1990).
  * Suzuki (1991) M. Suzuki, Journal of Mathematical Physics 32, 400 (1991).
  * Berry  _et al._ (2007) D. W. Berry, G. Ahokas, R. Cleve, and B. C. Sanders, Communications in Mathematical Physics 270, 359 (2007).
  * Childs and Su (2019) A. M. Childs and Y. Su, Physical review letters 123, 050503 (2019).
  * Childs  _et al._ (2021) A. M. Childs, Y. Su, M. C. Tran, N. Wiebe, and S. Zhu, Physical Review X 11, 011020 (2021).
  * Low  _et al._ (2023) G. H. Low, Y. Su, Y. Tong, and M. C. Tran, PRX Quantum 4, 020323 (2023).
  * Zeng  _et al._ (2022) P. Zeng, J. Sun, L. Jiang, and Q. Zhao, arXiv preprint arXiv:2212.04566 (2022), 10.48550/arXiv.2212.04566.
  * Cho  _et al._ (2022) C. H. Cho, D. W. Berry, and M.-H. Hsieh, arXiv preprint arXiv:2210.11281 (2022), 10.48550/arXiv.2210.11281.
  * Tomesh  _et al._ (2021) T. Tomesh, K. Gui, P. Gokhale, Y. Shi, F. T. Chong, M. Martonosi, and M. Suchara, in  _2021 International Conference on Rebooting Computing (ICRC)_ (IEEE, 2021) pp. 1–13.
  * Aharonov and Ta-Shma (2003) D. Aharonov and A. Ta-Shma, in  _Proceedings of the thirty-fifth annual ACM symposium on Theory of computing_ (2003) pp. 20–29.
  * Childs and Wiebe (2012) A. M. Childs and N. Wiebe, arXiv preprint arXiv:1202.5822 (2012), 10.26421/QIC12.11-12.
  * Leadbeater  _et al._ (2023) C. Leadbeater, N. Fitzpatrick, D. M. Ramo, and A. J. Thom, arXiv preprint arXiv:2304.07917 (2023), 10.48550/arXiv.2304.07917.
  * Fitzpatrick  _et al._ (2021) N. Fitzpatrick, H. Apel, and D. M. Ramo, arXiv preprint arXiv:2106.03985 (2021), 10.48550/arXiv.2106.03985.
  * Mc Keever and Lubasch (2023) C. Mc Keever and M. Lubasch, Physical Review Research 5, 023146 (2023).
  * Campbell (2017) E. Campbell, Physical Review A 95, 042306 (2017).
  * Childs  _et al._ (2019) A. M. Childs, A. Ostrander, and Y. Su, Quantum 3, 182 (2019).
  * Faehrmann  _et al._ (2022) P. K. Faehrmann, M. Steudtner, R. Kueng, M. Kieferová, and J. Eisert, Quantum 6, 806 (2022).
  * Wan  _et al._ (2022) K. Wan, M. Berta, and E. T. Campbell, Physical Review Letters 129, 030503 (2022).
  * Wang  _et al._ (2023) S. Wang, S. McArdle, and M. Berta, arXiv preprint arXiv:2302.01873 (2023), 10.48550/arXiv.2302.01873.
  * Gong  _et al._ (2023) W. Gong, Y. Kharkov, M. C. Tran, P. Bienias, and A. V. Gorshkov, arXiv preprint arXiv:2307.13028 (2023), 10.48550/arXiv.2307.13028.
  * Campbell (2019) E. Campbell, Physical review letters 123, 070503 (2019).
  * Ouyang  _et al._ (2020) Y. Ouyang, D. R. White, and E. T. Campbell, Quantum 4, 235 (2020).
  * Chen  _et al._ (2021) C.-F. Chen, H.-Y. Huang, R. Kueng, and J. A. Tropp, PRX Quantum 2, 040305 (2021).
  * Nakaji  _et al._ (2023) K. Nakaji, M. Bagherimehrab, and A. Aspuru-Guzik, arXiv preprint arXiv:2302.14811 (2023), 10.48550/arXiv.2302.14811.
  * Kiss  _et al._ (2023) O. Kiss, M. Grossi, and A. Roggero, Quantum 7, 977 (2023).
  * Hagan and Wiebe (2022) M. Hagan and N. Wiebe, arXiv preprint arXiv:2206.06409 (2022), 10.48550/arXiv.2206.06409.
  * Huggins  _et al._ (2021) W. J. Huggins, S. McArdle, T. E. O’Brien, J. Lee, N. C. Rubin, S. Boixo, K. B. Whaley, R. Babbush, and J. R. McClean, Physical Review X 11, 041036 (2021).
  * Berry  _et al._ (2020) D. W. Berry, A. M. Childs, Y. Su, X. Wang, and N. Wiebe, Quantum 4, 254 (2020).
  * Pocrnic  _et al._ (2023) M. Pocrnic, M. Hagan, J. Carrasquilla, D. Segal, and N. Wiebe, arXiv preprint arXiv:2306.16572 (2023), 10.48550/arXiv.2306.16572.
  * Aspuru-Guzik  _et al._ (2005) A. Aspuru-Guzik, A. D. Dutoi, P. J. Love, and M. Head-Gordon, Science 309, 1704 (2005).
  * Albash and Lidar (2018) T. Albash and D. A. Lidar, Reviews of Modern Physics 90, 015002 (2018).
  * Su  _et al._ (2021) Y. Su, D. W. Berry, N. Wiebe, N. Rubin, and R. Babbush, PRX Quantum 2, 040332 (2021).
  * Reiher  _et al._ (2017) M. Reiher, N. Wiebe, K. M. Svore, D. Wecker, and M. Troyer, Proceedings of the national academy of sciences 114, 7555 (2017).
  * Bauer  _et al._ (2020) B. Bauer, S. Bravyi, M. Motta, and G. K.-L. Chan, Chemical Reviews 120, 12685 (2020).
  * McArdle  _et al._ (2020) S. McArdle, S. Endo, A. Aspuru-Guzik, S. C. Benjamin, and X. Yuan, Reviews of Modern Physics 92, 015003 (2020).
  * Motta and Rice (2022) M. Motta and J. E. Rice, Wiley Interdisciplinary Reviews: Computational Molecular Science 12, e1580 (2022).
  * Kim  _et al._ (2022) I. H. Kim, Y.-H. Liu, S. Pallister, W. Pol, S. Roberts, and E. Lee, Physical Review Research 4, 023019 (2022).
  * Berry  _et al._ (2015) D. W. Berry, A. M. Childs, R. Cleve, R. Kothari, and R. D. Somma, Physical review letters 114, 090502 (2015).
  * Berry  _et al._ (2014) D. W. Berry, A. M. Childs, R. Cleve, R. Kothari, and R. D. Somma, in  _Proceedings of the forty-sixth annual ACM symposium on Theory of computing_ (2014) pp. 283–292.
  * Low and Chuang (2017) G. H. Low and I. L. Chuang, Physical review letters 118, 010501 (2017).
  * Low and Chuang (2019) G. H. Low and I. L. Chuang, Quantum 3, 163 (2019).
  * Low and Wiebe (2018) G. H. Low and N. Wiebe, arXiv preprint arXiv:1805.00675 (2018), 10.48550/arXiv.1805.00675.
  * Low  _et al._ (2016) G. H. Low, T. J. Yoder, and I. L. Chuang, Physical Review X 6, 041067 (2016).
  * Kikuchi  _et al._ (2023) Y. Kikuchi, C. M. Keever, L. Coopmans, M. Lubasch, and M. Benedetti, arXiv preprint arXiv:2303.05533 (2023), 10.48550/arXiv.2303.05533.
  * Berry and Childs (2009) D. W. Berry and A. M. Childs, arXiv preprint arXiv:0910.4157 (2009), 10.26421/QIC12.1-2.
  * Childs (2010) A. M. Childs, Communications in Mathematical Physics 294, 581 (2010).
  * Buluta and Nori (2009) I. Buluta and F. Nori, Science 326, 108 (2009).
  * Blatt and Roos (2012) R. Blatt and C. F. Roos, Nature Physics 8, 277 (2012).
  * Aspuru-Guzik and Walther (2012) A. Aspuru-Guzik and P. Walther, Nature physics 8, 285 (2012).
  * Bloch  _et al._ (2012) I. Bloch, J. Dalibard, and S. Nascimbene, Nature Physics 8, 267 (2012).
  * Georgescu  _et al._ (2014) I. M. Georgescu, S. Ashhab, and F. Nori, Reviews of Modern Physics 86, 153 (2014).
  * Moses and Others (2023) S. Moses and Others, arXiv preprint arXiv:2305.03828 (2023), https://doi.org/10.48550/arXiv.2305.03828.
  * Wei  _et al._ (2023) Z.-Y. Wei, D. Malz, and J. I. Cirac, Phys. Rev. Research 5, L022037 (2023), https://doi.org/10.1103/PhysRevResearch.5.L022037.
  * Bravyi and Gosset (2016) S. Bravyi and D. Gosset, Physical review letters 116, 250501 (2016).