Hamiltonian dynamics on digital quantum computers without discretization error

Abstract

We introduce an algorithm to compute expectation values of time-evolved observables on digital quantum computers that requires only bounded average circuit depth to reach arbitrary precision, i.e. produces an unbiased estimator with finite average depth. This finite depth comes with an attenuation of the measured expectation value by a known amplitude, requiring more shots per circuit. The average gate count per circuit for simulation time t is ${\mathcal{O}}({t}^{2}{\mu }^{2})$ with μ the sum of the Hamiltonian coefficients, without dependence on precision, providing a significant improvement over previous algorithms. With shot noise, the average runtime is ${\mathcal{O}}({t}^{2}{\mu }^{2}{\epsilon }^{-2})$ to reach precision ϵ. The only dependence in the sum of the coefficients makes it particularly adapted to non-sparse Hamiltonians. The algorithm generalizes to time-dependent Hamiltonians, appearing for example in adiabatic state preparation. These properties make 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^1,2. It is believed to be exponentially faster to perform on quantum computers and finds multiple applications, either directly or indirectly as subroutines of other quantum algorithms^3,4,5,6,7,8. Given a Hamiltonian H and a simulation time t, there exist several ways to implement 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 elementary gates $U(t)\,\approx\, {e}^{it{H}_{1}}\ldots {e}^{it{H}_{n}}$^2,9,10. These methods have been generalized and improved in many ways^{11,12,13,14,15,16,17,18,19,20,21,22,23,24}. In particular, the introduction of randomness^{25,26,27,28,29,30} such as in the qDRIFT algorithm^{31,32,33,34,35,36,37,38,39} has particularly good theoretical performance for non-sparse Hamiltonians H, i.e. when the number of terms in the Hamiltonian grows faster than ${\mathcal{O}}(L)$ where L is the number of qubits. 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 with infinitely many gates. This is particularly problematic for adiabatic state preparation (since Trotter errors generically lead to an effective heating)^5,40,41 and in applications like quantum chemistry^{42,43,44,45,46,47} where extreme precision is desired. More sophisticated algorithms have better theoretical scalings, such as linear combination of unitaries^19,48,49 (LCU), quantum signal processing^{50,51,52,53,54}, or quantum walks^55,56, but require a large resource overhead. As a consequence, at least in the near term, exact continuous-time dynamics are considered to be the prerogative of non-gate-based, analog quantum simulators^{57,58,59,60,61}.

In this Article we introduce an algorithm to compute Hamiltonian dynamics of observables $\langle {\mathcal{M}}(t)\rangle =\langle 0| {e}^{-iHt}{\mathcal{M}}{e}^{iHt}| 0\rangle$ on digital quantum computers without discretization error even with a finite number of gates, for any Hamiltonian H and unitary ${\mathcal{M}}$. This means that the algorithm produces an unbiased estimator for expectation values of time-evolved observables, with bounded average number of gates. It also applies to amplitude probabilities such as 〈0∣e^iHt∣0〉. 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 factor. One can choose freely the angle τ of the rotations entering the circuit, only modifying the multiplicative factor and not the precision of the result. It allows one to decrease the gate count at the cost of an increased number of shots on a quantum computer, yielding a natural noise-mitigation technique. The optimal choice of gate angle yields a shot overhead of order ${\mathcal{O}}(1)$ with an average gate count per circuit ${\mathcal{O}}({t}^{2}{\mu }^{2})$ where μ is the sum of the coefficients of the Hamiltonian, without any dependence on the precision desired, in contrast with all previous algorithms. Since the gate count per circuit depends only on μ, the algorithm is particularly adapted to non-sparse Hamiltonians. It also straightforwardly generalizes to time-dependent Hamiltonians without any overhead. Taking into account the shot noise, the runtime of our algorithm is thus ${\mathcal{O}}({t}^{2}{\mu }^{2}{\epsilon }^{-2})$ to reach precision ϵ on the expectation value measured. This is in contrast with all previous shot-based algorithms where the dependence in ϵ is at least ${\epsilon }^{-2}\log 1/\epsilon$, due to the gate count per circuit that depends on ϵ. We demonstrate the algorithm with numerical simulations of the electronic structure of the stretched water molecule and on a 2D Ising model where it outperforms both Trotter formulas as well as other randomised compilation techniques.

Results

Simulating small-angle gates with large-angle gates

We first introduce the basic mechanism that underpins our algorithm. For a real number 0 < p ≤ 1, an angle 0 ≤ τ ≤ π/2 and O any operator such that O² = I, we have the equality

$$(1-p)I+p{e}^{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 is obtained from the fact that ${e}^{i\tau O}=I\cos \tau +iO\sin \tau$ for any τ, when O² = I. This relation can be interpreted as follows: the random unitary operator that is equal to the rotation e^iτO with probability p and to identity with probability 1 − p is equal to the (non-random, deterministic) unitary operator ${e}^{i{\tau }^{{\prime} }O}$, up to a rescaling factor λ. This factor λ will be called attenuation factor. Inverting the relation (2), we can instead see τ and ${\tau }^{{\prime} }$ as parameters and deduce the function $p(\tau ,{\tau }^{{\prime} })$ as

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

(3)

Let us check that for $0\le {\tau }^{{\prime} }\le \tau \le \pi /2$, this quantity is between 0 and 1. We have

$$\begin{array}{ll}{\partial }_{{\tau }^{{\prime} }}p\,=\,\frac{1+{\tan }^{2}{\tau }^{{\prime} }}{\sin \tau +(1-\cos \tau )\tan {\tau }^{{\prime} }}\left(1-\frac{(1-\cos \tau )\tan {\tau }^{{\prime} }}{\sin \tau +(1-\cos \tau )\tan {\tau }^{{\prime} }}\right).\end{array}$$

(4)

Since $\sin \tau$ and $(1-\cos \tau )\tan {\tau }^{{\prime} }$ are both positive in this regime, we see that the fraction in the parenthesis is smaller than 1. Hence, ${\partial }_{{\tau }^{{\prime} }}p\ge 0$. Moreover we have $p(\tau ,{\tau }^{{\prime} }=0)=0$ and $p(\tau ,{\tau }^{{\prime} }=\tau )=1$, we conclude that for $0\le {\tau }^{{\prime} }\le \tau \le \pi /2$ we always have 0 ≤ p ≤ 1.

It follows that in this regime $p(\tau ,{\tau }^{{\prime} })$ can be interpreted as a probability. Let us now explain how (1) can be used to express a circuit with small gate angles in terms of circuits with large gate angles drawn randomly, through the following minimal example. We consider the state $\vert \psi \rangle ={e}^{i{\tau }^{{\prime} }O}\vert 0\rangle$ with $\vert 0\rangle$ a given initial state, and an observable ${\mathcal{M}}$, with $0\le {\tau }^{{\prime} }\le \tau$. Using (1), the expectation value of ${\mathcal{M}}$ within $\vert \psi \rangle$ can be written as

$$\begin{array}{ll}\langle \psi | {\mathcal{M}}| \psi \rangle \,=\,{\lambda }^{-2}{(1-p)}^{2}\langle 0| {\mathcal{M}}| 0\rangle \\ \qquad\qquad+\,{\lambda }^{-2}p(1-p)\langle 0| {e}^{-i\tau O}{\mathcal{M}}| 0\rangle +{\lambda }^{-2}p(1-p)\langle 0| {\mathcal{M}}{e}^{i\tau O}| 0\rangle \\ \qquad\qquad+\,{\lambda }^{-2}{p}^{2}\langle 0| {e}^{-i\tau O}{\mathcal{M}}{e}^{i\tau O}| 0\rangle .\end{array}$$

(5)

This equation means that $\langle 0| {e}^{-i{\tau }^{{\prime} }O}{\mathcal{M}}{e}^{i{\tau }^{{\prime} }O}| 0\rangle$ can be obtained in average by deleting with probability 1 − p the rotations e^±iτO in $\langle 0| {e}^{-i\tau O}{\mathcal{M}}{e}^{i\tau O}| 0\rangle$, up to a rescaling factor λ^−2. This factor λ^−2 increases the statistical shot noise when measuring these quantum expectation values on a quantum computer. To reach precision ϵ on $\langle 0| {e}^{-i{\tau }^{{\prime} }O}{\mathcal{M}}{e}^{i{\tau }^{{\prime} }O}| 0\rangle$, one has to do around $1/{(\epsilon {\lambda }^{2})}^{2}$ shots on each of the quantum expectation values on the right-hand side of (5). We note that similar mechanisms have been used recently to implement arbitrary rotation angles using predefined rotation angles⁶².

Let us write more formally this simple example in a general context. We consider a circuit that contains G many rotations ${e}^{i{\tau }^{{\prime} }O}$ and that produces the wave function $\vert \psi \rangle$. For a subset S of these rotations, we denote $\vert {\psi }_{S}\rangle$ the wave function obtained with the same circuit but with deleting the rotations ${e}^{i{\tau }^{{\prime} }O}$ contained in S and replacing the remaining rotations ${e}^{i{\tau }^{{\prime} }O}$ not in S by e^iτ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 ,$$

(6)

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 rotations are randomly removed with probability 1 − p. The attenuation factor λ^2G increases by a factor λ^−4G the number of shots to run on the quantum computer to reach a given precision.

Continuous time limit

Let us now apply the previous mechanism to Hamiltonian simulation. We consider a Hamiltonian H that we write as a linear combination of unitaries

$$H=\mathop{\sum }\limits_{n=1}^{N}{c}_{n}{O}_{n},$$

(7)

with ${O}_{n}^{2}=I$, which is always possible by writing it as a sum of product of Pauli matrices. Up to replacing O_n by − O_n, we can always assume c_n > 0. We would like to compute the expectation value of a unitary ${\mathcal{M}}$ within the wave function $\vert \psi (t)\rangle ={e}^{iHt}\vert \psi (0)\rangle$. To that end, we use the Trotter product formula to write

$${e}^{itH}=\mathop{\lim }\limits_{{\tau }^{{\prime} }\to {0}^{+}}{({e}^{i{\tau }^{{\prime} }{c}_{1}{O}_{1}}\ldots {e}^{i{\tau }^{{\prime} }{c}_{N}{O}_{N}})}^{t/{\tau }^{{\prime} }}.$$

(8)

As ${\tau }^{{\prime} }\to 0$, this formula involves ${\mathcal{O}}(1/{\tau }^{{\prime} })$ rotations ${e}^{i{\tau }^{{\prime} }{c}_{n}{O}_{n}}$ with an angle that goes to 0. It requires circuits increasingly deeper with increasingly smaller angles. Let us apply the previous mechanism to this equation. For each n = 1, …, N, we choose 0 < τ_n ≤ π/2 and implement ${e}^{i{\tau }^{{\prime} }{c}_{n}{O}_{n}}$ using the previously explained protocol with angle τ_n, for ${\tau }^{{\prime} } > 0$ such that ${\tau }^{{\prime} }\le {\tau }_{n}$ for all n. We introduce here a different gate angle τ_n for every term in the Hamiltonian, but for simplicity it could be taken the same for all n. In the sequence of rotations appearing on the right-hand side of (8), instead of each rotation ${e}^{i{\tau }^{{\prime} }{c}_{n}{O}_{n}}$ we apply ${e}^{i{\tau }_{n}{O}_{n}}$ with probability p_n given by

$${p}_{n}=\frac{{\tau }^{{\prime} }{c}_{n}}{\sin {\tau }_{n}}+{\mathcal{O}}({({\tau }^{{\prime} })}^{2}),$$

(9)

and identity with probability 1 − p_n. The corresponding attenuation factor λ_n is

$${\lambda }_{n}=1-{\tau }^{{\prime} }{c}_{n}\tan ({\tau }_{n}/2)+{\mathcal{O}}({({\tau }^{{\prime} })}^{2}).$$

(10)

We are thus writing the time evolution operator e^itH as a statistical expectation value over random unitaries as

$${e}^{itH}=\mathop{\lim }\limits_{{\tau }^{{\prime} }\to 0}\left(\mathop{\prod }\limits_{n\,=\,1}^{N}{\lambda }_{n}^{-M}\right){\mathbb{E}}[{V}^{(1)}\ldots {V}^{(M)}],$$

(11)

where we introduced $M=t/{\tau }^{{\prime} }$ assumed to be an integer, and the random unitary operator V^(m) defined as

$${V}^{(m)}={U}_{1}^{(m)}\ldots {U}_{N}^{(m)}.$$

(12)

Here, ${U}_{n}^{(m)}$ are independent random unitary operators that are equal to I with probability 1 − p_n, and equal to ${e}^{i{\tau }_{n}{O}_{n}}$ with probability p_n. The limit when ${\tau }^{{\prime} }\to 0$ of the product of attenuation factors is given by

$$\mathop{\lim }\limits_{{\tau }^{{\prime} }\to 0}\mathop{\prod }\limits_{n=1}^{N}{\lambda }_{n}^{M}=\exp \left(-t\mathop{\sum }\limits_{n=1}^{N}{c}_{n}\tan ({\tau }_{n}/2)\right)\equiv \Lambda .$$

(13)

In V^(m), for finite ${\tau }^{{\prime} } >\, 0$, there is a non-zero probability that two unitaries ${U}_{n}^{(m)}$ are not identity. We write

$${V}^{(m)}={\tilde{V}}^{(m)}+\delta {V}^{(m)},$$

(14)

with ${\tilde{V}}^{(m)}$ the random unitary operator equal to I with probability $1-{{\prod }_{i}(1-p_i)}\times{{\sum }_{n}{p}_{n}/(1-{p}_{n})}$, and for all n = 1, …, N, equal to ${e}^{i{\tau }_{n}{O}_{n}}$ with probability ${p}_{n}/(1-{p}_{n})\times {{\prod }_{i}(1-p_i)}$. The operator δV^(m) is thus different from 0 with probability

$$\begin{array}{ll}q({\tau }^{{\prime} })\,=\,1-(1-{p}_{1})\ldots (1-{p}_{N})\left(1+\mathop{\sum }\limits_{n=1}^{N}\frac{{p}_{n}}{1-{p}_{n}}\right)\\ \qquad=\,{({\tau }^{{\prime} })}^{2}{\mathop{\sum }\limits_{n < n'}\frac{{c}_{n}c_{n'}}{\sin {\tau }_{n}\sin {\tau }_{n'}}}+{\mathcal{O}}({({\tau }^{{\prime} })}^{3}).\end{array}$$

(15)

When it is not 0, the norm-2 of δV^(m) is bounded by 2. Hence $| | {\mathbb{E}}[\delta {V}^{(m)}]| | \le 2q({\tau }^{{\prime} })$. It follows that we have

$${\mathbb{E}}[{V}^{(1)}\ldots {V}^{(M)}]={\mathbb{E}}[{\tilde{V}}^{(1)}\ldots {\tilde{V}}^{(M)}]+\delta U,$$

(16)

with

$$\begin{array}{ll}| | \delta U| | \,\le \,\mathop{\sum }\limits_{m=1}^{M}\left(\begin{array}{c}M\\ m\end{array}\right){(2q)}^{m}={(1+2q)}^{M}-1\\ \qquad\quad=\,2t{\tau }^{{\prime} }{\mathop{\sum }\limits_{n < n'}\frac{{c}_{n}c_{n'}}{\sin {\tau }_{n}\sin {\tau }_{n'}}}+{\mathcal{O}}({({\tau }^{{\prime} })}^{2}).\end{array}$$

(17)

Hence the term δU vanishes when taking the Trotter limit ${\tau }^{{\prime} }\to 0$. The sequence of $M=t/{\tau }^{{\prime} }$ unitaries ${\tilde{V}}^{(m)}$ are all independent and are each of them different from I with probability

$$\frac{{\sum }_{n = 1}\frac{{p}_{n}}{1-{p}_{n}}}{1+{\sum }_{n = 1}\frac{{p}_{n}}{1-{p}_{n}}}={\tau }^{{\prime} }\mathop{\sum }\limits_{n=1}^{N}\frac{{c}_{n}}{\sin {\tau }_{n}}+{\mathcal{O}}({({\tau }^{{\prime} })}^{2}).$$

(18)

In probability theory, a process where events occur randomly and independently with probability ηδt in each time window δt converges by definition to a so-called Poisson process with rate η when δt → 0. Hence, in the limit ${\tau }^{{\prime} }\to 0$, the occurrence of unitaries ${\tilde{V}}^{(m)}$ that are different from I exactly converges to a Poisson process with rate $\mathop{\sum }\nolimits_{n = 1}^{N}\frac{{c}_{n}}{\sin {\tau }_{n}}$, applied during a time t. Since when ${\tilde{V}}^{(m)}\,\ne\, I$, the conditional probability that ${\tilde{V}}^{(m)}={e}^{i{\tau }_{n}{O}_{n}}$ is

$$\frac{{p}_{n}/(1-{p}_{n})}{\mathop{\sum }\nolimits_{i = 1}^{N}{p}_{i}/(1-{p}_{i})}=\frac{{c}_{n}/\sin {\tau }_{n}}{\mathop{\sum }\nolimits_{i = 1}^{N}{c}_{i}/\sin {\tau }_{i}}+{\mathcal{O}}({\tau }^{{\prime} }),$$

(19)

we obtain that in the limit ${\tau }^{{\prime} }\to 0$, the application of gate ${e}^{i{\tau }_{n}{O}_{n}}$ exactly follows a Poisson process with rate $\frac{{c}_{n}}{\sin {\tau }_{n}}$ during a time t. Since there is at most one gate per Trotter step in ${\tilde{V}}^{(m)}$, the rotations ${e}^{i{\tau }_{n}{O}_{n}}$ drawn from the Poisson process for different n’s are ordered only according to their time of occurrence, not according to n.

Namely, for each rotation ${e}^{i{\tau }_{n}{O}_{n}}$ we obtain a sequence of gate times $0 < {t}_{n}^{(1)} < \ldots < {t}_{n}^{({m}_{n})}\, <\, t$ drawn from a Poisson process with rate ${c}_{n}/\sin {\tau }_{n}$. From the Poisson process, this number of times m_n follows a Poisson distribution with parameter $\frac{t{c}_{n}}{\sin {\tau }_{n}}$, and once m_n has been drawn randomly, the random times ${t}_{n}^{(j)}$ are real numbers drawn uniformly and independently between 0 and t. For a given realization T of all these times for different n’s, we denote $\vert {\psi }_{T}\rangle$ the wave function obtained by applying the rotations ${e}^{i{\tau }_{n}{O}_{n}}$ ordered by time of occurrence ${t}_{n}^{(i)}$, irrespective of n. Such a configuration $\vert {\psi }_{T}\rangle$ is depicted schematically in Fig. 1. Then, denoting ${\mathbb{E}}$ the expectation value with respect to two independent configurations $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 }^{2}}.$$

(20)

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 discretization error, by only applying rotations with application time τ_n > 0 fixed. The number of rotations N_rot in each circuit is a random variable. Its average is

$${\mathbb{E}}[{N}_{{\rm{rot}}}]=2t\mathop{\sum }\limits_{n=1}^{N}\frac{{c}_{n}}{\sin {\tau }_{n}},$$

(21)

which is bounded independently of the precision wanted on the result, and only controlled by τ_n. We will study more in details the statistical distribution of this random variable in a section below.

Fig. 1: Example of a random circuit generated by the algorithm.

Each box represents a unitary operator U_O = e^iτO as appearing in the time evolution of a transverse-field Ising Hamiltonian, with time represented vertically and qubits horizontally. 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.

Statement of the algorithm

We state here the algorithm deduced from the previous explanations as a list of instructions that can be implemented by the reader. It takes as input a Hamiltonian H decomposed as in (7), a unitary observable ${\mathcal{M}}$, a time t and an initial state $\vert \psi (0)\rangle$, and outputs $\langle \psi (t)| {\mathcal{M}}| \psi (t)\rangle$. It takes 0 < τ_n ≤ π/2 as N parameters.

• 1. Prepare the system in the initial state $\vert \psi (0)\rangle$.

• 2. For each n ∈ {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 ∈ {1, …, m_n}, uniformly at random between 0 and t. These $M\equiv \mathop{\sum }\nolimits_{n = 1}^{N}{m}_{n}$ numbers are collected into a set T.

• 3. Deduce then the sequence k₁, …, k_M ∈ {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 rotation ${e}^{i{\tau }_{{k}_{m}}{O}_{{k}_{m}}}$ on the system.

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

• 5. Repeat steps (2) and (3) with τ_n replaced by − τ_n.

• 6. Measure the overlap with the initial state and divide the result by Λ² given in (13). On a quantum computer, this step requires performing a Hadamard test.

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

Measurement of the time evolution operator

The algorithm presented above can be adapted to computing the amplitude 〈ψ(0)∣e^iHt∣ψ(0)〉. It is often called Loschmidt amplitude in the condensed matter literature and can be used to determine ground state energies as well as microcanical expectation values^7,63. With the same notations as in Eq. (20), we have $\langle \psi (0)| {e}^{iHt}| \psi (0)\rangle ={\mathbb{E}}[\langle \psi (0)| {\psi }_{T}\rangle ]/\Lambda$. One thus skips steps (4) and (5) in the algorithm above, and replace the attenuation factor Λ² by Λ in step (6).

Sampling cost

Computing the expectation value $\langle \psi (t)| {\mathcal{M}}| \psi (t)\rangle$ through formula (20) involves a statistical average of quantum expectation values. When performing the measurement on a quantum computer, one has to choose a number of shots S to do per configuration, i.e. how to distribute the resources between number of configurations and accuracy of each configuration. If we neglect compilation time, the optimal number of shots can be seen to be S = 1 (see Supplementary Information). In that case, denoting o ∈ { + 1, − 1} the measurement outcome for every single-shot circuit, the expectation values and variances are

$$\begin{array}{ll}{\mathbb{E}}[o]\qquad\qquad\,=\,{\Lambda }^{2}\langle \psi (t)| {\mathcal{M}}| \psi (t)\rangle \\ {\mathbb{E}}[{o}^{2}]-{\mathbb{E}}{[o]}^{2}\,\,=\,1-{\Lambda }^{4}{|\langle \psi (t)| {\mathcal{M}}| \psi (t)\rangle |}^{2}\le 1.\end{array}$$

(22)

Hence, after dividing by Λ² the measured values, denoting $\bar{o}\in \{+{\Lambda }^{-2},-{\Lambda }^{-2}\}$ the random variable that is an unbiased estimator for $\langle \psi (t)| {\mathcal{M}}| \psi (t)\rangle$, we have

$$\begin{array}{ll}{\mathbb{E}}[\bar{o}]\qquad\qquad\,=\,\langle \psi (t)| {\mathcal{M}}| \psi (t)\rangle \\ {\mathbb{E}}[{\bar{o}}^{2}]-{\mathbb{E}}{[\bar{o}]}^{2}\,\,=\,{\Lambda }^{-4}-{|\langle \psi (t)| {\mathcal{M}}| \psi (t)\rangle| }^{2}\le {\Lambda }^{-4}.\end{array}$$

(23)

This means that to reach precision ϵ on the result, one has to perform around 1/(ϵ²Î›⁴) shots at most. These expressions can moreover be used to derive confidence intervals on the expectation value of the observable ${\mathcal{M}}$ using Hoeffding’s inequality. Our algorithm can thus provide guarantees that the expectation value of a time-evolved observable has a given high probability of being comprised in a window of some size ϵ.

Optimal angle on a perfect quantum computer

Formula (20) holds for any choice of gate angle τ_n. Increasing the gate angle τ_n decreases the number of rotations per circuit ${\mathbb{E}}[{N}_{{\rm{rot}}}]$, but also decreases the attenuation factor Λ², and so increases the number of shots Λ^−4 that have to be performed to reach a given precision. The total number of rotations $R(\{{\tau }_{n}\})={\mathbb{E}}[{N}_{{\rm{rot}}}]{\Lambda }^{-4}$ to perform to get a precision of order ${\mathcal{O}}(1)$ on the result is thus a non-trivial function of the gate angles τ_n, that reaches a minimum at some ${\tau }_{n}^{* }$ that can be computed numerically. At large t, we find analytically that the optimal angle is

$${\tau }_{n}^{* }=\frac{1}{2t\mu },$$

(24)

for all n = 1, …, N, where $\mu =\mathop{\sum }\nolimits_{n = 1}^{N}{c}_{n}$ is the sum of the coefficients of the Hamiltonian (see Supplementary Information). In that limit we find Λ⁴ = e^−1 which remains finite and ${\mathbb{E}}[{N}_{{\rm{rot}}}]=4{t}^{2}{\mu }^{2}$. Hence the average gate count per circuit scales as

$${\mathcal{O}}({t}^{2}{\mu }^{2}).$$

(25)

The number of shots to perform to reach precision ϵ scales as ${\mathcal{O}}({\epsilon }^{-2})$, and so the runtime (number of shots times number of gates per shot) to reach precision ϵ scales as ${\mathcal{O}}({t}^{2}{\mu }^{2}{\epsilon }^{-2})$. But remarkably, the precision ϵ does not enter the average gate count per circuit, which remains finite even when ϵ → 0. This differs from usual algorithms such as Trotter, qDRIFT or LCU, whose gate count per circuit we recall in Table 1.

Table 1 Algorithms’ gate counts

Optimal angle on a noisy quantum computer

The optimal angle value (24) holds for a perfect, noiseless quantum computer. We now consider the large-t optimal angle for a reasonable model of hardware noise, that consists in assuming that the application of each gate ${e}^{i{\tau }_{n}{O}_{n}}$ comes with a signal attenuation ${e}^{-{r}_{n}}$ with some r_n > 0. This kind of noise model is known to be a rough first approximation of the effect of the noise on current hardware. Moreover, other types of noise can be converted into this global depolarizing noise through some noise mitigation techniques^6,64,65. In state-of-the-art trapped-ion devices typical two-qubit gate fidelities are e^−r ~ 99.8%, in which case r ~ 2 ⋅ 10^{−3 66}. There are on average $\frac{t{c}_{n}}{\sin {\tau }_{n}}$ such rotations per configuration. The damping due to hardware imperfection is thus

$$q=\exp \left(-2t\mathop{\sum }\limits_{n=1}^{N}\frac{{r}_{n}{c}_{n}}{\sin {\tau }_{n}}\right).$$

(26)

The total attenuation of the signal is Λ²q. Assuming r_n small, the optimal times ${\tau }_{n}^{* }$ that minimize the total gate count per circuit at large t are now (see Supplementary Information)

$${\tau }_{n}^{* }=\sqrt{2{r}_{n}}.$$

(27)

The total attenuation Λ²q scales then as

$${\Lambda }^{2}q=\exp \left(-2t\mathop{\sum }\limits_{n=1}^{N}\sqrt{2{r}_{n}}{c}_{n}\right).$$

(28)

To reach precision ϵ on the result, our algorithm requires thus $1/{({\Lambda }^{2}q\epsilon )}^{2}$ shots. Importantly, the dependence in the precision is only polynomial. As a comparison, let us consider the same noise model for a Trotter implementation of the time evolution. For simplicity, we take r_n = r constant. In the Trotter algorithm, absent hardware error, a Trotter step of size τ has error N²tCτ for the first order Trotter formula, with some coefficient C. To reach precision ϵ, we thus necessarily have to take τ < ϵ/(N²tC), hence run a circuit with more than N³t²C/ϵ gates. For local Hamiltonians, this scaling improves to N²t²C/ϵ. Then, under the same noise model, the imperfections in the hardware incur a signal attenuation of order ${e}^{-{N}^{3}{t}^{2}Cr/\epsilon }$. The number of shots M_Trotter thus has to be at least

$${M}_{{\rm{Trotter}}} \sim \frac{{e}^{2{N}^{3}{t}^{2}Cr/\epsilon }}{{\epsilon }^{2}}.$$

(29)

In contrast with our algorithm the dependence in the precision is exponential. In particular, our algorithm requires exponentially fewer shots for a given precision ϵ when

$$\epsilon \,<\, \frac{t{N}^{3}C\sqrt{r}}{2\sqrt{2}\mathop{\sum }\limits_{n=1}^{N}{c}_{n}}.$$

(30)

For local Hamiltonians, the bound would be N² instead of N³. In both cases, we see that our algorithm always outperforms Trotter at high precision in presence of noise.

Distribution of the number of gates per circuit

We stated in Eq (25) that the average gate count in the random circuits generated by our algorithm is ${\mathcal{O}}({t}^{2}{\mu }^{2})$, without any dependence on the precision ϵ achieved on the result. In our algorithm, the number of rotations is a random variable

$${N}_{{\rm{rot}}}=\pi \left(\mathop{\sum }\limits_{n=1}^{N}\frac{2t{c}_{n}}{\sin {\tau }_{n}}\right),$$

(31)

with π(λ) a Poisson distribution with parameter λ. Indeed, the number of each rotation ${e}^{i{\tau }_{n}{O}_{n}}$ follows a Poisson distribution with parameter $\frac{t{c}_{n}}{\sin {\tau }_{n}}$, and the sum of Poisson random variables follows a Poisson distribution with parameter given by the sum of their parameters. The variance of this random variable is

$${\mathbb{E}}[{N}_{{\rm{rot}}}^{2}]-{\mathbb{E}}{[{N}_{{\rm{rot}}}]}^{2}=2t\mathop{\sum }\limits_{n=1}^{N}\frac{{c}_{n}}{\sin {\tau }_{n}}.$$

(32)

At the optimal angle on a perfect quantum computer (24), this is

$${\mathbb{E}}[{N}_{{\rm{rot}}}^{2}]-{\mathbb{E}}{[{N}_{{\rm{rot}}}]}^{2}=4{t}^{2}{\mu }^{2}+{\mathcal{O}}(1).$$

(33)

Hence the gate count per circuit is in average ${\mathcal{O}}({t}^{2}{\mu }^{2})$ with standard deviation ${\mathcal{O}}(t\mu )$.

A particularity of our algorithm is that, because the Poisson distribution is unbounded, the random number of gates in the random circuits is unbounded. Using Chernov’s bound, we have that the probability that the number of rotations N_rot is k + 1 times the average ${\mathbb{E}}[{N}_{{\rm{rot}}}]$ is⁶⁷

$${\mathbb{P}}[{N}_{{\rm{rot}}}\ge (k+1){\mathbb{E}}[{N}_{{\rm{rot}}}]]\le 2\exp \left(-\frac{{k}^{2}}{2(k+1)}{\mathbb{E}}[{N}_{{\rm{rot}}}]\right).$$

(34)

The number of gates is thus strongly concentrated around the average value. Let us assume that one imposes a maximal number of rotations that can be run on the quantum computer, namely one discards all circuits with more than $(k+1){\mathbb{E}}[{N}_{{\rm{rot}}}]$ rotations for a given k. We note first that for all practical purposes, taking say $k=7/\sqrt{{\mathbb{E}}[{N}_{{\rm{rot}}}]}=3.5/(t\mu )$ at optimal gate angle and large tμ ensures a negligible discard probability < 10^−10. More theoretically, one will have to draw ${\mathcal{O}}({e}^{{k}^{2}{\mathbb{E}}[{N}_{{\rm{rot}}}]/(2(k+1))})$ circuits before discarding one circuit and biasing the result. Before that, one achieves a precision ${\mathcal{O}}({e}^{-{k}^{2}{\mathbb{E}}[{N}_{{\rm{rot}}}]/(4(k+1))})$ on the expectation value measured. For this to be smaller than some required precision η, with ${\mathbb{E}}[{N}_{{\rm{rot}}}]=4{t}^{2}{\mu }^{2}$ at the optimal gate angle, one can take $k={\mathcal{O}}(\sqrt{\log (1/\eta )}/(t\mu ))$ at large tμ and small η, in the regime $\sqrt{\log (1/\eta )}/(t\mu )\to 0$. It follows that the maximal number of rotations that one needs to include to guarantee precision η is

$${\mathcal{O}}({t}^{2}{\mu }^{2})+{\mathcal{O}}(t\mu \sqrt{\log 1/\eta }).$$

(35)

We thus see that the maximal gate count per circuit of our algorithm possesses a dependence in the precision. However, the average gate count per circuit is independent of the precision, contrary to algorithms like LCU where the gate count per circuit is deterministic with a dependence in the precision, for all circuits. This translates into the fact that the runtime of our algorithm, i.e. number of shots times number of gates per shot, scales as

$${\mathcal{O}}\left(\frac{{t}^{2}{\mu }^{2}}{{\epsilon }^{2}}\right)$$

(36)

to reach precision ϵ on the expectation value measured. This is in contrast with previously known algorithms whose runtime would be the gate count per circuit in Table 1 times ϵ^−2. Our algorithm has thus a better performance in runtime with respect to the precision ϵ.

Time-dependent Hamiltonian

The algorithm can be generalized straightforwardly to dynamics with a time-dependent Hamiltonian H(t). Let us decompose

$$H(t)=\mathop{\sum }\limits_{n=1}^{N}{c}_{n}(t){O}_{n},$$

(37)

with time-dependent coefficients c_n(t). For simplicity and without loss of generality, we will assume c_n(t) > 0 for all times. We denote $\vert \psi (t)\rangle ={\mathcal{T}}\exp (i\mathop{\int}\nolimits_{0}^{t}{\rm{d}}sH(s))\vert \psi (0)\rangle$, namely the state $\vert \psi (t)\rangle$ that satisfies the differential equation

$${\partial }_{t}\left\vert \psi (t)\right\rangle =iH(t)\left\vert \psi (t)\right\rangle .$$

(38)

As in the time-independent case, this time-ordered exponential can be obtained as the limit of a product formula

$${\mathcal{T}}\exp \left(i\mathop{\int}\nolimits_{0}^{t}{\rm{d}}sH(s)\right)=\mathop{\lim }\limits_{{\tau }^{{\prime} }\to {0}^{+}}\mathop{\prod }\limits_{j=1}^{t/{\tau }^{{\prime} }}\left({e}^{i{\tau }^{{\prime} }{c}_{1}(j{\tau }^{{\prime} }){O}_{1}}\ldots {e}^{i{\tau }^{{\prime} }{c}_{N}(j{\tau }^{{\prime} }){O}_{N}}\right).$$

(39)

For each n = 1, …, N, we choose 0 < τ_n ≤ π/2 and implement ${e}^{i{\tau }^{{\prime} }{c}_{n}(j{\tau }^{{\prime} }){O}_{n}}$ using (1) with angle τ_n. In the sequence of rotations appearing on the right-hand side of (39), instead of each rotation ${e}^{i{\tau }^{{\prime} }{c}_{n}(s){O}_{n}}$ we thus apply ${e}^{i{\tau }_{n}{O}_{n}}$ with probability p_n given by ${p}_{n}=\frac{{\tau }^{{\prime} }{c}_{n}(s)}{\sin {\tau }_{n}}+{\mathcal{O}}({({\tau }^{{\prime} })}^{2})$, that now depends on time s. As in the time-independent case, this probability is proportional to ${\tau }^{{\prime} }$ as ${\tau }^{{\prime} }\to 0$, and this random choice is tried ${\mathcal{O}}(1/{\tau }^{{\prime} })$ times. The limit ${\tau }^{{\prime} }\to 0$ thus exactly defines a Poisson process but with a time-dependent rate given by $\frac{{c}_{n}(s)}{\sin {\tau }_{n}}$. This Poisson process with a time-dependent rate can be mapped to an ordinary Poisson process after a mapping of the time. Namely, we introduce the function

$${z}_{n}(u)={\int \nolimits_{0}^{u}}{c}_{n}(s){\rm{d}}s,$$

(40)

as well as ${z}_{n}^{-1}(u)$ its reciprocal, i.e. such that ${z}_{n}^{-1}({z}_{n}(u))=u$. The random times generated by the Poisson process with time-dependent rate $\frac{{c}_{n}(s)}{\sin {\tau }_{n}}$ are obtained by applying the function ${z}_{n}^{-1}$ to the times generated by a Poisson process with constant rate $\frac{1}{\sin {\tau }_{n}}$ during a time z_n(t).

As for the product of the attenuation factors λ_n’s, it becomes

$$\mathop{\lim }\limits_{{\tau }^{{\prime} }\to 0}\prod _{m}{\lambda }_{m}=\exp \left(-\mathop{\sum }\limits_{n=1}^{N}\tan ({\tau }_{n}/2)\mathop{\int}\nolimits_{0}^{t}{c}_{n}(s){\rm{d}}s\right)\equiv \Lambda .$$

(41)

To summarize, the algorithm for the time-independent case is modified as follows. We replace step (2) of the time-independent case by

• 2. For each n ∈ {1, …, N}, draw an integer m_n from a Poisson distribution with parameter $\frac{1}{\sin {\tau }_{n}}$. Then draw m_n real numbers ${\tilde{t}}_{n}^{(i)}$, where i ∈ {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 \mathop{\sum }\nolimits_{n = 1}^{N}{m}_{n}$ numbers ${t}_{n}^{(i)}$ are collected into a set T.

Then in step (3), the rotation applied is ${e}^{i{\tau }_{{k}_{m}}{O}_{{k}_{m}}}$ where $s\in \{{t}_{n}^{(i)}\}$ is the time at which the rotation is applied. For the backward propagation in step (5) one has to reverse the order of the gates (this was not necessary in the time-independent case as the Poisson process was invariant by time reversal). The total attenuation Λ of step (6) is replaced by (41).

This algorithm produces the exact expectation value of an observable under the time-dependent Hamiltonian H(t), without errors arising from discretizing the continuous function c_n(t) or Trotter errors.

Background Hamiltonian

In our algorithm, one can choose an arbitrary different gate angle τ_n for every rotation ${e}^{i{\tau }_{n}{O}_{n}}$. This possibility of choosing different gate angles can be used to improve the algorithm the following way. 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}.$$

(42)

The gate angle τ_p of the rotations corresponding to these terms involved in (20) can be chosen freely. Let us investigate what happens if we take the limit τ_p → 0 for ${O}_{p}\in {\mathcal{O}}$. We note that this limit is not the Trotter limit ${\tau }^{{\prime} }\to 0$ appearing in (8). Here, the limit ${\tau }^{{\prime} }\to 0$ has been taken first at fixed τ_p, and we then take the limit τ_p → 0. In this limit, the rate of the Poisson process corresponding to these gates $\frac{{c}_{p}}{\sin {\tau }_{p}}$ diverges to ∞, while the rates of the Poisson process of the other gates $\frac{{c}_{n}}{\sin {\tau }_{n}}$ not in ${\mathcal{O}}$ are not affected. Hence, when drawing the times at which the rotations are applied, almost all of them will correspond to gates in ${\mathcal{O}}$. If we take two times t₁ < t₂ at which rotations of terms not in ${\mathcal{O}}$ are applied, there will be typically ${\mathcal{O}}(\frac{{c}_{p}}{{\tau }_{p}})$ rotations ${e}^{i{\tau }_{p}{O}_{p}}$ applied in between when τ_p → 0. Hence, since all these terms commute among themselves, when we take τ_p → 0, these rotations will combine into ${e}^{i({t}_{2}-{t}_{1}){H}_{{\rm{background}}}}$.

Therefore the algorithm can be modified as follows. In step (2) we only draw times corresponding to terms not in ${\mathcal{O}}$, and use them to produce an ordered sequence of gates as in step (3). However, between each application of these rotations ${e}^{i{\tau }_{n}{O}_{n}}$ and ${e}^{i{\tau }_{m}{O}_{m}}$ at times s_n < s_m, we apply ${e}^{i\delta s{H}_{{\rm{background}}}}$ with δs = s_m − s_n the difference between the two instants at which the rotations are applied. The Hamiltonian H_background thus plays the role of a “background" Hamiltonian that is always applied on the system, but on top of which rotations ${e}^{i{\tau }_{n}{O}_{n}}$ with ${O}_{n}\,\notin\, {\mathcal{O}}$ are applied.

The upshot of this improvement is that it reduces the attenuation factor Λ, which is now computed only including the terms not in ${\mathcal{O}}$

$$\Lambda =\exp \left(-t\sum _{\begin{array}{c}n=1,\ldots ,N\\ {O}_{n}\notin {\mathcal{O}}\end{array}}{c}_{n}\tan ({\tau }_{n}/2)\right).$$

(43)

Comparison with previous work

After completion and online publication of this work, the recent paper⁶⁸ was brought to our attention, where a simulation algorithm called URCC was introduced. It uses a completely different decomposition and sampling strategy, based on Dyson expansion and leading order rotations, and is distinct from ours. However, it has a gate count per circuit ${\mathcal{O}}({t}^{2}{\mu }^{2})$ independent of the precision like this present work.

Let us perform a comparison of runtime and gate count between this algorithm and ours⁶⁸. The URCC algorithm takes as a parameter a number of segment N_seg. The variance of the measured expectation value after M shots is

$${V}_{{\rm{URCC}}}=\frac{f{\left(\frac{t\mu }{{N}_{{\rm{seg}}}}\right)}^{4{N}_{{\rm{seg}}}}}{M},$$

(44)

with a function f(x) that behaves at small x as $f(x)=1+{x}^{2}+{\mathcal{O}}({x}^{3})$. To ensure that the variance remains of order ${\mathcal{O}}(1)$ at large tμ, one needs to choose ${N}_{{\rm{seg}}}={\mathcal{O}}({t}^{2}{\mu }^{2})$. The random circuits generated contain always ${{\mathcal{N}}}_{{\rm{URCC}}}=2{N}_{{\rm{seg}}}$ rotations (possibly of angle π/2). In our case, taking τ_n = τ independent of n, the variance of the measured expectation value is

$${V}_{{\rm{ours}}}=\frac{{e}^{4t\mu \tan (\tau /2)}}{M}.$$

(45)

The number of rotations in our case is a random variable that is unbounded, whose average is ${{\mathcal{N}}}_{{\rm{ours}}}=\frac{2t\mu }{\sin \tau }$. To ensure that the variance remains of order ${\mathcal{O}}(1)$, one needs to choose $\tau ={\mathcal{O}}({(t\mu )}^{-1})$. Choosing N, τ so as to equate the variances V_ours = V_URCC at large tμ gives τ = 2tμ/N_seg. Then one obtains an average number of gates per circuit ${{\mathcal{N}}}_{{\rm{ours}}}={{\mathcal{N}}}_{{\rm{URCC}}}/2$. Although the gate count is a random variable that is unbounded in our case, the average gate count is twice as small in our algorithm, at same overhead in terms of shots. The interpretation of this difference is that our algorithm allows for optimizing the gate angles so as to minimize the average gate count. Although the decomposition in terms of rotations and Pauli strings (which are rotations with angle π/2) used in the URCC algorithm cleverly guarantees bounded gate count, it is suboptimal in terms of average gate count.

Numerical tests

We now show numerical simulations of our algorithm. We start with a noisy implementation on a non-sparse Hamiltonian as routinely encountered in quantum chemistry. Hamiltonians describing molecular electronic structure problems are typically decomposed in a basis set of fermionic orbitals, where they read

$$H=\mathop{\sum }\limits_{ij=1}^{N}{h}_{ij}{c}_{i}^{\dagger }{c}_{j}+\mathop{\sum }\limits_{ijkl=1}^{N}{h}_{ijkl}{c}_{i}^{\dagger }{c}_{j}^{\dagger }{c}_{k}{c}_{l},$$

(46)

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 per circuit for any Trotter implementation of the dynamics on near-term hardware. However, the sum of the coefficients of the Hamiltonian when written in terms of Pauli strings has a lower scaling with N, which makes our algorithm attractive^69,70,71. We consider a H₂O water molecule in a STO-6G basis with 14 orbitals, and use our algorithm to study the Loschmidt amplitude

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

(47)

where $\vert 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 Å, to have a significant difference between the exact ground state and $\vert HF\rangle$. On a noisy hardware, we propose to compute the ratio

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

(48)

with ℜ and ℑ denoting the real and imaginary parts. This ratio is not sensitive to the depolarizing part of the hardware noise, but still contains information about the energy levels in the oscillation frequencies. In panel (a) of Fig. 2 we show a noisy simulation of our algorithm that includes a depolarizing channel with probability 0.002 after each 2-qubit gate, and observe very good agreement.

Fig. 2: Examples of numerical implementations.

a: chemistry Hamiltonian (see text), with gate angles (27), with 2000 circuits with 20 shots each, noiseless (teal) and noisy (orange) simulations, and exact value (black). b: adiabatic state preparation (see text), noiseless simulations with 10⁵ circuits and τ = 0.04 (teal, the shade indicating one standard deviation), and exact value (black).

Next, we show numerical results for the time-dependent case of our algorithm. For this case, we consider the 2D square lattice Ising model in a transverse field h, whose Hamiltonian is

$$H=-\sum _{\langle i,j\rangle }{Z}_{i}{Z}_{j}-h\sum _{j}{X}_{j},$$

(49)

with periodic boundary conditions. We implement the adiabatic preparation of the ground state at h_f = 2.5, starting from h = 0. To connect these magnetic field values adiabatically we use the schedule $h(t)={h}_{{\rm{f}}}\sin {(\frac{\pi }{2}\sin {(\frac{\pi t}{2{T}_{{\rm{f}}}})}^{2})}^{2}$ for 0 ≤ t ≤ T_f. In panel (b) of Fig. 2 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 τ, as imposed by (13).

Finally, we compare our algorithm to first-order Trotter and qDRIFT in the case of the 2D square lattice Ising model (49) with magnetic field h = 2, which is a typical model studied in condensed matter. In Fig. 3 we plot the Loschmidt amplitude 〈0∣e^iHt∣0〉 starting from the Z = + 1 product state $\vert 0\rangle$, using a same gate angle τ = 0.1 for our algorithm and for qDRIFT, and with a same time step 0.1 for Trotter. We take 10⁵ shots per point (which is necessary to have a precision < 0.01 with high confidence), and distribute them over 10³ random circuits with 10² shots per circuit for our case and for qDRIFT. On the time points in Fig. 3, the average difference in absolute value between exact expectation value and simulated is 0.064 for Trotter, 0.135 for qDRIFT and 0.009 for our method, at equal runtime. We see that our algorithm performs better than both qDRIFT and Trotter. This demonstrates our algorithm can perform well in problems of scientific relevance.

Fig. 3: Comparison of performance with previous algorithms.

Real part of the Loschmidt amplitude 〈0∣e^itH∣0〉 in the 2D Ising model at h = 2 in size 3 × 4 with $\vert 0\rangle$ the Z = + 1 product state, comparing Trotter, qDRIFT and our algorithm, all with the same gate angle and Trotter step equal to τ = 0.1, and all with 10⁵ shots per point. For qDRIFT and our algorithm, the 10⁵ shots are distributed over 10³ random circuits with 10² shots per circuit.

Between quantum evolution and path integral representation

In formula (20), the gate angle τ of the rotations is a free parameter that does not modify the precision of the result, but that modifies the attenuation factor Λ. Let us comment on two interesting limits of this gate angle τ, that are τ → 0 and τ = π/2. When τ → 0, the different rotations ${e}^{i\tau {O}_{n}}$ commute at order τ and there are increasingly more rotations drawn from the Poisson process. With probability 1 in this limit, a single configuration $\vert {\psi }_{T}\rangle$ is equal to the exact time evolution ${e}^{iHt}\vert 0\rangle$, and the attenuation factor Λ is equal to 1. This corresponds to a “purely quantum" limit of (20), where no averaging over circuits is required and where the exact quantum state is prepared.

On the other hand, when τ = π/2, we have ${e}^{i\tau {O}_{n}}=i{O}_{n}$, and so when O_n is a string of Pauli operators, the rotations can be implemented with only one-qubit gates. Each configuration $\vert {\psi }_{T}\rangle$ remains thus a product state in the Z basis and describes a classical trajectory in discrete time. It is thus a “purely classical" limit. For that value of τ, formula (20) becomes thus 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 thus an interesting interpretation of our algorithm as an interpolation between a single quantum trajectory for τ = 0, and exponentially many classical trajectories for τ = π/2. Intermediate values of τ near π/2 allow one to decrease the entanglement present in each configuration $\vert {\psi }_{T}\rangle$ compared to the exact time-evolved state ${e}^{iHt}\vert 0\rangle$, but at the cost of increasing the number of samples required through the attenuation factor (13). This could find applications in the classical simulation of quantum dynamics.

Simulation of Clifford circuits + T gates

It is well-known that Clifford circuits (i.e. circuits built out of S, H, CNOT gates) can be simulated in polynomial time⁷². These circuits are however not universal, in the sense that it is not possible to approximate arbitrary unitary operators with arbitrary precision with Clifford circuits⁷³. To obtain a universal gate set, one has to include as well T gates, defined as $T=\left(\begin{array}{cc}1&0\\ 0&{e}^{i\pi /4}\end{array}\right)$. It is then natural to ask the complexity of simulating circuits that contain Clifford gates and t-many T gates.

The basic mechanism of gate angle reduction in (1) can be applied as well here. Let us consider a unitary U that we decompose as

$$U={U}_{1}T{U}_{2}T{U}_{3}\ldots {U}_{t}T{U}_{t+1},$$

(50)

where the T gates are applied on some qubits (whose index we omit for readability), and each unitary U_j can be implemented with only Clifford gates. We would like to evaluate the runtime for computing 〈0∣U∣0〉 classically. Writing

$$T=\frac{{e}^{i\pi /8}}{2\cos (\pi /8)}(I+{e}^{-i\pi /4}S),$$

(51)

we can express

$$\begin{array}{ll}\langle 0| U| 0\rangle \,=\,\frac{{e}^{it\pi /8}}{\cos {(\pi /8)}^{t}}\\ \qquad\quad\,\times \,\frac{1}{{2}^{t}}\mathop{\sum}\limits _{{V}_{1},\ldots ,{V}_{t}\in \{I,S\}}{e}^{-in\pi /4}\langle 0| {U}_{1}{V}_{1}{U}_{2}\ldots {U}_{t}{V}_{t}{U}_{t+1}| 0\rangle ,\end{array}$$

(52)

with n the number of times S is chosen in the sum. Each circuit U₁V₁U₂…U_tV_tU_t+1 is a Clifford circuit and so can be simulated in polynomial time. To compute exactly the sum, one has to compute 2^t such Clifford circuits. However, one can also approximate probabilistically the sum with the 1/2^t prefactor by drawing each term I or S with probability 1/2. To obtain precision ϵ on the amplitude 〈0∣U∣0〉, because of the factor $1/\cos {(\pi /8)}^{t}$, one will have to include of order $1/({\epsilon }^{2}\cos {(\pi /8)}^{2t})$ terms. Hence the cost in simulating the Clifford circuit with t-many T gates up to precision ϵ is exponential in t, with exponent ${e}^{2| \log \cos (\pi /8)| t}$. The exponent turns out to match the currently best known exponent⁷⁴.

DISCUSSION

We presented a quantum algorithm to compute expectation values of time-evolved observables that requires an average circuit depth that is independent of the desired precision ϵ, contrary to many algorithms that require an infinite depth when ϵ → 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 discretization error while having bounded-depth circuits in average. This amplification factor increases the number of shots required to reach a given precision. The algorithm is particularly adapted to non-sparse Hamiltonians as the gate count per circuit only depends on the sum of the coefficients of the Hamilotnian and not the number of terms. The gate count per circuit for simulation time t is ${\mathcal{O}}({t}^{2}{\mu }^{2})$, and with shot cost ${\mathcal{O}}({\epsilon }^{-2})$ to reach precision ϵ. We also showed that the algorithm generalizes to time-dependent Hamiltonian, achieving again no time-discretization error or no Trotter error with a finite gate count per circuit. We finally showed that the same mechanism can be used to simulate classically Clifford circuits with T gates, with a runtime as good as the best known simulation algorithms.

We expect this algorithm for implementing time evolution on a quantum computer to have particular impact on certain applications. For example, the absence of discretization error is particularly adapted to quantum chemistry applications where extreme precision is required. It is also particularly relevant for the quantum adiabatic algorithm due to the absence of time discretization errors. One shortcoming of the approach is that it is not adapted to perform sampling (instead of computing expectation values), as occurs for example in using quantum computers to solve classical optimization problems. It would be valuable to see if an algorithm with similar characteristics could be devised for performing sampling. We will study more in depth these applications in further works.