Phase estimation of local Hamiltonians on NISQ hardware

Here we recap our theoretical abstraction of the hardware capabilities of NISQ-era devices [CBC21]. We describe three related concepts in this section: a circuit model, an error model, and a cost model. These concepts are summarized by definitions 1 and 2 and equation (1). The circuit model can describe any quantum computing architecture in which discrete gates are generated by continuous-time interactions between qubits, and where these interactions can effectively be described by a Hamiltonian—or effective Hamiltonian—coupling those qubits.

Many current matter-based quantum computing architectures fit this description, including superconducting circuits, most ion trap architectures and nuclear magnetic resonance (NMR) quantum computing (liquid-state or solid-state). Although this error model is a theoretical idealisation work has recently been published which supports the validity of this model and demonstrates the advantages of exploiting analogue control [Ste+21].

We call this circuit model the subcircuit model and it is defined as follows.

Definition 1  (Sub-circuit model).Given a set of qubits Q, a set $I \subseteq Q \times Q$ specifying which pairs of qubits may interact, a fixed two qubit interaction Hamiltonian h, and a minimum switching time t_min, a sub-circuit pulse-sequence C is a quantum circuit of L pairs of alternating layer types $C = \prod_l^L U_l V_l$ with $U_l = \prod_{i \in Q} u_i^l$ being a layer of arbitrary single qubit unitary gates, and $V_l = \prod_{ij \in \Gamma_l} v_{ij} (t_{ij}^l )$ being a layer of non-overlapping, variable time, two-qubit unitary gates:

$$\begin{align*} v_{ij}(t)=\mathrm{e}^{\mathrm{i} t h_{ij}}\end{align*}$$

with the set $\Gamma_l \subseteq I$ containing no overlapping pairs of qubits, and $t \geqslant t_\mathrm{min}$. Throughout this paper we assume $h_{ij} = Z_iZ_j$. As all $\sigma_i \sigma_j$ are equivalent to $Z_i Z_j$ up to single qubit rotations this can be left implicit and so we take $h_{ij} = \sigma_i \sigma_j$.

This defines the gate-set we have access to. We also refer to these gates as 'sub-circuit' as they exist one level of abstraction below the digital gate-set of the standard circuit model. We note that any standard quantum circuit can be expressed in this 'sub-circuit' form as up to single qubit rotations a controlled NOT or CNOT gate is equivalent to $\mathrm{e}^{-{i} \frac{\pi}{4}ZZ}$. This means that we can cost standard circuits in our framework as sub-circuit gates with $t = \mathrm{O}(1)$. We will use this fact in section 4.3. We will express all quantum circuits in terms of sub-circuit gates throughout this work. We do this so that we can cost all circuits under the same noise model.

The sub-circuit model is useful as it can distinguish circuits under an error model where errors accumulate proportionally to the time for which interactions are switched on. This is in contrast to traditional error models, where errors accumulate according to gate depth and all gates are costed equally. This leads us to the following error model. We will assume that errors accumulate proportionally to the time interactions in the system are switched on for. Formally, we assume that for Q qubits every layer l of sub-circuit gates in a circuit $C = \prod_{l}U_l V_l$ is interspersed with a single qubit depolarizing noise channel $\mathcal{E}$. For each qubit in layer l the channel is then approximated as

$$\begin{align} \mathcal{E}_l(\rho): \rho \rightarrow (1-q|t^l|) \rho + q |t^l| \frac{I}{2}, \end{align} \tag{ 1 }$$

for some noise parameter $q \in [0,1]$ and where we have upper bounded each individual pulse-time in the lth-layer as $|t^l|: = \text{max}_{ij \in \Gamma_l}\left(|t^l_{ij}|\right)$.

Finally, we define the following cost model. We have assumed that errors accumulate based on the length of time two-qubit interactions are switched on for and therefore we cost circuits according to the following definition.

Definition 2  (Sub-circuit cost).The physical run-time of a sub-circuit pulse-sequence C is defined as

$$\begin{align*} \mathcal{T}(C) := \sum_l^L \max_{ij \in \Gamma_l}\left(t_{ij}^l\right). \end{align*}$$

The run-time is normalised to the physical interaction strength, so that $\vert h\vert = 1$. This cost metric assumes that the single qubit layers contribute a negligible amount to the total time duration of the circuit.

This cost-metric will help us estimate what circuits can and cannot be implemented for a given depolarizing noise rate q if we are willing to tolerate a target error $\epsilon_{\text{tar}}$.

This is captured by the following. We can calculate a trivial error bound by considering the probability that no error occurs at all. Let $\mathcal{U}(\rho)$ be the channel which implements circuit C without noise, and $\mathcal{U^{^{\prime}}}(\rho)$ the channel which implements circuit C with noise. Then for any measured observable $\mathrm{O}$ we get

$$\begin{align} \epsilon := \mathrm{Tr} \left[\left(\mathcal{U}(\rho)-\mathcal{U^{^{\prime}}}(\rho)\right)\mathrm{O}\right] \leqslant \epsilon_{\text{tar}}, \end{align} \tag{ 2 }$$

as we want to determine the maximum allowable run-time which keeps this stochastic error ε below a target error $\epsilon_{\text{tar}}$.

This stochastic error ε is upper bounded by the probability that any error occurs in the circuit. For a Q-qubit circuit the noise channel defined above this gives us

$$\begin{align} \epsilon &\leqslant 1 - \prod_{l}(1-q |t^l|)^Q \end{align} \tag{ 3 }$$

$$\begin{align} & \phantom{0}\kern1.9pt \leqslant 1 - (1-q)^{V} \leqslant \epsilon_{\text{tar}} , \end{align} \tag{ 4 }$$

where $V = Q \times \mathcal{T}(C)$ is the volume of the circuit and $\mathcal{T}(C) = \sum_{l} |t^l|$. Inverting this, we say that given an error parameter $q \in [0,1]$, we can only implement circuits whose run-time is bounded as

$$\begin{align} \mathcal{T}(C) \leqslant \frac{\log(1-\epsilon_{\text{tar}})}{Q \log(1-q)}, \end{align} \tag{ 5 }$$

assuming we are willing to tolerate a stochastic error rate at most $\epsilon_\mathrm{tar}$. We will use this in section 4.3 and particularly in figure 1. These are the key assumptions we make in our circuit error model. In summary, we assume there is sufficient analogue control to express circuits in the sub-circuit form of definition 1 and that depolarizing noise in the device can be described by equation (1). This in turn makes the cost-metric of definition 2 useful as then the trivial stochastic error bound equation (5) can help us estimate requisite depolarizing error rates q for a given task. Our analysis and benchmarking in section 4 will follow this logic.

Zoom In Zoom Out Reset image size

Having established this circuit error model in our previous work [CBC21], we also derived novel exact gate synthesis techniques to minimize the sub-circuit cost $\mathcal{T}$ in this error model. These results can be summarized as follows. For a k-local Pauli interaction P_k , there exists a sub-circuit pulse-sequence $C: = \prod_l^{L} U_l V_l$ which implements the evolution operator $\mathrm{e}^{-\mathrm{i} t P_k}$. Most importantly, for any target time t the run-time of that circuit is bounded as

$$\begin{align} \mathcal{T} \left(C \right) \leqslant \mathrm{O}\left( |t|^{\frac{1}{k-1}}\right), \end{align} \tag{ 6 }$$

according to the notion of run-time established in definition 2. Note that this is the scaling in terms of $|t|$, with k assumed to be a constant. We see from this that we can only expect an advantage for smaller values of k. The details of the decompositions which accomplish this can be found in our previous work [CBC21, appendix B]. Here we only note that they are exact and minimize $\mathcal{T}(C)$ as $t \rightarrow 0$ at the expense of increased circuit depth.

Let f(t) be an integrable function. For consistency with past literature we adopt the convention for the Fourier transform and its inverse to be

$$\begin{align} F(w) = \mathcal F[{f\,}](w) = \frac{1}{2\pi} \int_{\mathbb{R}} f(t) \mathrm{e}^{-\mathrm{i} w t} \mathrm{d} t, \quad\text{and}\quad \mathcal F^{-1}[F](t) = \int_{\mathbb{R}} F(w) \mathrm{e}^{\mathrm{i} w t} \mathrm{d} w. \end{align} \tag{ 7 }$$

For a periodic function f one can express f(t) as a Fourier series via

$$\begin{align} f(t) = \sum_{w=-\infty}^\infty F_w \mathrm{e}^{\mathrm{i} w t} \quad\text{with coefficients}\quad F_w := F(w). \end{align} \tag{ 8 }$$

With these conventions, the Fourier inversion theorem reads $\mathcal F^{-1} = \mathcal F\mathcal R$ for the reflection functional $R[\,f\,](t) : = f(-t)$.

QPE describes a family of quantum algorithms for retrieving information about the phases φ_j associated with the eigenvalues $\mathrm{e}^{\mathrm{i} \phi_j}$ of a unitary operator U. Such a unitary operator is generated by a Hermitian operator H via $U = \exp(-i H) = \sum_{n} \exp(-i \lambda_n) \vert \phi_n\rangle\langle \phi_n\vert$, where $\{\lambda_n \}$ are the eigenvalues of H—which are equivalent to the phases φ_j of U modulo 2π—and where the $\{\vert \phi_n\rangle \}$ are the (same) eigenvectors of H and U, respectively.

All variants of phase estimation may be thought of in the following terms. Consider some input state ρ and a time evolution operator U(t) with phases $\{\phi_j\}$ and eigenstates $\{\vert \phi_j\rangle\}$. Prepare some time series data

$$\begin{equation} g(t) := \sum_n \mathrm{e}^{\mathrm{i} \phi_n t} \langle \phi_n\vert \rho \vert \phi_n\rangle, \end{equation} \tag{ 9 }$$

and extract information about the spectrum of this data

$$\begin{equation} G(w):= \mathcal{F}[g] = \sum_i \delta(w-\phi_i) \langle \phi_i\vert \rho \vert \phi_i\rangle. \end{equation} \tag{ 10 }$$

In the original phase estimation protocol—illustrated in figure 2—this process is performed coherently. First one prepares the desired input state $\rho = \rho_I$ on an input register I, and a coherent superposition of time states $\sum_{t = 1}^T\vert t\rangle$ on a control register C. The time register then controls the unitary time evolution on register I, to obtain

$$\begin{equation} \rho_{CI} \propto \sum_{tt^{^{\prime}}} \vert t\rangle\langle t^{^{\prime}}\vert_C \otimes U(t) \rho_I U(t^{^{\prime}})^\dagger. \end{equation} \tag{ 11 }$$

Tracing out the input register then yields time series data, namely

$$\begin{equation} \mathrm{Tr}_I(\rho_{CI}) =: \rho_C \propto \sum_n \sum_{tt^{^{\prime}}} \vert t\rangle\langle t^{^{\prime}}\vert_C \mathrm{e}^{-\mathrm{i} \phi_n t} \langle \phi_n\vert\rho_I\vert \phi_n\rangle \mathrm{e}^{\mathrm{i} \phi_n t^{^{\prime}}} = \sum_{t t^{^{\prime}}}g(t^{^{\prime}}-t) \vert t\rangle\langle t^{^{\prime}}\vert_C . \end{equation} \tag{ 12 }$$

The QFT is applied, and we obtain

$$\begin{equation} \textrm{QFT}(\rho_C) \propto \sum_{w w^{^{\prime}}}\left(\sum_{t t^{^{\prime}}} \mathrm{e}^{-\mathrm{i} w t + \mathrm{i} w^{^{\prime}}t^{^{\prime}}} g(t^{^{\prime}}-t)\right) \vert w\rangle\langle w^{^{\prime}}\vert_C. \end{equation} \tag{ 13 }$$

Measuring the control qubits in the computational basis is equivalent to sampling from G(w) over possible n-bit values of w.

Zoom In Zoom Out Reset image size

Applying QPE to problems involving local Hamiltonians differs significantly from its application to factoring in two important ways. First, it is unlikely that one can construct a polynomial depth circuit decomposition of U(T) for exponentially large T [AA17].

Second, the size of a problem involving a local Hamiltonian goes as the number of qubits. As the number of qubits grows, the number of eigenvalues increase exponentially, but the spectral range of the Hamiltonian increases polynomially ⁴ . It is the spectral range of the Hamiltonian which determines the bandwidth of g(t), and so it is not clear that the inclusion of a QFT should yield significant gains for such problems. If on the other hand one wishes to resolve fine-grained details of G(w), then from a classical signal processing perspective the parameter which dictates this resolution is the maximum simulated time T, not the number of samples of g(t).

This motivates an alternative method, which avoids the use of an expensive QFT. Instead of preparing and processing a coherent superposition of g(t) sampled at an exponential number of time steps, g(t) is directly measured for a polynomial number of time steps up to a maximum time T, and this data is processed classically. The procedure for measuring g(t) (illustrated in figure 2) is a straightforward phase kickback protocol, where the control register is prepared in the $\vert +\rangle$ state, the input register is prepared with the input state ρ_I , and a controlled time evolution is applied

$$\begin{equation} \rho_{CI} = \frac{1}{2} \sum_{ij=0}^1 \vert i\rangle\langle j\vert \otimes U(t)^i \rho_I U(-t)^j. \end{equation} \tag{ 14 }$$

Measuring the expectation value of X + iY on the control register yields:

$$\begin{align} \langle X_C \rangle + i \langle Y_C \rangle = \mathrm{Tr}[\rho_I U(t)] = \sum_i \mathrm{e}^{- i \phi_i t} \langle \phi_i\vert\rho_I \vert \phi_i\rangle = g^*(t) = g(-t). \end{align} \tag{ 15 }$$

Equipped with a quantum black box for values of g(t) it would appear that the remainder of the problem is merely a classical signal processing problem, about which a rich body of literature exists. However the existing classical signal processing literature typically assumes that the signal is composed of a polynomial number of frequencies, which is not generally true for a signal as in equation (9) composed of up to exponentially many distinct eigenvalues of the Hamiltonian of interest.

As outlined in the last section, resolving all exponentially many eigenvalues of a Hamiltonian is costly, as it would require us to simulate an exponential time evolution under H. Yet it is conceivable that one can obtain approximate information about the spectrum in a more efficient manner. Such an algorithm was introduced by [Som19] and we recap it here. The QEEP—as defined by [Som19]—is just such an approximation. In contrast to the ambitious goal of obtaining exponentially many eigenvalues for an operator H, the range of possible eigenvalues is divided up into bins, and the aim is to estimate the fraction of eigenvalues in each bin up to some smoothing on the boundaries of the bin.

Quantum eigenvalue estimation problem (QEEP)
Input:	A Hamiltonian H with associated eigenvalues $\{\lambda_i\}$ and eigenvectors $\{\vert \lambda_i\rangle\}$.
	An input quantum state ρ. A bin width η > 0, precision parameter ε > 0 and confidence level c < 1 5 . An indicator function $f(w): \mathbb{R} \rightarrow [0,1]$ s.t. $\forall w \not\in [-\eta, \eta]: \; f(w) = 0$ and $\forall w \in [0, \eta]: \; f(w)+f(w-\eta) = 1$.
Promise:	$||H||_{\infty} \leqslant 1/2$.
Question:	A vector $\textbf{q} = (q_0, q_1,{\ldots}, q_M)$ with $M = \lfloor 1/\eta \rfloor$ s.t. with probability $\geqslant c$, $||\textbf{q}-\textbf{p}||_1 \leqslant \epsilon$. Here $p_j = p(w_j): = \sum_{i} f(\lambda_i-w_j) \langle \lambda_i\vert\rho\vert \lambda_i\rangle$, and $w_j = j\eta-1/2$.

More specifically if one sets the input state to be the maximally mixed state this becomes $p_j = p(w_j): = \sum_{i} f(\lambda_i-w_j) 2^{-Q}$, where Q is the number of qubits.

The indicator function f in QEEP may be thought of as delineating the 'shape' of a bin. The values p_j can be thought of as the probability of observing the initial state ρ within a certain energy band, centred around $j \eta - 1/2$, whose resolution and sharpness is determined by the indicator function f. Equivalently, one may understand the values p_j as evenly spaced samples of G(w) convolved with f:

$$\begin{equation} p(w)=[G * f\,](w). \end{equation} \tag{ 16 }$$

Depending on the choice of f, solutions to a QEEP instance might have varying applications, and may be easier or harder to solve. A general scheme for constructing a function f satisfying the conditions outlined in the problem statement is to take the ideal indicator function

$$\begin{equation} \textrm{rect}_\eta(w) = \begin{cases}1 & w \in [-\eta/2, -\eta/2] \\ 0 & \textrm{otherwise} \end{cases} \end{equation} \tag{ 17 }$$

and convolve it with a normalized smoothing kernel $h(w) \in [0, 1]$ with support only in $[-\eta/2, \eta/2]$ ⁶ to obtain

$$\begin{equation} f(w) = [\textrm{rect}_\eta * h](w). \end{equation} \tag{ 18 }$$

The algorithm that Somma proposes for resolving QEEP relies on the following signal processing principle. For a well chosen smoothing kernel h, the Fourier transform $P(t): = \mathcal{F}[p](t)$ of p(w) may be made to fall off rapidly beyond a maximum time T, by smoothing out the sharply varying spectral function G(w). Furthermore, since p(w) has bounded support, P(t) need only be sampled at a finite rate. Thus, one need only sample g(t) at a finite rate to a maximum time T to retrieve a good approximation to the values p_j .

More precisely, we may consider that since G(w) and f(w) have bounded support they may be taken to be 2π periodic ⁷ and decomposed into Fourier series expansions:

$$\begin{align} G(w) =\frac{1}{2 \pi} \sum_{t = - \infty}^{\infty} g_{-t} \mathrm{e}^{\mathrm{i} t w} \quad \text{with} \quad g_t = g(t) \end{align} \tag{ 19 }$$

$$\begin{align} f(w) = \sum_{t = - \infty}^{\infty} F_t \mathrm{e}^{\mathrm{i} t w} \quad \text{with} \quad F_t = \mathcal{F}[\,f\,](t) \end{align} \tag{ 20 }$$

and thus:

$$\begin{equation} p(w) = [G * f\,](w) = \int_{-\infty}^{\infty} G(w^{^{\prime}}) f(w-w^{^{\prime}}) \mathrm{d}w^{^{\prime}} = \sum_{t = - \infty}^{\infty}\mathrm{e}^{\mathrm{i} t w} F_t g_{-t}. \end{equation} \tag{ 21 }$$

For fixed values $w_j = -1/2+j\eta$ we retrieve

$$\begin{equation} p_j = \sum_{t = - \infty}^{\infty} F_t^j g_{-t} \;,\; F_t^j=\mathrm{e}^{\mathrm{i} t w_j} F_t. \end{equation} \tag{ 22 }$$

In Somma's construction the indicator function f(w) is defined as in equation (18), with a smoothing kernel $h(w) = a \exp(-1/(1-(cw)^2) )$. With appropriate scaling (correct choice of c and normalization a ≈ 2.25). It is then argued that for any bin width η there exists some constant α > 1 such that for all $|t| \geqslant \alpha/\eta$

$$\begin{equation} |F_t| \leqslant \eta \exp\left(- \sqrt{|t| \eta/2}\right). \end{equation} \tag{ 23 }$$

Thus truncating the sum in equation (22) up to some maximum value $T \geqslant \alpha/\eta$ yields an approximation

$$\begin{equation} p_j^{^{\prime}} = \sum_{t = - T}^{T} F\,_t^j g_{-t} \end{equation} \tag{ 24 }$$

such that $\| \mathbf{p^{^{\prime}}} - \mathbf{p}\| \leqslant \frac{\epsilon}{2}$. We include an explicit computation of this maximum value T in appendix B. However we emphasize that this is an asymptotic result, and that finding α is a non-trivial task. For this reason we also consider a different indicator function in section B.1.

In addition to the source of error from this truncation, there is a source of error stemming from imperfect measurements of g(t). It is assumed that the T-dimensional vector of actual measurements $\tilde{\textbf{g}}$ satisfies $\| \tilde{\textbf{g}}-\textbf{g} \|_1 \leqslant \epsilon/2$ with probability greater than the confidence parameter c. Thus one retrieves a vector q given by

$$\begin{equation} q_j = \sum_{t = - T}^{T} F\,_t^j \tilde{g}_{-t} \end{equation} \tag{ 25 }$$

such that $\| \textbf{q}-\textbf{p} \|_1 \leqslant \epsilon$ with probability greater than c. Further details concerning this can be found in [Som19]. Importantly c depends on the number of repetitions of the circuit depicted in figure 2; and not the maximum evolution time T.

Choosing the 1-norm yields bounds on the overall deviation from the ideal vector p. One may consider other norms, such as the max norm, which bounds the maximum deviation, and yields a milder variant of Somma's error bound. For completeness and for the purpose of deriving explicit constants in the error bounds we walk through Somma's error derivation in detail in appendix B.

As discussed in the introduction, a fundamental numerical primitive is to ask for eigenvalues for a given normal or Hermitian operator A. More specifically, given a Hermitian matrix A, the task is to return its spectrum—the set of its eigenvalues—up to some bit precision b, i.e. numerical accuracy $\xi = 2^{-b}$. Instead of speaking of the numeric precision with which results are obtained, one can take the viewpoint that such a function sorts A's eigenvalues into 2^b many bins of width $\propto \xi$. Such a function is e.g. implemented in matlab's eig or Mathematica's Eigenvalues functions. Assuming no other numerical errors are present during the calculation ⁸ , we expect the result to be such that eigenvalues are rounded to their nearest bin.

As currently formulated, the QEEP contains overlapping indicator functions $f^j(\omega) : = f(\omega - \omega_j)$, so that an eigenvalue in the support of both f^j and $f\,^{j+1}$ will contribute to both p_j and $p_{j+1}$. There is an important reason for formulating QEEP in this manner. Consider what a 'perfect' choice of non-overlapping $f^j(\omega)$ might look like. As described in [Som19], the natural choice would be an indicator function such as

$$\begin{align} f(\omega) = \begin{cases} 1 & - \eta/2 \leqslant \omega < + \eta/2 \\ 0 & \text{otherwise}. \end{cases} \end{align} \tag{ 26 }$$

This choice would resolve the spectrum perfectly; and furthermore remove $f(\omega)$ from the input of the problem statement. However, as explained by the author of [Som19] the slow decay of the Fourier coefficients would mean one would have to simulate $\mathrm{e}^{-\mathrm{i} T H}$ for impractically large T. Due to this one must choose a smooth $f(\omega)$ which approximates this ideal solution through having smoothly decaying support over $-\eta \leqslant \omega\lt + \eta$. That is to say the indicator functions overlap each other so that $\forall w \in [0, \eta]: \; f(w)+f(w-\eta) = 1$.

To address these considerations we formulate a problem which we call the rQEEP for short, which is close in spirit to Matlab's eig or Mathematica's Eigenvalues function. For a given H rQEEP outputs a binned estimate of the spectrum up to a probabilistic guarantee on the likelihood that eigenvalues close to bin boundaries have been misplaced; and which replaces the issue overlapping bins with randomly defined bins. We will now define this problem formally and then propose an algorithm to solve rQEEP via a reduction to QEEP.

Randomized quantum eigenvalue estimation problem (rQEEP)
Input:	Local Hamiltonian $H = \sum_i h_i$. Maximal point spacing ξ > 0, deviation $\Delta\geqslant0$, and confidence $c \in [0, 1]$.
Promise:	$\operatorname{spec}(H) \subseteq [-1/2,1/2]$.
Output:	$m : = 2/\xi$. List of coordinates $(x_i, y_i)$ for $i = 1, \ldots, (m+1)$ with $x_0 : = -1/2$ and $x_{m+1} : = 1/2$, such that
	1. xi are randomly but uniformly spaced, i.e. for all $x\in[-1/2,1/2]$ there exists a xi s.t. $|x-x_i| \lt \xi/4$, and the distance between neighbouring points is bounded by $|x_i-x_{i-1}| \lt \xi$.
	2. The total absolute difference between the estimated (yi ) and true (ni ) number of eigenvalues of H within the interval $B_i = [x_{i-1}, x_i)$ is upper bounded—with probability greater than c—by the deviation Δ, i.e. $\sum_i |y_i - n_i| \leqslant \Delta$.

We first remark that the above problem is well-defined: indeed, we can always first obtain the entire spectrum of H to sufficient precision, and choose the abscissas such that the intervals B_i never have an eigenvalue on their boundaries (which is possible even with the restriction of $n = \mathrm{poly} 1/\xi$, as the spectrum of H is finite); and then define y_i as the number of eigenvalues within that interval. This answers the problem even for deviation $\Delta = 0$ and confidence c = 1.

Naturally, what we will be interested in is whether we can answer rQEEP efficiently, and for which sets of parameters we can do so.

Given a local Hamiltonian $H = \sum_i h_i$ acting on Q qubits (which we assume to be rescaled such that the spectrum lies within $[-1/2,1/2]$), we can solve the rQEEP by using QEEP as a subroutine and given access to a QEEP solver with a sub-circuit cost of $\mathcal{T}(Q,\eta,\epsilon,c)$. We do this via the following algorithm.

• (a)  
  Given a maximal point spacing ξ > 0, deviation $\Delta\geqslant0$, and confidence $c \in [0, 1]$ uniformly partition the interval $[-1/2,1/2]$ into $m: = 2/\xi$ many segments

  $$\begin{align} \left[ -\frac12, -\frac12 + \frac\xi2 \right) \cup \left[ -\frac12 + \frac\xi2, -\frac12 + \xi \right) \cup \ldots \cup \left[ \frac12 - \frac\xi2, \frac12 \right] =: S_1 \cup \ldots \cup S_m. \end{align} \tag{ 27 }$$
• (b)  
  With confidence c run the following subroutine $\lceil \log_{(1-\frac{c}{2})}(1-c) \rceil$ times, and collect the results in a list:

  • 1.  
    Within each segment, choose x_i uniformly at random from S_i . Set $x_0 : = -1/2$, and $x_{m+1} : = 1/2$.
  • 2.  
    Solve a QEEP instance with a suitable indicator function f, input state $\rho = I/2^Q$, bin width $\eta : = \Delta/(6\times 2^{Q+1}(m+1)^2)$, precision $\epsilon = \Delta /(8 \times 2^{Q+1})$, and confidence c. Denote the output with $\mathbf q = (q_0,\ldots,q_M)$ where $M = \lfloor 1/\eta \rfloor$, for the bins centred around $w_j = -1/2 + j\eta$ for $j = 0,\ldots,M$.
  • 3.  
    For all $i = 0,\ldots,m$, define lower and upper envelopes via

    $$\begin{align} y^\mathrm{lwr}_i &:= 2^Q \!\!\!\! \sum_{w_j \in R^\mathrm{lwr}_i} \!\!\!\! q_j &\text{and}&& y^\mathrm{upr}_i &:= 2^Q \!\!\!\! \sum_{w_j \in R^\mathrm{upr}_i} \!\!\!\! q_j \end{align} \tag{ 28 }$$

    $$\begin{align} R^\mathrm{lwr}_i &:= [x_i + \eta, x_{i+1} - \eta ] &\text{and}&& R^\mathrm{upr}_i &:= [x_i - \eta, x_{i+1} + \eta ] \end{align} \tag{ 29 }$$

    and where we implicitly assume $R^\mathrm{lwr}_i = \emptyset$ if $x_i + \eta \gt x_{i+1} - \eta$.
  • 4.  
    Return $y_i : = (y^\mathrm{lwr}_i + y^\mathrm{upr}_i) / 2$.
• (c)  
  Return the result from the subloop with the smallest deviation of lower and upper bounds.

This is illustrated in figure 3. In appendix A we prove the above algorithm has the following quantum overhead

Theorem 3. Given access to a QEEP solver with a subcircuit runtime $\mathcal{T}(Q,\eta,\epsilon,c)$ then the following holds. For a point spacing ξ > 0 (or equivalently $m : = 2/\xi$), deviation $\Delta\geqslant0$, and confidence c, the overall subcircuit runtime to answer rQEEP using QEEP is

$$\begin{align} \mathcal{T}\left(Q,\frac{1}{12 (m+1)^2} \times \frac{\Delta}{2^Q},\frac{1}{16} \times \frac{\Delta}{2^Q},c \right), \end{align} \tag{ 30 }$$

and the additional classical overhead is $\log_{(1-\frac{c}{2})}(1-c)$.

Zoom In Zoom Out Reset image size

Using theorem 3 we can solve for m and Δ in terms of ε, η and Q and find the following dependence

$$\begin{align} m &= \sqrt{\frac{4\epsilon}{3\eta}} -1 \end{align} \tag{ 31 }$$

$$\begin{align} \Delta &= 2^Q \times 16 \epsilon. \end{align} \tag{ 32 }$$

Equations (31) and (32) tells us the number of bins and deviation with which we can expect when solving an instance of rQEEP given access to a QEEP solver which can achieve and bin width of η and a precision ε. The question is whether this now allows us to solve a problem of interest. The total deviation Δ refers to maximum number of misplaced eigenvalues. From the above we see that the total fraction of misplaced eigenvalues is proportional to ε, even though the total number scales as 2^Q .

The runtime $\mathcal{T}$ depends on this η and ε. If we assume that we can solve QEEP efficiently for $\eta, \epsilon = \text{poly}(Q^{-1})$, then we can expect to efficiently solve rQEEP for a polynomially scaling number of bins $m = \text{poly}(Q)$ and a polynomially scaling fraction of misplaced eigenvalues, $\Delta/2^Q = \text{poly}(Q^{-1})$.

From [HMH04, HMH05] we know that for typical Hamiltonians the energy distribution will become Gaussian as $Q \rightarrow \infty$. Therefore we do not expect to be able to resolve features in the low energy subspace in this limit. Additionally, from [BFS11] we know that the computational difficulty of calculating the density of states for local Hamiltonians is in the counting version of QMA (known as #BQP) and so we cannot use rQEEP to do this efficiently.

However, we could still determine important features of the spectrum. The existence of a Mott insulator transition in Fermi-Hubbard models, and the Hamiltonian parameters for which this transition occurs, are of widespread interest [IFT98]. As the Mott insulator phase features a spectral gap we could try to observe this with rQEEP for a fixed number of qubits Q but while varying parameters such as the strength of on-site and hopping interactions.

In this section we propose an alternative protocol to obtain g_t using several geometrically local control ancillas in place of a single globally connected ancilla as described in [Som19].

We consider a reasonably general case. Suppose we wish to solve QEEP for a k-local qubit Hamiltonian, $H = \sum_{i} \alpha_i h_i$. As we are interested in the NISQ era algorithms we will assume that when laid out on hardware, every Pauli interaction h_i acts non-trivially on geometrically local qubits. Furthermore we assume that the device has sufficient connectivity so that we can introduce a geometrically local ancilla qubit for each of these interactions. We denote these ancillas by dashed indices, iʹ, in contrast to the data qubits with undashed indices. This is illustrated for a two-local Hamiltonian with two terms in figure 4.

Zoom In Zoom Out Reset image size

Over the joint Hilbert space we define an augmented Hamiltonian

$$\begin{align} H^{^{\prime}} := \sum_{i=1} \alpha_i Z_{i^{^{\prime}}} \otimes h_i . \end{align} \tag{ 33 }$$

To obtain time series data g_t from this augmented Hamiltonian, we first prepare a cat state on the control ancillary space, C. For $\rho_I = \vert \psi\rangle\langle \psi\vert_I$ the overall input is

$$\begin{align} \vert \psi(0)\rangle_{CI} = \frac{1}{\sqrt{2}} \big( \vert 00 \cdots 0\rangle_C + \vert 11 \cdots 1\rangle_C\big)\otimes \vert \psi\rangle_I. \end{align} \tag{ 34 }$$

On this initial state, the cat state induces a sign flip on the augmented Hamiltonian H^ʹ, conditioned on the cat state to be in state $\vert 1\rangle_C$; more precisely, applying the evolution operator $U^{^{\prime}}_c(t) : = \exp(\mathrm{i} t H^{^{\prime}}/2)$ we get

$$\begin{align} \vert \psi(t)\rangle_{CI} &= \frac{1}{\sqrt{2}} \left( \vert 00 \cdots 0\rangle_C \mathrm{e}^{\mathrm{i} t H/2} \vert \psi\rangle_I + \vert 11 \cdots 1\rangle_C \mathrm{e}^{- \mathrm{i} t H/2} \vert \psi\rangle_I \right), \end{align} \tag{ 35 }$$

as can be readily verified. By measuring the expectation value $g_t = \langle X_{1^{^{\prime}}} + \mathrm{i} Y_{1^{^{\prime}}} \rangle_{\rho(t)}$ for $\rho(t) = \vert \psi(t)\rangle\langle \psi(t)\vert_{CI}$ we obtain the time series data g_t , as desired.

In this section we will establish the run-time of the single most costly experiment involved in solving QEEP for $H = \sum_i \alpha_i h_i$. That is we will calculate the run-time—specified by definition 2—associated with the protocol to obtain g_T outlined in section 4.1.

In previous work [CBC21] we establish tighter bounds on the run-time of implementing $U(T) = \mathrm{e}^{- \mathrm{i} T H}$ where $H = \sum_i \alpha_i h_i$ is a k-local Hamiltonian with $\alpha_i \leqslant 1$ and $\|h_i\| = 1$. First we decompose U(T) via a Trotter approximation, that is set $U(T) \approx (\prod_i \mathrm{e}^{\mathrm{i} \delta \alpha_i h_i})^{T/\delta} + \mathrm{O}(T \delta)$, and then use the sub-circuit synthesis techniques given in [CBC21] to implement each local Trotter step $\mathrm{e}^{\mathrm{i} \delta \alpha_i h_i}$ with sub-circuit run-time $\mathcal{T}(\mathrm{e}^{\mathrm{i} \delta \alpha_i h_i})\approx \mathrm{O}(\delta^{1/(k-1)})$ ⁹ .

In order to perform QEEP with local control as described in the previous section we must instead implement $U^{^{\prime}}_{c}(T) = \mathrm{e}^{- \mathrm{i} T H^{^{\prime}}}$ where $H^{^{\prime}} = \sum_i \alpha_i Z_{i^{^{\prime}}} h_i$ in the manner described above. Now the run-time of each local Trotter step is $\mathcal{T}(\mathrm{e}^{\mathrm{i} \delta \alpha_i Z_{i^{^{\prime}}} h_i})\approx \mathrm{O}(\delta^{1/k})$. As there are $T/\delta$ Trotter steps this means

$$\begin{align} \mathcal{T}(U^{^{\prime}}_{c}(T)) \leqslant \mathrm{O}\left( T \delta^{\frac{1-k}{k}} \right). \end{align} \tag{ 36 }$$

This diverges as $\delta \rightarrow 0$. Therefore, we must choose the largest δ for a given Trotter error ε ¹⁰ . That is we choose the value of δ which saturates the inequality $\epsilon \geqslant \mathrm{O}(T \delta)$. Substituting this into equation (36) gives

$$\begin{align} \mathcal{T}(U^{^{\prime}}_{c}(T)) \leqslant \mathrm{O} \left(T^{1 + \frac{k-1}{k}} \epsilon^{-\frac{k-1}{k}} \right). \end{align} \tag{ 37 }$$

As discussed in section B.1, T will depend in the bin-width η and precision. We will use the bound equation (69) to determine this, meaning

$$\begin{align} \mathcal{T}(U^{^{\prime}}_{c}(T)) \leqslant \mathrm{O} \left(T^{1 + \frac{k-1}{k}} \epsilon^{-\frac{k-1}{k}} \right) \quad \text{with} \quad T \rightarrow \frac{\pi}{\eta} \times \frac{\mathrm{e}^{\pi \eta \epsilon/2}}{\sqrt{\mathrm{e}^{\pi \eta \epsilon}-1}}. \end{align} \tag{ 38 }$$

This run-time does not take any connectivity assumptions into account, however we will address the issue of connectivity in our numeric benchmarking of this approach in section 4.3. As we are interested in the NISQ regime we will focus on a numerical analysis of a specific example Hamiltonian, the 2D spin Fermi-Hubbard Hamiltonian, rather than an asymptotic analysis of the general algorithm.

We will take the following problem instance; a rQEEP input Hamiltonian which is a two-dimensional L × L lattice Fermi-Hubbard model with spin. This Hamiltonian is made up of on-site interactions and nearest neighbour hopping terms

$$\begin{align} H &:= \sum_{i=1}^{L^2} h_{\text{on-site}}^{(i)} \ + \sum_{i < j,\sigma} h_{\text{hopping}}^{(i,j,\sigma)} \end{align} \tag{ 39 }$$

$$\begin{align} \kern84.4pt &:= u \sum_{i=1}^{L^2} a^{\dagger}_{i \uparrow} a_{i \uparrow} a^{\dagger}_{i \downarrow} a_{i \downarrow} + v \sum_{i < j,\sigma} \left(a^{\dagger}_{i \sigma} a_{j \sigma} + a^{\dagger}_{j\sigma} a_{i \sigma}\right) \end{align} \tag{ 40 }$$

which describe electrons; and with spin $\sigma = \ \uparrow$ or $\downarrow$ hopping between neighbouring sites on a lattice, with an on-site interaction between opposite-spin electrons at the same site. We will assume that this fermionic Hamiltonian has been encoded into a qubit Hamiltonian using the compact fermion encoding [Der+21]. This first requires us to decide on the layout of qubits, as shown in figure 5. This figure shows the layout of a single spin sector for L = 3 on a device with planar connectivity. First the on-site interaction becomes,

$$\begin{align} h_{\text{on-site}}^{(i)} \rightarrow \frac{u}{4} \left(I - Z_{i \uparrow}\right)\left(I - Z_{i \downarrow}\right). \end{align} \tag{ 41 }$$

Zoom In Zoom Out Reset image size

The exact expressions for hopping interactions are more complicated and depend on whether two nearest neighbour fermionic sites are horizontally or vertically connected on the lattice [Der+21]. We also need ancilla qubits labelled with a dashed index and shown in blue in figure 5. They are encoded as

$$\begin{align} h_{\text{hopping}}^{(i,j,\sigma)} \rightarrow \begin{cases}\frac{v}{2} \left(X_{i,\sigma}X_{j,\sigma}Y_{f_{ij}^{^{\prime}},\sigma}+Y_{i,\sigma}Y_{j,\sigma}Y_{f_{ij}^{^{\prime}},\sigma}\right) & \text{horizontal} \\ \frac{v}{2} (-1)^{g(i,j)}\left( X_{i,\sigma}X_{j,\sigma}X_{f_{ij}^{^{\prime}},\sigma}+Y_{i,\sigma}Y_{j,\sigma}X_{f_{ij}^{^{\prime}},\sigma}\right) & \text{vertical} \end{cases}. \end{align} \tag{ 42 }$$

In this notation i labels the data qubit for lattice site i and σ its spin lattice. Dashed indices such as $f_{ij}^{^{\prime}}$ refer to ancillary qubits for the fermionic encoding. The function $g(i,j)$ is either zero or one and depends on the details of the compact encoding [Der+21].

The locality of these Hamiltonian terms will increase when we introduce further ancilla qubits required to perform QEEP. These are shown in purple in figure 5. With planar connectivity this layout avoids the introduction of any SWAP gates, if one were only interested in one spin sector. To avoid SWAP overheads we therefore need a device with 3D connectivity. We must assume the device consists of two stacked and connected 2D grids of qubits with nearest neighbour connectivity, where one of these 2D grids is shown in figure 5.

To work out the individual costs of these now three and four local gates, we reference [CBC21, appendix B], and for an application of either term to time δ we obtain

$$\begin{align} \kern29.7pt \text{three-local on-site:}&\quad g_3(u, \delta) = g_3(u\delta) = \frac{|u \delta|}{2} + 2 \sqrt{2 |u \delta|} \end{align} \tag{ 43 }$$

$$\begin{align} \text{four-local hopping:}&\quad g_4(v, \delta) = g_4(v\delta) = 12 \sqrt[3]{2 |v \delta|} \end{align} \tag{ 44 }$$

for a first order Trotter formula. We also remark that the compact encoding requires a single layer of the on-site terms, but four layers for the hopping terms; the total cost for a single Trotter step with time δ we have $g_3(u\delta) + 4g_4(v\delta)$. Additionally we assume the noise model explained in section 2.1 applies to this device and use equation (5) to calculate the run-time region where the combined stochastic and analytic error remains bounded by $\epsilon_\mathrm{tar}$ for a variety of noise parameters q.

Now we will proceed to analyse QEEP and then rQEEP for a specific problem instance assuming this device connectivity. We will assume that our device has a given run-time budget before it is overcome by noise. We will determine numerically how the minimum obtainable bin-width η goes as a function of this run-time budget for various values of the target precision. That is we will perform the numerical equivalent of the asymptotic analysis of section 4.2 but solving equation (38) to get $\eta = \eta(\mathcal{T}, \epsilon)$.

In order to solve QEEP with a bin-width of η we must rescale H in order to apply the bound equation (73) or equation (69). For H on an L × L lattice, we know that $\|H\| \leqslant 5\Lambda$, where Λ is an upper bound on the norm of any of the groups of commuting interactions in the Trotter decomposition. We thus need to rescale the maximum evolution time by a factor $1/10\Lambda$. The maximum evolution time T is then bounded by

$$\begin{align} \text{Equation (73)}:&\quad T \geqslant \frac{\pi}{10\eta\Lambda} \times \frac{\mathrm{e}^{\pi \eta \epsilon/2}}{\sqrt{\mathrm{e}^{\pi \eta \epsilon}-1}}, \end{align} \tag{ 45 }$$

$$\begin{align} \kern36.6pt \text{Equation (69)}:&\quad T \geqslant \frac{1}{5\eta\Lambda}\left(1+\operatorname{plog}_{-1}\left( - \frac{\epsilon\eta}{32e} \right)\right)^2, \end{align} \tag{ 46 }$$

where we are considering two possible choices of $f(\omega)$ with each choice leading to a different bound on the maximum evolution time as derived in appendix B. We assume $U^{^{\prime}}_c(T)$ is implemented by a pth-order Trotter formula using five layers. Each layer is made up of mutually commuting interactions which are also disjoint. We assume the NISQ device can implement these interactions in parallel (see [CBC21] for more details). The function $\eta = \eta(\mathcal{T}, \epsilon)$ has a further dependence on which order Trotter formula is used to implement $U^{^{\prime}}_c(T)$ as this determines $\mathcal{T}$. However in our numerical analysis we will simply show the minimum η obtained by varying over $p = 1, 2$ or 4. Furthermore we will use the numerical Trotter bounds obtained in [CBC21], based on extrapolating optimal Trotter step sizes from lattice sizes up to $3\times3$.

To benchmark the impact of our sub-circuit synthesis techniques we perform the same analysis but restricting ourselves to decomposing the local unitaries into CNOT gates and single qubit rotations. We cost these gates in the same framework, meaning that in this case the run-time of each local Trotter step equals the circuit depth (disregarding single-qubit gates) of the decompositions of the three- and four-local interactions using the lowest-depth gate decomposition method (via conjugation), which yields

$$\begin{align} \kern5.5pt \text{three-local on-site:}&\quad g^{^{\prime}}_3(u, \delta) = g^{^{\prime}}_3 = \frac{5 \pi}{4} \end{align} \tag{ 47 }$$

$$\begin{align} \text{four-local hopping:}&\quad g^{^{\prime}}_4(v, \delta) = g^{^{\prime}}_4 = \frac{5 \pi}{2} \end{align} \tag{ 48 }$$

and again that a single layer has a runtime $g^{^{\prime}}_3 + 4g^{^{\prime}}_4$. The results of these numerics are presented in figure 1 and show that–within the parameter range considered–our asymptotically suboptimal choice of indicator function (solid lines) and synthesis techniques (red lines) yield improvements in the minimal achievable bin-width η for a given depolarizing noise rate and total runtime $\mathcal{T}$.

We have built on our previous work, recapped in section 2, where we defined a non-standard circuit error model from an abstraction of the capabilities of NISQ hardware and then introduced gate synthesis techniques tailored to this framework. Here we have applied this framework to the problem of QEEP—introduced in [Som19]—to determine whether our gate synthesis techniques yield any improvements within this circuit error model.

In order to do this we have outlined a simple implementation of QEEP in section 4 aimed at minimizing runtime by performing controlled unitary evolution using only subcircuit gates. This QEEP implementation preserves any potential subcircuit gains, shown in [CBC21], by exchanging the global control ancilla of the QEEP implementation in [Som19] with local control ancillas and the creation of a global CAT state, in order to highlight its runtime minimizing properties within our minimize the runtime in the per-time error model. Note that we do not include the cost associated with measurement in our work, however this overhead will be a constant addition to both standard and subcircuit approaches.

It is essentially true by the definition of our per-time error model that our sub-circuit synthesis techniques will have an asymptotic advantage over standard synthesis techniques in this framework, as is shown in section 4.2. In addition to this, asymptotic analysis is of limited use within a NISQ context, as there may be large but hidden constants subsumed into big O notation. For both of these reasons we have anchored our conclusions in a numerical analysis of the 2D spin Fermi-Hubbard model. We have considered how implementing QEEP with local control relates to the connectivity of a NISQ device (see figure 5) and then perform numerics to obtain the results shown in figure 1. These plots demonstrate a numeric advantage for our synthesis techniques within a per-time error model. Modulo the applicability of a per-time error model, these techniques may allow one to perform a two bin instance of QEEP on an L = 5 spin Fermi-Hubbard model on hardware with a noise rate of $q = 10^{-5}$ while keeping the total error under $10\%$.

This encapsulates the first avenue of research explored in this paper, a development and application of our work in [CBC21] to a assess the possibility of performing QEEP in the NISQ-era.

Secondly, we have reframed the QEEP protocol presented in [Som19] as an algorithm which allows one to sample from a convolution between the spectrum of a Hamiltonian, $G(\omega)$ and a freely chosen function $f(\omega)$ ¹¹ . This is explained in section 2.3 and specifically via equation (16). By recasting the problem and algorithm in this way we hope to clarify which aspects we can modify and adapt to different computational problems.

Most trivially our analysis decouples the precision parameter ε from the bin-width η. However we have also explored the possibility of tailoring our choice of $f(\omega)$ to problem and parameter regime of interest, in our case an instance of QEEP with a Fermi-Hubbard Hamiltonian in a NISQ regime.

We define one such alternative indicator function in appendix B and compare the resultant lower bound on the maximum evolution time T with the original indicator function of [Som19] in figure 6. We emphasize that the original function gives asymptotically superior results, as it was intended to. Despite this, within the non-asymptotic parameter range of figure 1 we can see that for both standard and subcircuit synthesis techniques our bound gives improvements in the minimum obtainable η. For example, we estimate that the depolarizing single qubit error rate required to implement a minimum two bin example of QEEP—with a $5\times5$ Fermi-Hubbard model and up to a precision of ε = 0.1—can be reduced from 10⁻⁷ to 10⁻⁵.

Zoom In Zoom Out Reset image size

As we see from the above, the optimal choice of indicator function may depend on the task at hand. It is also possible that other applications might benefit from taking non-uniform samples from the time series g(t) and employing more sophisticated classical signalling processing methods to reconstruct $p(\omega) = [G*f](\omega)$. By formulating QEEP as a protocol for sampling from a $p(\omega)$ we hope to make it easier to explore these questions in future work.

Finally, we suggest an application of QEEP having thus reformulated it as a general tool. We consider a problem analogous to matlab's eig or Mathematica's Eigenvalues. This was outlined in section 3.1 and in section 3.2 we demonstrate how to solve it using QEEP as a subroutine. We also incorporate this problem into our numerics as shown in figure 1 and find that it is more demanding than performing QEEP alone—as expected, due to the superior guarantees rQEEP gives on the placement of eigenvalues within each bin.

Computing an explicit lower bound for T using Somma's original indicator function is a non-trivial task. We include an asymptotic calculation in appendix B.2. However for near term applications, a simpler indicator function may yield better analytic bounds in the non-asymptotic regime. For this reason we consider a different indicator function, which allows us to derive an explicit bound on T without further restrictions when said bound holds. To this end—and re-using the symbols f_j and F_j in the following—we define

$$\begin{equation} f(w) := \left\{ \begin{array}{ll} \cos^2\left(\frac{\pi w}{2\eta}\right) & -\eta < w < \eta \\ 0 & \text{otherwise}. \end{array} \right. \end{equation} \tag{ 65 }$$

It is easy to check that these indicator functions satisfy the conditions of QEEP. Furthermore, we can explicitly derive its Fourier transform (and thus the Fourier series coefficients) from equation (20), namely

$$\begin{align} F(t) = \left\{ \begin{array}{ll} \displaystyle\frac{\pi \sin(t\eta)}{2 t \pi^2 - 2 t^3 \eta^2} & t \neq 0 \\[5mm] \displaystyle\frac{\eta}{2\pi} & \text{otherwise}\end{array} \right. \end{align} \tag{ 66 }$$

and thus for $t\gt\pi/\eta$ we have

$$\begin{align} |F(t)| \leqslant \frac{1}{2\pi} \times \frac{1}{t^3\eta^2/\pi^2 - t}. \end{align} \tag{ 67 }$$

Together with equation (24) we then get

$$\begin{align} |p_j - p^{^{\prime}}_j| \leqslant \frac{1}{2\pi} \sum_{|t|>T} \frac{\pi^2}{t^3\eta^2 - t \pi^2} \leqslant \frac{1}{2\pi} \int_{T}^\infty \frac{\pi^2}{t^3\eta^2 - t \pi^2} \mathrm dk = \frac{1}{4\pi} \log \left(\frac{\eta ^2 T^2}{\eta ^2 T^2-\pi ^2}\right) \overset{!}{\leqslant} \frac{\epsilon}{2(M+1)} \end{align} \tag{ 68 }$$

which together with $M+1 \leqslant 2/\eta$ yields

$$\begin{equation} T \geqslant \frac{\pi}{\eta} \times \frac{\mathrm{e}^{\pi \eta \epsilon/2}}{\sqrt{\mathrm{e}^{\pi \eta \epsilon}-1}}. \end{equation} \tag{ 69 }$$

Both for equations (69) and (73) we gain a factor of M in case we are interested in the max-norm instead of the one-norm, as mentioned. We compare this new bound with Somma's bound in figure 6.

In this section we find a lower bound on the Fourier series cutoff T in equation (24), when using Somma's indicator function, as described in section 2.3. We wish to determine the minimal T such that

$$\begin{align} \|\mathbf{p} - \mathbf{p}^{^{\prime}}\| \leqslant \frac{\epsilon}{2}, \end{align} \tag{ 70 }$$

where p is given in equation (22). By equations (24) and (22), we have

$$\begin{align} |p_j - p^{^{\prime}}_j| \leqslant \sum_{|t| > T} |F\,_t^j| |g_t|. \end{align} \tag{ 71 }$$

Since $|g_t| \leqslant 1$, it suffices to find a bound $F_t^{\,j}$. Assuming that T is chosen large enough such that equation (23) holds gives us

$$\begin{equation} |p_j - p^{^{\prime}}_j| \leqslant \sum_{|t|>T} \eta \left|\exp\left(- \sqrt{|t| \eta/2 }\right)\right| \leqslant 2 \eta \int_{t=T}^\infty \exp\left(- \sqrt{|T| \eta/2 }\right) \mathrm dt = 8 \mathrm{e}^{-\sqrt{T\eta/2}}\left(1 + \sqrt{T \eta/2} \right). \end{equation} \tag{ 72 }$$

Distributing $\epsilon/2$ equally amongst the $M+1 \leqslant 2M$ bins gives the bound on the Fourier truncation T for $T\geqslant \alpha/\eta$

$$\begin{equation} \mathrm{e}^{-\sqrt{T\eta/2}}\left(1 + \sqrt{T \eta/2} \right) \overset{!}{\leqslant} \frac{\epsilon}{32M} \quad\Longrightarrow\quad T \overset!\geqslant \frac2\eta \left( 1 + \operatorname{plog}_{-1}\left(-\frac{\epsilon\eta}{32e} \right) \right)^2. \end{equation} \tag{ 73 }$$

Here the product logarithm $\operatorname{plog}_{-1}(y)$ denotes the −1th solution for w in the equation $y = w \mathrm{e}^w$, this is also known as the Lambert W function. We emphasise that this is an asymptotic result, finding α is not straightforward.