$p_{\nu}=\sum_{\mu=2}^{n_{p}+1}\sum_{\nu\in
B_{\mu}}\frac{\lceil{\mathcal{M}(2^{\mu-2}/\left\lVert\nu\right\rVert)^{2}}\rceil}{\mathcal{M}2^{2\mu}2^{n_{\mu}+1}}.$
(106)

As discussed in (Babbush  _et al._ , 2019, Section II D), this is a constant
around $1/4$ provided the ceiling function is not used. In the case where the
ceiling function is used, then the probability of success is only increased.

To bound the error due to the ceiling function, we can write each weight that
is prepared as $1/\left\lVert\nu^{\prime}\right\rVert^{2}$, for a modified
$\nu^{\prime}$. Now it would appear that all coefficients are larger, because
of the ceiling function. However, when accounting for the normalization of the
state, that will be larger, so each coefficient is multiplied by a constant
that is smaller than 1, resulting in approximate coefficients that may be
larger or smaller than the exact coefficients. In fact, because we are
choosing the amplitude on the potential terms, we have freedom to adjust it to
minimise the error. Let us therefore call the constant factor $\alpha$. The
squared amplitudes in the state we have prepared correspond to approximately

$\frac{1}{2^{n_{p}+6}}\frac{1}{\left\lVert\nu\right\rVert^{2}}.$ (107)

Therefore, let us put the values of $\nu^{\prime}$ as

$\frac{1}{2^{n_{p}+6}}\frac{1}{\left\lVert\nu^{\prime}\right\rVert^{2}}=\alpha\frac{\lceil{\mathcal{M}2^{\mu-2}/\left\lVert\nu\right\rVert)^{2}}\rceil}{\mathcal{M}2^{2\mu}2^{n_{p}+2}},$
(108)

so

$\frac{1}{\left\lVert\nu^{\prime}\right\rVert^{2}}=16\alpha\frac{\lceil{\mathcal{M}2^{2\mu}/(16\left\lVert\nu\right\rVert^{2})}\rceil}{\mathcal{M}2^{2\mu}}.$
(109)

Using this approach the error in the eigenvalues due to the approximation of
the state preparation may be bounded as follows.

######  Lemma 3.

Given that the state preparation for $\nu$ is performed as described above
with a given $\mathcal{M}$, if $\lambda_{j}$ is the $j^{\rm th}$ eigenvalue of
the true Hamiltonian then there exists an eigenvalue of the approximated
Hamiltonian $\widetilde{\lambda}_{j}$ such that

$\displaystyle|\lambda_{j}-\widetilde{\lambda}_{j}|$
$\displaystyle\leq\frac{2\eta}{\pi\mathcal{M}\Omega^{1/3}}\left(\eta-1+2\lambda_{\zeta}\right)\left(7\times
2^{n_{p}+1}-9n_{p}-11-3\times 2^{-n_{p}}\right).$
$\displaystyle=\frac{2^{n_{p}}(28\eta)\left(\eta-1+2\lambda_{\zeta}\right)}{\pi\mathcal{M}\Omega^{1/3}}\left(1+\mathcal{O}(2^{-n_{p}})\right)$

###### Proof.

Again we use the fact that the norm of the difference in the operators can be
used to upper bound the difference in the eigenvalues. The error in the
approximation of $U$ can be bounded as

$\displaystyle\left\lVert U-\widetilde{U}\right\rVert$
$\displaystyle\leq\left\lVert\sum_{\nu\in
G_{0}}\sum_{\ell=1}^{L}\frac{\zeta_{\ell}}{\pi\Omega^{1/3}}\left(\frac{1}{\left\lVert\nu\right\rVert^{2}}-\frac{1}{\left\lVert\nu^{\prime}\right\rVert^{2}}\right)\sum_{j=1}^{\eta}e^{-\mathrm{i}k_{\nu}\cdot
R_{\ell}}\bigg{(}\sum_{p\in G}\lvert p-\nu\rangle\langle
p\rvert_{j}\bigg{)}\right\rVert$ $\displaystyle\leq\sum_{\nu\in
G_{0}}\frac{\eta\lambda_{\zeta}}{\pi\Omega^{1/3}}\left|\frac{1}{\left\lVert\nu^{\prime}\right\rVert^{2}}-\frac{1}{\left\lVert\nu\right\rVert^{2}}\right|$
$\displaystyle=\frac{\eta\lambda_{\zeta}}{\pi\Omega^{1/3}}\sum_{\mu=2}^{n_{p}+1}\sum_{\nu\in
B_{\mu}}\left|\frac{1}{\left\lVert\nu^{\prime}\right\rVert^{2}}-\frac{1}{\left\lVert\nu\right\rVert^{2}}\right|.$
(110)

Similarly, the error in the approximation of $V$ can be bounded as

$\displaystyle\left\lVert V-\widetilde{V}\right\rVert$
$\displaystyle\leq\left\lVert\sum_{\nu\in
G_{0}}\frac{1}{2\pi\Omega^{1/3}}\left(\frac{1}{\left\lVert\nu\right\rVert^{2}}-\frac{1}{\left\lVert\nu^{\prime}\right\rVert^{2}}\right)\sum_{\begin{subarray}{c}i,j=1\\\
i\neq j\end{subarray}}^{\eta}\left(\sum_{p,q\in G}\lvert p+\nu\rangle\langle
p\rvert_{i}\cdot\lvert q-\nu\rangle\langle q\rvert_{j}\right)\right\rVert$
$\displaystyle\leq\sum_{\nu\in
G_{0}}\frac{\eta(\eta-1)}{2\pi\Omega^{1/3}}\left|\frac{1}{\left\lVert\nu\right\rVert^{2}}-\frac{1}{\left\lVert\nu^{\prime}\right\rVert^{2}}\right|$
$\displaystyle=\frac{\eta(\eta-1)}{2\pi\Omega^{1/3}}\sum_{\mu=2}^{n_{p}+1}\sum_{\nu\in
B_{\mu}}\left|\frac{1}{\left\lVert\nu^{\prime}\right\rVert^{2}}-\frac{1}{\left\lVert\nu\right\rVert^{2}}\right|.$
(111)

Thus the total error in the potential component of the Hamiltonian can be
upper bounded as

$\left\lVert
U+V-(\widetilde{U}+\widetilde{V})\right\rVert\leq\frac{\eta}{2\pi\Omega^{1/3}}\left(\eta-1+2\lambda_{\zeta}\right)\sum_{\mu=2}^{n_{\mu}}\sum_{\nu\in
B_{\mu}}\left|\frac{1}{\left\lVert\nu^{\prime}\right\rVert^{2}}-\frac{1}{\left\lVert\nu\right\rVert^{2}}\right|.$
(112)

For given parameters, it is possible to bound the error by numerically
performing this sum. To give an analytic result, we take $\alpha=1$ and place
an upper bound on the sum. We obtain

$\displaystyle\sum_{\mu=2}^{n_{\mu}}\sum_{\nu\in
B_{\mu}}\left|\frac{1}{\left\lVert\nu^{\prime}\right\rVert^{2}}-\frac{1}{\left\lVert\nu\right\rVert^{2}}\right|$
$\displaystyle\leq\sum_{\mu=2}^{n_{p}+1}\sum_{\nu\in
B_{\mu}}\frac{16}{\mathcal{M}2^{2\mu}}$
$\displaystyle=\sum_{\mu=2}^{n_{p}+1}\left[(2^{\mu}-1)^{3}-(2^{\mu-1}-1)^{3}\right]\frac{16}{\mathcal{M}2^{2\mu}}$
$\displaystyle=\frac{4}{\mathcal{M}}\left(7\times
2^{n_{p}+1}-9n_{p}-11-3\times 2^{-n_{p}}\right).$ (113)

Thus the total error in the approximation of the Hamiltonian due to finite
$\mathcal{M}$ can be upper bounded by

$\left\lVert
U+V-(\widetilde{U}+\widetilde{V})\right\rVert\leq\frac{2\eta}{\pi\mathcal{M}\Omega^{1/3}}\left(\eta-1+2\sum_{\ell=1}^{L}\zeta_{\ell}\right)\left(7\times
2^{n_{p}+1}-9n_{p}-11-3\times 2^{-n_{p}}\right).$ (114)

This expression then provides the upper bound on the difference in eigenvalues
as required. ∎

As a result, if we are allowing an error due to the state preparation of
$\epsilon_{\mathcal{M}}$, then the value of $\mathcal{M}$ that should be
chosen is then about

$\mathcal{M}\approx\frac{28\eta
2^{n_{p}}}{\epsilon_{\mathcal{M}}\pi\Omega^{1/3}}\left(\eta+2\lambda_{\zeta}\right).$
(115)

Note that $\mathcal{M}$ is required to be a power of 2, so
$2^{n_{\mathcal{M}}}$. In numerical testing, we find that computing the actual
sum gives a result about $50\%$ of the upper bound. This is what would be
expected, because the ceiling will add $1/2$ on average, which is $50\%$ of
the upper bound of 1. If $\alpha\neq 1$ is used, then the error can be reduced
further to about $30\%$ of the upper bound, using $\alpha$ in the range
$1-3/2\mathcal{M}$ to $1-1/\mathcal{M}$.

The errors due to approximation of $\nu$ and approximation of $R_{\ell}$ can
be upper bounded by the sum of the two sources of error. Calling the $U$ with
both approximations $\overline{U}$, and that with just $R_{\ell}$ approximated
$\widetilde{U}$, we have
$\|\overline{U}-U\|\leq\|\overline{U}-\widetilde{U}\|+\|\widetilde{U}-U\|$.
Here $\|\widetilde{U}-U\|$ is upper bounded as in Lemma 2. For
$\|\overline{U}-\widetilde{U}\|$ exactly the same reasoning as in (II.5)
holds, except with $R_{\ell}$ replaced with $\widetilde{R}_{\ell}$. This means
that the same error bound holds and so the bounds on the error due to the two
sources can be added.

The other main source of error is the error in the phase estimation. With $m$
steps of phase estimation, the root-mean-square error in the estimated energy
due to the phase estimation is $\pi\lambda/(2m)$. It is possible to use $m$
that is not a power of 2 simply by using unary iteration on the control
register for the phase estimation. The effective value of $\lambda$ will vary
dependent on the relative sizes of the $T$, $U$, and $V$ operators.

As discussed above, if
$p_{\nu}\lambda_{T}=(1-p_{\nu})(\lambda_{U}+\lambda_{V})$, then in the case of
the failure of the inequality test for the preparation of $\nu$, we apply the
$\mathrm{SEL}$ operation for $T$, and the net value of $\lambda$ is
$\lambda_{T}+\lambda_{U}+\lambda_{V}$. If
$p_{\nu}\lambda_{T}>(1-p_{\nu})(\lambda_{U}+\lambda_{V})$, we apply
$\mathrm{SEL}$ for $T$ in the case of failure of the inequality test, _or_
$\lvert 0\rangle$ on an ancilla qubit, and again
$\lambda=\lambda_{T}+\lambda_{U}+\lambda_{V}$. In the case where
$p_{\nu}\lambda_{T}<(1-p_{\nu})(\lambda_{U}+\lambda_{V})$, we apply
$\mathrm{SEL}$ for $T$ in the case of failure of the inequality test, _and_
$\lvert 0\rangle$ on an ancilla qubit. In that case, the value of $\lambda$ is
$(\lambda_{U}+\lambda_{V})/p_{\nu}$. To account for the three cases, we have

$\lambda=\max\left[\lambda_{T}+\lambda_{U}+\lambda_{V},(\lambda_{U}+\lambda_{V})/p_{\nu}\right].$
(116)

We can alternatively perform amplitude amplification in the state preparation
for $\lvert\nu\rangle$, which triples that cost, but boosts the probability of
success to

$p_{\nu}^{\rm amp}=\sin^{2}\left(3\arcsin\left(\sqrt{p_{\nu}}\right)\right),$
(117)

which is close to 1. As a result the value of $\lambda$ can be given as

$\lambda=\max\left[\lambda_{T}+\lambda_{U}+\lambda_{V},(\lambda_{U}+\lambda_{V})/p_{\nu}^{\rm
amp}\right].$ (118)

A further subtlety is that we are testing $i\neq j$ in preparing the equal
superposition over $i$ and $j$ for $V$. In that case, if we apply $T$ in the
case where this test fails, then the effective $\lambda$ becomes

$\lambda=\max\left[\lambda_{T}+\lambda_{U}+\lambda_{V},[\lambda_{U}+\lambda_{V}/(1-1/\eta)]/p_{\nu}\right].$
(119)

The reason for this expression is that the preparation over all values of $i$
and $j$ gives an effective $\lambda_{V}$ which is proportional to $\eta^{2}$
rather than $\eta(\eta-1)$. This means the effective $\lambda_{V}$ is
$\lambda_{V}\eta^{2}/[\eta(\eta-1)]=\lambda_{V}/(1-1/\eta)$. Then, in the case
of amplitude amplification on $\lvert\nu\rangle$, the effective $\lambda$
would be

$\lambda=\max\left[\lambda_{T}+\lambda_{U}+\lambda_{V},[\lambda_{U}+\lambda_{V}/(1-1/\eta)]/p_{\nu}^{\rm
amp}\right].$ (120)

See Appendix D for a more in-depth derivation of this expression.

Another adjustment to the value of $\lambda$ is that needed to account for the
probabilities of failure in the preparation of the equal superposition states
for $w$, $i$ and $j$. These equal superpositions are required for both parts
of the Hamiltonian (kinetic and potential), and so the value of $\lambda$ is
divided by $P_{\rm eq}={\rm Ps}\left(3,8\right){\rm
Ps}\left(\eta+2\lambda_{\zeta},b_{r}\right){\rm
Ps}\left(\eta,b_{r}\right)^{2}$. That gives the resulting values of $\lambda$
without amplitude amplification

$\lambda=\max\left[\lambda_{T}+\lambda_{U}+\lambda_{V},[\lambda_{U}+\lambda_{V}/(1-1/\eta)]/p_{\nu}\right]/P_{\rm
eq},$ (121)

and with amplitude amplification

$\lambda=\max\left[\lambda_{T}+\lambda_{U}+\lambda_{V},[\lambda_{U}+\lambda_{V}/(1-1/\eta)]/p_{\nu}^{\rm
amp}\right]/P_{\rm eq}.$ (122)

Lastly we note that the $\lambda$-values are modified by our block encoding
method. We use $\lambda^{\prime}_{T}$ as given by (71), because of the state
preparation used for $T$. Moreover, because we use an imperfect state
preparation for $\nu$, the values of $\lambda_{U},\lambda_{V}$ are modified.
We define the modified $\lambda_{\nu}$ value

$\lambda_{\nu}^{\alpha}:=\alpha\sum_{\nu\in
G_{0}}\frac{\lceil{\mathcal{M}(2^{\mu-2}/\left\lVert\nu\right\rVert)^{2}}\rceil}{\mathcal{M}2^{2\mu-4}}.$
(123)

Then the $\lambda$-values for $U$ and $V$ are scaled in this way to give

$\lambda_{U}^{\alpha}:=\lambda_{\nu}^{\alpha}\lambda_{U}/\lambda_{\nu},\qquad\lambda_{V}^{\alpha}:=\lambda_{\nu}^{\alpha}\lambda_{V}/\lambda_{\nu}.$
(124)

As a result of these approaches, we have the following theorem for the
complexity.

######  Theorem 4.

The eigenvalue of the Hamiltonian given by (8) can be estimated to within RMS
error $\epsilon$ using a number of Toffoli gates

$\displaystyle\left\lceil\frac{\pi\lambda}{2\epsilon_{\rm
pha}}\right\rceil\big{\\{}2(n_{T}+4n_{\eta\zeta}+2b_{r}-12)+14n_{\eta}+8b_{r}-36+a[3n_{p}^{2}+15n_{p}-7+4n_{\mathcal{M}}(n_{p}+1)]$
$\displaystyle\quad+\lambda_{\zeta}+{\rm
Er}\left(\lambda_{\zeta}\right)+2(2n_{p}+2b_{r}-7)+12\eta
n_{p}+5(n_{p}-1)+2+24n_{p}+6n_{p}n_{R}+18$
$\displaystyle\quad+n_{\eta\zeta}+2n_{\eta}+6n_{p}+n_{\mathcal{M}}+16+\mathcal{O}(\log(1/\epsilon))\big{\\}}$
(125)

where “$a$” is equal to 3 if we perform amplitude amplification and 1
otherwise. For this complexity we require that a phase gradient state and a
$\lvert T\rangle$ state are given as a resource. The value of $\lambda$ is
given by

$\lambda=\max\left[\lambda^{\prime}_{T}+\lambda^{1}_{U}+\lambda^{1}_{V},[\lambda^{1}_{U}+\lambda^{1}_{V}/(1-1/\eta)]/p_{\nu}^{\rm
amp}\right]/P_{\rm eq}$ (126)

if we perform amplitude amplification, and

$\lambda=\max\left[\lambda^{\prime}_{T}+\lambda^{1}_{U}+\lambda^{1}_{V},[\lambda^{1}_{U}+\lambda^{1}_{V}/(1-1/\eta)]/p_{\nu}\right]/P_{\rm
eq}$ (127)

otherwise, where

$\displaystyle p_{\nu}$ $\displaystyle=\sum_{\mu=2}^{n_{p}+1}\sum_{\nu\in
B_{\mu}}\frac{\lceil{\mathcal{M}(2^{\mu-2}/\left\lVert\nu\right\rVert)^{2}}\rceil}{\mathcal{M}2^{2\mu}2^{n_{\mu}+1}},$
(128) $\displaystyle p_{\nu}^{\rm amp}$
$\displaystyle=\sin^{2}\left(3\arcsin\left(\sqrt{p_{\nu}}\right)\right),$
(129) $\displaystyle P_{\rm eq}$ $\displaystyle={\rm Ps}\left(3,8\right){\rm
Ps}\left(\eta+2\lambda_{\zeta},b_{r}\right){\rm
Ps}\left(\eta,b_{r}\right)^{2}.$ (130)

The value of $\epsilon_{\rm pha}$ is the allowable RMS error in phase
estimation. We require that $\epsilon_{\rm pha}>0$,
$n_{R},n_{\mathcal{M}},n_{T}\in\mathbb{N}$ are chosen such that

$\displaystyle\epsilon^{2}$ $\displaystyle\geq\epsilon_{\rm
pha}^{2}+(\epsilon_{\mathcal{M}}+\epsilon_{R}+\epsilon_{T})^{2}$ (131)
$\displaystyle\epsilon_{\mathcal{M}}$
$\displaystyle=\frac{2\eta}{2^{n_{\mathcal{M}}}\pi\Omega^{1/3}}\left(\eta-1+2\lambda_{\zeta}\right)\left(7\times
2^{n_{p}+1}-9n_{p}-11-3\times 2^{-n_{p}}\right)$ (132)
$\displaystyle\epsilon_{R}$
$\displaystyle=\frac{\eta\lambda_{\zeta}}{2^{n_{R}}\Omega^{1/3}}\sum_{\nu\in
G_{0}}\frac{1}{\left\lVert\nu\right\rVert}$ (133) $\displaystyle\epsilon_{T}$
$\displaystyle=\frac{\pi\lambda}{2^{n_{T}}}.$ (134)

###### Proof.

The complexity given is the sum of the costs for the block encoding given in
the list in Table 2, plus the number of qubits that need to be reflected upon
given in (98). The term $18$ is to account for the Toffolis used to produce
the qubits selecting between $T,U,V$, and the one flagging success of the
state preparation. That takes 9 Toffolis both in the preparation and the
inverse preparation. The term $\mathcal{O}(\log(1/\epsilon)$ is to account for
complexity of preparing the initial state to use as control for the phase
measurement, and the processing of this control at the end to obtain the phase
estimate.

The appropriate values of $\lambda$ in the cases with and without amplitude
amplification were given above in (122) and (121). In those expressions we
replace $\lambda_{T},\lambda_{U},\lambda_{V}$ with
$\lambda^{\prime}_{T},\lambda^{1}_{U},\lambda^{1}_{V}$ to account for our
block encoding, where we have taken $\alpha=1$ in order to use Lemma 3. The
error due to the phase estimation is independent from all other sources of
error, so we can upper bound the RMS error due to all sources as the sum of
squares of the error due to the phase estimation and the error from the other
sources. The bound $\epsilon_{\mathcal{M}}$ on the error due to finite
$\mathcal{M}$ is obtained from (115), which is based on Lemma 3. The bound
$\epsilon_{R}$ on the error due to the nuclear positions is obtained from
Lemma 2. In that expression we have substituted
$\delta_{R}=\Omega^{1/3}/2^{n_{R}+1}$, which comes from the cell length of a
side being $\Omega^{1/3}$, and the maximum error in the position being half of
$1/2^{n_{R}}$ times the length of the region. The bound $\epsilon_{T}$ on the
error due to the rotation of the qubit selecting between $T$ and $U+V$ is
given in (55).

The phase gradient state and $\lvert T\rangle$ state are used catalytically
throughout the procedure, with the phase gradient state being used to
implement phase rotations and the $\lvert T\rangle$ state used to implement
controlled Hadamards. These states can be prepared at the beginning of the
algorithm with negligible complexity compared to the rest of the algorithm. ∎

In order to choose parameters for the algorithm, we first take some maximum
allowable values for $\epsilon_{\mathcal{M}}$, $\epsilon_{R}$ and
$\epsilon_{T}$ as, for example, $\epsilon/10$. Then we substitute those values
into (132), (133), and (55) to find minimum values of $n_{\mathcal{M}}$,
$n_{R}$ and $n_{T}$. We can then use those values on the right-hand side (RHS)
of (132) and (133) to compute updated (smaller) values of
$\epsilon_{\mathcal{M}}$ and $\epsilon_{R}$. Those values can then be
substituted into (131) to solve for the maximum $\epsilon_{\rm pha}$ that
satisfies the inequality. These values of the parameters can then be
substituted into (4) to determine the total Toffoli cost. In order to optimise
the parameters, instead of using just the values of $n_{\mathcal{M}}$, $n_{R}$
and $n_{T}$ as given by this method, we can compute the cost for an array of
values of $n_{\mathcal{M}}$, $n_{R}$ and $n_{T}$, for example testing 4 below
and 4 above the numbers of bits given by this approach. We can then select
from this array the parameters that give the minimum Toffoli complexity.

Another improvement is to compute the exact bound on the state preparation
error, rather than just using the bound given in Lemma 3. That is, we use

$\epsilon_{\mathcal{M}}=\frac{\eta}{2\pi\Omega^{1/3}}\left(\eta-1+2\lambda_{\zeta}\right)\sum_{\mu=2}^{n_{\mu}}\sum_{\nu\in
B_{\mu}}\left|\frac{1}{\left\lVert\nu^{\prime}\right\rVert^{2}}-\frac{1}{\left\lVert\nu\right\rVert^{2}}\right|,$
(135)

with

$\frac{1}{\left\lVert\nu^{\prime}\right\rVert^{2}}=\frac{16\alpha}{\mathcal{M}2^{2\mu}}\left\lceil{\frac{\mathcal{M}2^{2\mu}}{16\left\lVert\nu\right\rVert^{2}}}\right\rceil.$
(136)

Choosing a finite number of bits for $n_{T}$ corresponds to giving a value of
$\alpha\neq 1$. Since the lowest $\epsilon_{\mathcal{M}}$ is obtained in the
region $\alpha=1-3/2\mathcal{M}$ to $1-1/\mathcal{M}$, we can choose $n_{T}$
in order to give $\alpha$ in that range and recompute
$\epsilon_{\mathcal{M}}$. The value of $n_{T}$ needed to obtain $\alpha$ in
that range is numerically found to be typically around $n_{\mathcal{M}}$ plus
4 to 8. With this choice of $n_{T}$ and $\alpha$, $\epsilon_{\mathcal{M}}$ is
about 30% of the bound proven in Lemma 3, and there is no additional error due
to finite $n_{T}$.

##  III The interaction picture based algorithm

In the interaction picture approach, one has a time-independent Hamiltonian
$H=A+B$, where the norm of $B$ is much smaller than that of $A$. The time
evolution under $A+B$ for time $\tau$ can then be given as (with $\mathcal{T}$
the time ordering operator)

$\displaystyle e^{-\mathrm{i}(A+B)\tau}$ $\displaystyle=e^{-\mathrm{i}\tau
A}\mathcal{T}\exp\left(-\mathrm{i}\int_{0}^{\tau}ds\,e^{\mathrm{i}sA}Be^{-\mathrm{i}sA}\right)$
(137)
$\displaystyle=\sum_{k=0}^{\infty}\frac{(-\mathrm{i})^{k}}{k!}\int_{0}^{\tau}d\tau_{1}\int_{0}^{\tau}d\tau_{2}\cdots\int_{0}^{\tau}d\tau_{k}\,e^{-\mathrm{i}(\tau-\tau^{\prime}_{k})A}Be^{-\mathrm{i}(\tau^{\prime}_{k}-\tau^{\prime}_{k-1})A}B\ldots
Be^{-\mathrm{i}(\tau^{\prime}_{2}-\tau^{\prime}_{1})A}Be^{-\mathrm{i}\tau_{1}A}$
$\displaystyle=\lim_{\begin{subarray}{c}K\rightarrow\infty\\\
M\rightarrow\infty\end{subarray}}\sum_{k=0}^{K}\frac{(-\mathrm{i}\tau)^{k}}{M^{k}k!}\sum_{m_{1}=0}^{M-1}\sum_{m_{2}=0}^{M-1}\cdots\sum_{m_{k}=0}^{M-1}e^{-\mathrm{i}\tau(M-1/2-m^{\prime}_{k})A/M}Be^{-\mathrm{i}\tau(m^{\prime}_{k}-m^{\prime}_{k-1})A/M}B\ldots$
$\displaystyle\quad\ldots
Be^{-\mathrm{i}\tau(m^{\prime}_{2}-m^{\prime}_{1})A/M}Be^{-\mathrm{i}\tau(m^{\prime}_{1}+1/2)A/M}$

where $\tau^{\prime}_{1},\ldots,\tau^{\prime}_{k}$ are sorted times from
$\tau_{1},\ldots,\tau_{k}$, and $m^{\prime}_{1},\ldots,m^{\prime}_{k}$ are
sorted integers from $m_{1},\ldots,m_{k}$. The approximation comes from both
truncating the sum at $K$ and discretising the integrals. The discretization
of the integral here is taken by using $\tau_{j}=(m_{j}+1/2)/M$, so the
midpoints of the intervals are used. This does not affect the final
expression, except that there are shifts by 1/2 on the beginning and ending
rotations.

Note that we do not immediately get a linear combination of unitaries from the
above Dyson-series expansion, as the operator $B$ from the Hamiltonian is not
unitary in general. Instead, we now describe the Hamiltonian input model for
the interaction picture algorithm in which $B$ is specified. We assume that
$B$ is given by a block encoding as

$\langle
0\rvert\mathrm{PREP}_{B}^{\dagger}\cdot\mathrm{SEL}_{B}\cdot\mathrm{PREP}_{B}\lvert
0\rangle=\frac{B}{\lambda_{B}}\,.$ (138)

Here, $\mathrm{PREP}_{B}$ can be viewed as some state preparation procedure
acting nontrivially only on an ancilla register initialized in state $\lvert
0\rangle$, and $\mathrm{SEL}_{B}$ is a selection unitary acting jointly on the
ancilla register and the register of the target system. By preparing the
ancilla state, performing the selection operation, and unpreparing the state,
we can implement operator $B$ on the target system with normalization factor
$\lambda_{B}\geq\left\lVert B\right\rVert$. We describe a quantum circuit in
Section III.3 that block-encodes $B$ for the first quantized quantum chemistry
simulation.

For the operator $A$, we assume that the matrix exponential
$e^{-\mathrm{i}sA}$ can be efficiently implemented for an arbitrary real
number $s$. This is the case for instance when $A$ is already diagonalized and
contains pairwise commuting terms, each of which can be directly
exponentiated. In particular, this holds for the kinetic operator $T$ when the
first quantized Hamiltonian is represented in the plane wave basis. See
Section III.2 for a detailed discussion of the circuit implementation.

Using the block encoding of $B$, we now rewrite the Dyson-series expansion as

$\displaystyle e^{-\mathrm{i}(A+B)\tau}$
$\displaystyle\approx\sum_{k=0}^{K}\frac{(-\mathrm{i}\lambda_{B}\tau)^{k}}{M^{k}k!}\langle
0\rvert^{\otimes k}\mathrm{PREP}_{B}^{\dagger\otimes
k}\sum_{m_{1},...,m_{k}=0}^{M-1}\bigg{(}e^{-\mathrm{i}\tau(M-1/2-m^{\prime}_{k})A/M}\mathrm{SEL}_{B}e^{-\mathrm{i}\tau(m^{\prime}_{k}-m^{\prime}_{k-1})A/M}\mathrm{SEL}_{B}\ldots$
$\displaystyle\quad\ldots\mathrm{SEL}_{B}e^{-\mathrm{i}\tau(m^{\prime}_{2}-m^{\prime}_{1})A/M}\mathrm{SEL}_{B}e^{-\mathrm{i}\tau(m^{\prime}_{1}+1/2)A/M}\bigg{)}\mathrm{PREP}_{B}^{\otimes
k}\lvert 0\rangle^{\otimes k}$ (139) $\displaystyle=\left(\langle
0\rvert\mathrm{PREP}_{B}^{\dagger}\right)^{\otimes
K}\\!\sum_{k=0}^{K}\frac{(-\mathrm{i}\lambda_{B}\tau)^{k}}{M^{k}k!}\\!\\!\\!\sum_{m_{1},...,m_{k}=0}^{M-1}\\!\bigg{(}e^{-\mathrm{i}\tau(M-1/2-m^{\prime}_{k})A/M}\mathrm{SEL}_{B}e^{-\mathrm{i}\tau(m^{\prime}_{k}-m^{\prime}_{k-1})A/M}\mathrm{SEL}_{B}\ldots$
$\displaystyle\quad\ldots\mathrm{SEL}_{B}e^{-\mathrm{i}\tau(m^{\prime}_{2}-m^{\prime}_{1})A/M}\mathrm{SEL}_{B}e^{-\mathrm{i}\tau(m^{\prime}_{1}+1/2)A/M}\bigg{)}\left(\mathrm{PREP}_{B}\lvert
0\rangle\right)^{\otimes K}.$

Here for each fixed value of $k$, the $k$ selection operators
$\mathrm{SEL}_{B}$ act on $k$ different ancilla registers and the same target
register; the last equality follows from the fact that the remaining $K-k$
ancilla states $\mathrm{PREP}_{B}\lvert 0\rangle$ cancel with $\langle
0\rvert\mathrm{PREP}_{B}^{\dagger}$ in pairs and so the number of their
occurrences can be made independent of $k$. Note that the resulting expression
is expressed as a linear combination of unitaries and can thus in principle be
implemented using quantum circuits.

Although we have written this as performing all the preparations at the
beginning and inverting them at the end, in practice we will perform the
preparations and inverse preparations for each $\mathrm{SEL}_{B}$ before the
next preparation. That avoids needing to use all qubits for all preparations
at once, which would have a much higher qubit cost. A further simplification
is that we can separate out $e^{-\mathrm{i}\tau A/2M}$ at the beginning and
the end. Then $e^{-\mathrm{i}\tau(M-1/2-m^{\prime}_{k})A/M}$ is replaced with
$e^{-\mathrm{i}\tau(M-1-m^{\prime}_{k})A/M}$. The advantage here is that
$M-1-m^{\prime}_{k}$ can be computed simply by performing NOT gates on all
qubits, with zero Toffoli cost.

To implement this linear combination of unitaries, we need to prepare quantum
states encoding the coefficients from the combination. For readability, we
describe a naive circuit implementation here, and leave a detailed discussion
of the subtleties and improvements in subsequent subsections. We first prepare

$\frac{1}{\sqrt{\beta}}\sum_{k=0}^{K}\sqrt{\frac{(-\mathrm{i}\lambda_{B}\tau)^{k}}{k!}}\lvert
k\rangle,$ (140)

where $\beta$ is the normalization factor and $\lvert k\rangle$ is encoded in
unary as $\lvert 1^{k}0^{K-k}\rangle$. This state may be prepared using the
naive approach with one rotation gate and $K$ controlled rotation gates.
However, we propose an improved approach based on the inequality test which
can significantly reduce the implementation cost; see Section III.4 for
details.

Next, we prepare the correct time indices used in the interaction picture
simulation. One possible method is to prepare an equal superposition over $k$
time indices $m_{1},...,m_{k}$ and quantumly sort them Kieferová  _et al._
(2019), giving

$\frac{1}{\sqrt{\beta}}\sum_{k=0}^{K}\sqrt{\frac{(-\mathrm{i}\lambda_{B}\tau)^{k}}{M^{k}k!}}\lvert
k\rangle\sum_{m_{1},...,m_{k}=0}^{M-1}\lvert
m_{1}^{\prime},...,m_{k}^{\prime}\rangle,$ (141)

where again $0\leq m_{1}^{\prime}\leq...\leq m_{k}^{\prime}\leq M-1$ are
sorted integers from $m_{1},..,m_{k}$. We then compute the time differences as

$\frac{1}{\sqrt{\beta}}\sum_{k=0}^{K}\sqrt{\frac{(-\mathrm{i}\lambda_{B}\tau)^{k}}{M^{k}k!}}\lvert
k\rangle\sum_{m_{1},...,m_{k}=0}^{M-1}\lvert
m_{1}^{\prime},m_{2}^{\prime}-m_{1}^{\prime},...,m_{k}^{\prime}-m_{k-1}^{\prime},M-1-m_{k}^{\prime}\rangle.$
(142)

These difference values are further used to control the matrix exponentials in
the Dyson-series expansion. Specifically, we have the decomposition

$\sum_{m=0}^{M-1}\lvert m\rangle\langle m\rvert\otimes e^{-\mathrm{i}\tau
mA/M}=\prod_{j=0}^{\log(M)}\left(\lvert 0\rangle\langle 0\rvert_{j}\otimes
I+\lvert 1\rangle\langle 1\rvert_{j}\otimes e^{-\mathrm{i}\tau
2^{j}A/M}\right),$ (143)

where each controlled exponential of $A$ can be directly implemented. In
practice, in the implementation it will convenient to apply the controlled
exponential jointly controlled by all bits of $m$, rather than controlling on
the bits of $m$ independently. On the other hand, we apply the selection
unitaries $\mathrm{SEL}_{B}$ controlled by the unary-encoded $k$ values. These
operations, together with the preparation and unpreparation procedures,
implement the Dyson-series expansion up to the normalization factor $\beta$.
Beyond the sorting method, previous work also suggested alternative methods
based on a “compressed encoding” Kieferová  _et al._ (2019) and “inequality
tests” Low and Wiebe (2018). We compare these methods and justify our choice
of implementation in Section III.1.

For $K$ sufficiently large, the normalization factor $\beta$ is approximately
given by

$\beta=\sum_{k=0}^{K}\frac{(\lambda_{B}\tau)^{k}}{k!}\approx\exp(\lambda_{B}\tau).$
(144)

In the original interaction picture algorithm, we choose
$\tau\approx\ln(2)/\lambda_{B}$ so that $\beta\approx 1/2$. Then, a single
step of oblivious amplitude amplification boosts the success amplitude close
to $1$. The ideal evolution can be simulated for time $t$ by repeating this
$\lambda_{B}t/\ln(2)$ times. However, for performing phase estimation, we
develop a new approach by applying qubitization on the pre-amplified
interaction picture subroutine, which significantly reduces the implementation
cost. We discuss our approach in detail in Section III.1.

###  III.1 Implementation choices for the interaction picture algorithm

There are two main choices for how to provide the ordered integers
$m^{\prime}_{j}$. One is to generate a sequence of differences according to a
negative exponential weighting. Then one directly has the differences
$m^{\prime}_{j+1}-m^{\prime}_{j}$, but it is necessary to add these
differences together to determine which difference gives a total over $M$.
That is used to determine $k$, rather than generating a register for $k$
directly. This method is called the “compressed encoding” in Kieferová  _et
al._ (2019). The preparation of the state with the negative exponential
weighting can be performed as described in Berry  _et al._ (2014a) with a
sequence of controlled rotations. That results in state preparation complexity
corresponding to a factor of $\log(1/\epsilon)$ times the number of qubits in
the time registers.

The other choice for preparing the state with the ordered integers
$m^{\prime}_{j}$ is to prepare a set of unordered integers, then sort them.
The type of sort needed is a quantum version of a sorting network, as
described in Sheng-Tzong Cheng and Chun-Yen Wang (2006); Beals  _et al._
(2013). There are algorithms for generating sorting networks, but in practice
for small numbers (as we have here) it is best to choose optimized sorting
networks. The number of comparator operations required for the best known
sorting networks for up to 16 items are as shown in Table 3 Codish  _et al._
(2014). Each comparator needed corresponds to an inequality test and a
controlled swap. The inequality test uses approximately a Toffoli for each
qubit of the registers, and the controlled swap uses a Toffoli for each qubit
of the registers. The number of Toffolis therefore corresponds to twice the
number of comparators multiplied by the quits used for each time register. For
the case of $n=16$, that is $7.5$ times the number of qubits in the time
registers, as compared to $\log(1/\epsilon)$ for the compressed encoding
approach. This will be a smaller multiplying factor, which is why we will use
the sorting approach instead of the compressed encoding.

$K$ | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9 | 10 | 11 | 12 | 13 | 14 | 15 | 16
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---
comparators | 1 | 3 | 5 | 9 | 12 | 16 | 19 | 25 | 29 | 35 | 39 | 45 | 51 | 56 | 60
Table 3: The number of comparators needed for a quantum sort of $K$ items.
Each comparator involves an inequality test and a controlled swap. The number
of comparators for a given $K$ is denoted ${\rm St}\left(K\right)$.

For the implementation of $e^{-\mathrm{i}(A+B)\tau}$ as a linear combination
of unitaries, the general method is to prepare a superposition state with
weights that are the square roots of the amplitudes in the linear combination
of unitaries, then apply controlled unitary operations. Oblivious amplitude
amplification is used in order to implement the step deterministically, and in
order to only need a single step of oblivious amplitude amplification it is
necessary to choose the parameters so the amplitude of success is $1/2$.

In the case where the sorting is performed, the sum of the amplitudes needed
for the approximation of $e^{-\mathrm{i}(A+B)\tau}$ as a linear combination of
unitaries is

$\sum_{k=0}^{K}\frac{\tau^{k}}{k!}\lambda_{B}^{k}\leq\exp(\tau\lambda_{B}).$
(145)

An inequality is written here due to the truncation at finite $K$, but the two
sides of this equation will be very close (with error corresponding to error
in the evolution), so it is sufficient to assume that they are equal for the
purpose of costing the algorithm. The amplitude for success in the linear
combination of unitaries is the inverse of this value, so in order to obtain
amplitude $1/2$, we need $\exp(\tau\lambda_{B})\approx 2$, which corresponds
to choosing $\tau\approx\ln 2/\lambda_{B}$. In order to perform an accurate
phase estimation we need long time evolution, so greater efficiency is
obtained by the time being broken into intervals $\tau$ that are as large as
possible.

An alternative approach is given by Low and Wiebe (2018), where instead of
performing a sort, inequality tests are performed in the time registers. When
performing the inequality test, the factor of $k!$ is removed in the
denominator, because that corresponds to the number of permutations of times.
As a result, the sum of amplitudes in the method of Low and Wiebe (2018) is

$\sum_{k=0}^{K}\tau^{k}\lambda_{B}^{k}\leq\frac{1}{1-\tau\lambda_{B}}.$ (146)

Again the two sides of the expression are approximately equal. For the
amplitude amplification, we need this amplitude to be $2$ (because the
amplitude for success is the inverse of the sum of amplitudes). That means
that in the approach of Low and Wiebe (2018) one needs $\tau\approx
1/(2\lambda_{B})$. This value of $\tau$ is about $72\%$ of the value of $\tau$
in the method of Kieferová  _et al._ (2019), which means that the overall
complexity is about $39\%$ larger. The method of Low and Wiebe (2018) saves
the complexity of the sort, but that complexity is trivial as compared to
other complexities in the algorithm.

So far we have been assuming that phase estimation is performed on the time
evolution, which requires performing amplitude amplification. That triples the
cost, because even a single step of oblivious amplitude amplification requires
a forward step, a reverse step, then a forward step again. An alternative
method can be obtained by simply performing phase estimation on the single
step of the evolution, as per the discussion in Berry  _et al._ (2018); Poulin
_et al._ (2017). The difference here is that we are block encoding the
evolution $e^{-\mathrm{i}(A+B)\tau}$ rather than the Hamiltonian, and the
evolution is not Hermitian so the analysis does not hold. However, we can
easily block encode
$\sin[(A+B)\tau]=\left[e^{\mathrm{i}(A+B)\tau}-e^{-\mathrm{i}(A+B)\tau}\right]/(2\mathrm{i})$.
All we need to do is use a single qubit that chooses between the block
encoding of $e^{\mathrm{i}(A+B)\tau}$ and $e^{-\mathrm{i}(A+B)\tau}$. All that
qubit needs to do is control the sign of the exponential in terms like
$e^{-\mathrm{i}\tau(m^{\prime}_{2}-m^{\prime}_{1})A/M}$, which can be achieved
with negligible complexity.

If the operation $\sin[(A+B)\tau]$ has eigenvalues $\mu_{j}=\sin(E_{j}\tau)$,
where $E_{j}$ are eigenvalues of $A+B$, then as discussed in Berry  _et al._
(2018) the qubitization procedure gives an operator with eigenvalues $\pm
e^{\pm\mathrm{i}\arcsin(\mu_{j}/\lambda)}$, where $\lambda$ is the sum of the
amplitudes in the block encoding of $\sin[(A+B)\tau]$. (In Berry  _et al._
(2018) $H$ is equivalent to $\sin[(A+B)\tau]$ here, and $E_{k}$ there is
equivalent to $\mu_{j}$ here.) Now the value of $\lambda$ in the block
encoding of $\sin[(A+B)\tau]$ is the same as that of
$e^{-\mathrm{i}(A+B)\tau}$, and so with the method of Kieferová  _et al._
(2019) is $\exp(\tau\lambda_{B})$. The magnitude of the eigenvalues in the
exponential are therefore $\arcsin[\sin(E_{j}\tau)/\exp(\tau\lambda_{B})]$. In
the case of small $E_{j}$ (as we expect here), the linearization is
approximately $E_{j}\tau/\exp(\tau\lambda_{B})$, so this is equivalent to
using an effective evolution time $\tau_{\rm eff}=\tau/\exp(\tau\lambda_{B})$.
If we were to choose $\tau\approx\ln 2/\lambda_{B}$, then we would have
$\tau_{\rm eff}\approx\ln 2/(2\lambda_{B})$. That is, the evolution time is
halved from $\tau$, which corresponds to doubling the cost for the phase
estimation. On the other hand, the oblivious amplitude amplification triples
the cost, so this qubitization method gives $2/3$ the cost for performing
phase estimation on the time evolution. However, we can further reduce the
cost by increasing $\tau$ to $1/\lambda_{B}$. That gives an effective time
$\tau_{\rm eff}\approx 1/(e\lambda_{B})$ where $e$ is Euler’s number. That
means the complexity is then $63\%$ of the method via time evolution from
Kieferová  _et al._ (2019), or $45\%$ of the complexity using the method of
Low and Wiebe (2018). Therefore, we see that even though these individual
improvements are only minor, there is total improvement by more than a factor
of 2 over the method of Low and Wiebe (2018).

For this effective time, we are assuming that $E_{j}$ is near zero, so
$\arcsin[\sin(E_{j}\tau)/e]$ (where we have taken $\tau=1/\lambda_{B}$) is
approximately $E_{j}\tau/e$. In cases where $E_{j}$ is not near zero, then the
uncertainty in estimating $\arcsin[\sin(E_{j}\tau)/e]$ will correspond to
uncertainty in estimating $E_{j}$ with a different constant factor. However,
we have freedom in adding a phase shift to the operator we are implementing.
Instead of implementing $e^{\pm\mathrm{i}(A+B)\tau}$, we can implement
$e^{\pm\mathrm{i}(A+B-E_{0})\tau}$, simply by performing an additional phase
shift. That means we are then estimating eigenvalues
$\arcsin[\sin((E_{j}-E_{0})\tau)/e]$. Provided we have an initial estimate of
$E_{j}$, we can ensure that $E_{j}-E_{0}$ is near zero so we are in the region
where the sine and arcsine can be linearized.

To be more specific, the slope at $E_{j}-E_{0}$ is increased to

$\frac{\tau}{e}\left(1+\frac{(e^{2}-1)(E_{j}-E_{0})^{2}\tau^{2}}{2e^{2}}+\mathcal{O}((E_{j}-E_{0})^{4}\tau^{4})\right).$
(147)

This means that there is a relative increase in the number of steps needed to
obtain a given accuracy of $\mathcal{O}((E_{j}-E_{0})^{2}\tau^{2})$. There is
also a slight increase in the estimation error due to the nonlinearity of this
function. This relative increase turns out to be on the order of the standard
deviation of the phase estimation, so scales as $\mathcal{O}(1/\mathcal{N})$
for phase estimation on $\mathcal{N}$ steps. This higher-order increase in the
error can be accounted for by increasing $\mathcal{N}$ by $\mathcal{O}(1)$. We
do not need to take account of these factors when considering phase estimation
for the qubitization, because there the effect of the nonlinearity is to
reduce the error.

To summarise, there are two main choices for the algorithm. First, there is
the choice of how to prepare the time registers, for which there are three
main alternatives.

  1. 1.

The “compressed encoding” method from Kieferová  _et al._ (2019), which
requires controlled rotations by arbitrary angles on the qubits for the time
registers, and therefore has a complexity scaling like the logarithm of the
inverse error for each qubit of the time registers.

  2. 2.

The sorting method from Kieferová  _et al._ (2019), which requires a number of
Toffolis for each qubit that can be determined from Table 3, and is expected
to be around 8 (so less than for compressed encoding).

  3. 3.

The method of Low and Wiebe (2018), which involves postselecting on success of
inequality tests, rather than a sort. That reduces the complexity of the state
preparation for the time registers, but means the number of steps is increased
by a factor of $2\ln 2\approx 1.38$ over the methods of Kieferová  _et al._
(2019). This is expected to result in an increased overall cost, because the
state preparation for the time registers is a small part of the overall cost
in each step.

The other choice is between the following two alternatives.

  1. 1.

One can simulate the time evolution under the Hamiltonian, and perform phase
estimation on the time evolution. Then the number of steps scales as
$\lambda_{B}T/\ln 2$ for the methods of Kieferová  _et al._ (2019), or
$2\lambda_{B}T$ for the method of Low and Wiebe (2018).

  2. 2.

One can instead perform a qubitization of $\sin(H\tau)$, and perform phase
estimation on that. That reduces the complexity because it removes the need
for amplitude amplification for time evolution (and the corresponding error).
The number of steps would be $e\lambda_{B}T$ for the methods of Kieferová  _et
al._ (2019) or $2e\ln 2\lambda_{B}T$ for the method of Low and Wiebe (2018).
The number of steps is increased, but the complexity of each step is reduced
by a factor of 3 as compared to the time evolution approach.

We will focus on the case where the time registers are prepared using a sort,
and we perform qubitization of $\sin(H\tau)$, as this gives the best
performance. As in the case of the qubitization approach, we need to ensure
that our $\mathrm{SEL}$ procedure is self-inverse. We will use an ancilla to
control $e^{-\mathrm{i}(A+B)\tau}$ versus $e^{\mathrm{i}(A+B)\tau}$. To ensure
that it is self-inverse, it will be controlled on the ancilla qubit as

$\lvert 0\rangle\langle 1\rvert\otimes e^{-\mathrm{i}(A+B)\tau}+\lvert
1\rangle\langle 0\rvert\otimes e^{\mathrm{i}(A+B)\tau}.$ (148)

This ancilla qubit is prepared as $\lvert+\rangle$, ensuring we obtain the sum
of $e^{-\mathrm{i}(A+B)\tau}$ and $e^{\mathrm{i}(A+B)\tau}$ in the block
encoding. However, we do need to ensure that the implementation of
$e^{-\mathrm{i}(A+B)\tau}$ is the inverse of $e^{\mathrm{i}(A+B)\tau}$. To
achieve that, we apply the following method.

  * •

Use the qubit as control for a phase gate on each qubit in the unary encoding
of $k$ to control between $(-\mathrm{i})^{k}$ and $(\mathrm{i})^{k}$ in the
Dyson series. There is no extra Toffoli cost.

  * •

Use the qubit to control whether the sort of the time registers sorts them in
ascending or descending order. One can simply use the qubit as a control with
the qubits flagging the result of the inequality test as target. There is
again no Toffoli cost.

  * •

Use the qubit to control taking the absolute values of the time differences
(see Appendix E). This is the only part with a Toffoli cost.

  * •

Use this qubit to control between adding and subtracting $\nu$, which has no
extra Toffoli cost.

  * •

Use the qubit to control positive versus negative phasing by the kinetic
energy, which can be performed without Toffolis by controlling between adding
and subtracting the phase into the phase gradient state.

###  III.2 Exponentiating the kinetic operator

It is estimated in Babbush  _et al._ (2019) that the cost of simulating the
first quantized electronic structure Hamiltonian is dominated by the
exponentiation $e^{-\mathrm{i}\tau T}$. The reason for that is the kinetic
energy needs to be computed in the computational basis, which requires
squaring $3\eta$ numbers (for the 3 components of the momenta of each
electron) and adding them together. Squaring an $n_{p}$-bit number can be
performed with approximately $n_{p}^{2}$ Toffolis, giving an overall
complexity of approximately $3\eta n_{p}^{2}$.

However, it is possible to greatly reduce the complexity, by keeping the sum
of squares in an ancilla, instead of recomputing it every time. Then when we
do the block encoding of the $U+V$, we can update this sum of squares. The
cost of updating the sum of squares is only logarithmically dependent on
$\eta$. To be more specific, for $V$ we have $p$ replaced with $p+\nu$ and $q$
replaced with $q-\nu$. For the first the change in the sum of squares is
$\|p+\nu\|^{2}-\|p\|^{2}=2p\cdot\nu+\|\nu\|^{2}$, and for the second the
difference in squares is $\|q-\nu\|^{2}-\|q\|^{2}=-2q\cdot\nu+\|\nu\|^{2}$.
The quantity $\|\nu\|^{2}$ has already been computed in preparing the state
with amplitudes $1/\|\nu\|$, so does not need to be recomputed.

This method also enables us to take account of the energy offset $E_{0}$ we
are using. We aim to subtract all energies by $E_{0}$ to put us in the linear
region of the function we are estimating. We can simply subtract $E_{0}$ from
the kinetic energy register at the beginning of the procedure. That means the
energy of every step will be shifted by $E_{0}$ with no more gate complexity
(other than the initial subtraction).

Therefore the majority of the cost is in the computation of $p\cdot\nu$ and
$q\cdot\nu$. The cost for computing the product of the _absolute value_ of
each component is $2n_{p}(n_{p}-1)-n_{p}=2n_{p}^{2}-3n_{p}$, and the sign of
the product can be computed with a CNOT. Rather than adding all three
components then adding into the energy register, it is convenient to add them
one at a time. This is because we have signed integers, and adding the three
signed integers would incur an additional cost for converting to two’s
complement. When adding the signed integer into the energy register, we can
instead use the sign bit to control addition or subtraction into the energy
register, which incurs no additional cost.

The number of bits needed to represent the sum of momenta is
$n_{\eta}+2n_{p}$. We are guaranteed that the net result of all the arithmetic
will result in no overflow qubit, but because we are adding terms in
succession there may be overflow in intermediate steps. Therefore the cost of
addition of $\|\nu\|$ will be $n_{\eta}+2n_{p}$. On the other hand, because we
have twice $p\cdot\nu$ and $q\cdot\nu$, the least significant bit is 0, and
the cost of the addition of each component is reduced by 1 to
$n_{\eta}+2n_{p}-1$.

We will have qubits that control whether $U$ and $V$ are implemented, and
these will need to control whether each of the additions are performed. For
$\|\nu\|^{2}$ we have a cost of the control each time of $2n_{p}$, for a total
of $4n_{p}$. To eliminate the cost of controlling the addition of $p\cdot\nu$
and $q\cdot\nu$, we can instead make the swaps of these momenta into the
ancilla controlled on the qubits which control whether $U$ and $V$ are
implemented. That control on the unary iteration on $i$ and $j$ takes only one
extra Toffoli each time, for a total cost of 4 Toffolis there. However, if the
swap was not performed, then $p\cdot\nu$ or $q\cdot\nu$ was computed with a
momentum of zero, giving zero. Therefore these additions do not need to be
made controlled.

There is an important nuance in that the phase needed will not be exactly a
multiple of $2\pi/2^{b_{\rm grad}}$ for some integer $b_{\rm grad}$. The
kinetic energy is given by $2\pi^{2}/\Omega^{2/3}$ times the sum of squares of
the integers stored in the momentum registers, and this energy is multiplied
by $\tau=1/(\lambda_{U}+\lambda_{V})$ to give the phase required. The ideal
situation is that $2\tau\pi^{2}/\Omega^{2/3}$ (the smallest phase shift) is
exactly $2\pi/2^{b_{\rm grad}}$. Then one can simply take the sum of squares
of momenta and add it into the phase gradient register.

On the other hand, say that for example $2\tau\pi^{2}/\Omega^{2/3}$ is instead
$3\times 2\pi/2^{b_{\rm grad}}$. One can then instead use the modified phase
gradient state

$\frac{1}{\sqrt{2^{b_{\rm grad}}}}\sum_{\ell=0}^{2^{b_{\rm
grad}}-1}e^{-6\pi\mathrm{i}\ell/2^{b_{\rm grad}}}\lvert\ell\rangle.$ (149)

Then it is easy to see that adding an integer to the register will give a
phase shift that is a multiple of $3\times 2\pi/2^{b_{\rm grad}}$. Typically,
$2\tau\pi^{2}/\Omega^{2/3}$ would not be of the form of an integer multiple of
$2\pi/2^{b_{\rm grad}}$, but we have freedom to adjust the time $\tau$. Say
for example that we are using a multiple with 10 bits. We will use $b_{T}$ for
the number of bits in the multiple. Then the worst case is when
$2\tau\pi^{2}/\Omega^{2/3}$ with $\tau=1/(\lambda_{U}+\lambda_{V})$ is about
$512.5$ times $2\pi/2^{b_{\rm grad}}$. Then rounding to 513 gives the largest
relative error, whereas rounding numbers like, for example, $599.5$ to $600$
would result in smaller relative error. (Numbers below 512 could be
represented with fewer bits.) Adjusting the value of $\tau$ to make the
multiple exactly 513, means that the value of $\tau/\exp(\tau\lambda_{B})$
(which governs the number of steps needed in the phase estimation) would be
changed by less than 1 part in a million. In general, as discussed in Appendix
F, one can take

$b_{\rm
grad}=b_{T}-\left\lceil\log\left(\frac{\pi}{(\lambda_{U}+\lambda_{V})\Omega^{2/3}}\right)\right\rceil,$
(150)

and the increase in the $\lambda$ will be no more than $1+1/2^{2b_{T}+1}$.

To summarise, the costs associated with the phasing by $T$ are as follows.

  1. 1.

We need to compute time differences from two times of $n_{t}$ qubits with no
carry bit, which has cost $n_{t}-1$.

  2. 2.

There is cost $2n_{t}(n_{\eta}+2n_{p})-n_{t}$ for the multiplication of the
sum of momenta by the difference of times.

  3. 3.

The addition into the phase gradient register has cost $b_{\rm grad}-2$
Toffolis.

  4. 4.

The cost of computing all the products of components of $p$ or $q$ with $\nu$
is $6\times(2n_{p}^{2}-3n_{p})=12n_{p}^{2}-18n_{p}$ altogether.

  5. 5.

The cost of the addition of all these components into the energy register is
$6(n_{\eta}+2n_{p}-1)$.

  6. 6.

The cost of two additions of $\|\nu\|^{2}$ into the energy is
$2(n_{\eta}+2n_{p})$.

  7. 7.

The cost of making each addition of $\|\nu\|^{2}$ controlled is $4n_{p}$.

  8. 8.

For the kinetic energy phasing for time $\tau/2M$ at the beginning and the end
we simply need to add the (appropriately bit-shifted) kinetic energy into the
phase gradient state, with cost $b_{\rm grad}-2$ Toffolis each time.

These costs are needed different numbers of times. For the time difference,
there is only Toffoli cost where we compute differences of times given in
quantum registers, which is done $K-1$ times. This is because
$M-1-m^{\prime}_{k}$ can be computed with NOT gates (see Appendix E for more
detail). The multiplication in step 2 and the phasing in step 3 are both
needed $K+1$ times for all of the phasings. The remaining costs are updating
the energy register, which only needs to be done between phasings, so we have
$K$ of these costs.

The arithmetic for the time differences and multiplications can be performed
in place with ancillae, so the arithmetic can be inverted without Toffolis.
This gives a complexity for each time difference $n_{t}-1$. There is also cost
$n_{t}-1$ to take the absolute value of the time, as controlled by the qubit
selecting between forward and reverse time evolution (see Appendix E). With
$K-1$ differences, the total cost is

$2(K-1)(n_{t}-1).$ (151)

For the phase shift the total Toffoli cost is

$(K+1)[2n_{t}(n_{\eta}+2n_{p})-n_{t}+b_{\rm grad}-2]+2(b_{\rm grad}-2),$ (152)

where we have included $2(b_{\rm grad}-2)$ for the initial and final phasings.
For updating the sum of momenta

$\displaystyle
K\\{[12n_{p}^{2}-18n_{p}]+[6(n_{\eta}+2n_{p}-1)]+[2(n_{\eta}+2n_{p})]+[4n_{p}]\\}$
$\displaystyle=K(12n_{p}^{2}+2n_{p}+8n_{\eta}).$ (153)

Note that the complexity is now only logarithmically dependent on $\eta$.
There is still linear dependence on $\eta$ for selecting momentum registers in
block encoding $U+V$.

###  III.3 Implementing $\mathrm{SEL}$ and $\mathrm{PREP}$

In this section, we discuss the costing of the block encoding of $U+V$ for the
interaction picture. As part of this block encoding, we will update the sum of
squares of momenta, as described above. We have previously described the
complexity for block encoding $T+U+V$. The complexity of block encoding $U+V$
is very similar, except the parts needed to block encode $T$ are omitted, and
we need to update the register storing the sum of momenta. Therefore our
previous costing of the block encoding from Table 2 can be modified as in
Table 4. We will make the controlled swaps of momenta into the ancillas
controlled on qubits flagging whether $U$ and $V$ are to be performed. That
overall control takes one extra Toffoli for each of the controlled swaps, for
a total of 4. The final item in this table is the new cost for updating the
sum of momenta.

Procedure |  Toffoli cost
---|---
Preparing the superposition on the qubit selecting between $U$ and $V$, but not $T$; see (58). |  $2(4n_{\eta\zeta}+2b_{r}-9)$
Preparing equal superpositions over $\eta$ values of $i$ and $j$ in unary; see (62). |  $14n_{\eta}+8b_{r}-36$
The state preparation cost for the $w,r,s$ registers used for $T$ is eliminated. |  $0$
Controlled swaps of the $p$ and $q$ registers into and out of ancillae (which is used for $U$ and $V$), with another 4 Toffolis to give overall control; see (72). |  $12\eta n_{p}+4\eta-4$
The $\mathrm{SEL}$ cost for $T$ is eliminated. |  $0$
The preparation of the state with amplitude $1/{\left\lVert\nu\right\rVert}$ and inversion; see (90). |  $3n_{p}^{2}+15n_{p}-7+4n_{\mathcal{M}}\big{(}n_{p}+1\big{)}$
The QROM for $R_{\ell}$, combined with the state preparation with amplitudes $\sqrt{\zeta_{\ell}}$ (which may be implicit); see (92). |  $\lambda_{\zeta}+{\rm Er}\left(\lambda_{\zeta}\right)$
The additions and subtractions of $\nu$ into the momentum registers; see (93). |  $24n_{p}$
Phasing by $-e^{-\mathrm{i}k_{\nu}\cdot R_{\ell}}$; see (97). |  $6n_{p}n_{R}$
Updating the sum of momenta; see (III.2). |  $12n_{p}^{2}+2n_{p}+8n_{\eta}$
Table 4: The costs involved in the block encoding of the Hamiltonian for the
interaction picture, following the same entries as in Table 2.

For the total complexity, in the Dyson series the most efficient method is to
qubitize a single step of the Dyson series, and perform phase estimation on
it, rather than to use amplitude amplification to obtain the explicit time
evolution. As a result, the effective time for a single step is approximately
$1/[e\lambda_{B}]$. The RMS error in the energy due to phase estimation is
then $\pi e\lambda_{B}/(2\mathcal{N})$, where $\mathcal{N}$ is the number of
steps, so one should take

$\mathcal{N}=\left\lceil\frac{\pi e\lambda_{B}}{2\epsilon_{\rm
pha}}\right\rceil,$ (154)

where $\epsilon_{\rm pha}$ is the allowable error due to phase estimation. For
the application to the electronic structure problem, we take
$\lambda_{B}=\lambda_{U}+\lambda_{V}$. As described above, we will not be
using exactly the optimal size of time step, but the increase in the cost is
upper bounded by $1+1/2^{2b_{T}+1}$, giving

$\mathcal{N}=\left\lceil(1+1/2^{2b_{T}+1})\frac{\pi
e\lambda_{B}}{2\epsilon_{\rm pha}}\right\rceil.$ (155)

###  III.4 Total cost of the interaction picture simulation

In the block encoding of the Dyson series, there are a number of parts that
need to be costed.

  1. 1.

Preparing the superposition over $k$ in unary.

  2. 2.

The superpositions on the time registers will be prepared with controlled
Hadamards.

  3. 3.

The $K$ registers can be sorted with a number of comparators as given in Table
3.

The contributions to the error are as follows.

  1. 1.

The finite cutoff on the Dyson series, $K$.

  2. 2.

The finite angles of the rotations for preparing the superposition over $k$.

  3. 3.

The finite number of bits used to represent the times, $n_{t}$.

For the time steps we are using for the qubitized Dyson series approach, the
time interval is $\tau=1/\lambda_{B}$, which means that the size of term $k$
in the Dyson series is simply $1/k!$, and the error is approximately
$1/(K+1)!$. More precisely, the error is

$e-\sum_{k=0}^{K}\frac{1}{k!}.$ (156)

This is error in block-encoding $\sin((A+B)\tau)=\sin((A+B)/\lambda_{B})$, and
therefore corresponds to error in the estimation of the energy of

$\lambda_{B}\left(e-\sum_{k=0}^{K}\frac{1}{k!}\right)\approx\frac{\lambda_{B}}{(K+1)!}.$
(157)

Note that here we are assuming that the eigenvalue is near zero, so the sine
function can be linearized in estimating the error. In the case where $E_{j}$
(the energy eigenvalue to be measured) is not near zero, it is possible to
just apply a phase shift so we are block encoding $\sin((A+B-E_{0})\tau)$ with
$E_{0}$ near $E_{j}$.

For the error in preparing the superposition over $k$, we note that because
the amplitudes we need are precisely $1/\sqrt{k!}$, the state can be written
as proportional to

$\sum_{k=0}^{K}\sqrt{\frac{K!}{k!}}\lvert k\rangle,$ (158)

where the amplitudes are the square roots of integers. The $-\mathrm{i}$
factors in the Dyson series can be implemented just with $S$ gates on the
qubits with the unary encoding of $k$, so do not require non-Clifford gates.
We will also be controlling between forward and reverse time evolution, which
means that we will be controlling between $\pm\mathrm{i}$. That can be
achieved simply with a controlled phase gate, and again there is no non-
Clifford cost. We can therefore prepare an equal superposition over a number
of integers $\Sigma(0)\approx K!e$, where we define

$\Sigma(k)=\sum_{\ell=k}^{K}\frac{K!}{\ell!},$ (159)

and use inequality testing to exactly obtain the required state. The equal
superposition over the integers is not performed exactly, but the failure,
which has small probability, is flagged so it only corresponds to slightly
reducing the size of the block encoded operation, rather than introducing any
error into that block encoded operation.

For the complexity, the contributions are as follows.

  1. 1.

First there is the preparation of an equal superposition over a number of
integers that is upper bounded by $K!e$. The number of bits can be taken to be
$n_{k}=\lceil\log(\Sigma(0))\rceil$, so the cost of preparing the equal
superposition is (according to the usual formula)

$3n_{k}+2b_{r}-9.$ (160)

  2. 2.

To select $k=0$ or $k=1$, we can perform an inequality test between the equal
superposition register and $\Sigma(2)$. The Toffoli cost of this inequality
test is $n_{k}-1$. To distinguish $k=0$ and $k=1$, we simply use a Hadamard on
an ancilla qubit and use a Toffoli to set the qubit.

  3. 3.

For $1<k<K$, we perform an inequality test with the number $\Sigma(k+1)$. To
save complexity, we can only perform an inequality test with a number of bits
necessary to store $\Sigma(k)$, conditioned on already passing the inequality
test with $\Sigma(k)$ in the previous step. The Toffoli complexity may then be
given as $\lceil\log(\Sigma(k))\rceil-1$ for the inequality test, and another
Toffoli is needed for the conditioning.

The first contribution here is a complexity that must be paid again in the
inverse preparation, giving a total cost

$2(3n_{k}+2b_{r}-9).$ (161)

For contributions 2 and 3, the Toffolis and inequality tests may be inverted
with Clifford gates, giving a Toffoli cost

$n_{k}+\sum_{k=2}^{K-1}\lceil\log(\Sigma(k))\rceil.$ (162)

In order to ensure that these inequality tests and Toffolis can be inverted
with Cliffords, the same number of qubits are kept throughout the calculation.
For a more detailed explanation of the costing, see Appendix E.

The other contributor to the error listed above is the finite number of bits
used to represent the time, $n_{t}$. In prior work, this error was estimated
based on a bound on the derivative of the Hamiltonian. Here we use a tighter
bound on the error based on the second derivative. This is a lengthy
derivation, so is relegated to Appendix H. The result is that, for time
interval $\tau$, the error in the time-ordered exponential can be bounded by

$\left(\frac{2\|A\|}{\|B\|}+\frac{\|A\|^{2}}{\|B\|^{2}}\right)\left[2^{n_{t}}\left(e^{\|B\|\tau/2^{n_{t}}}-1\right)-\|B\|\tau(1+\|B\|\tau/2^{n_{t}+1})\right].$
(163)

Here we take $\tau=1/(\lambda_{U}+\lambda_{V})$, and replace $\|B\|$ with
$\lambda_{U}+\lambda_{V}$ and $\|A\|$ with $\lambda_{T}$. That gives the bound

$\left(\frac{2\lambda_{T}}{\lambda_{U}+\lambda_{V}}+\frac{\lambda_{T}^{2}}{(\lambda_{U}+\lambda_{V})^{2}}\right)\left[2^{n_{t}}\left(e^{2^{-n_{t}}}-1\right)-(1+2^{-(n_{t}+1)})\right].$
(164)

Again, this is error in block encoding of $\sin((A+B)/\lambda_{B})$, and
corresponds to error in the estimation of the energy that is multiplied by
$\lambda_{B}=\lambda_{U}+\lambda_{V}$, so the error in the energy is

$\left(2\lambda_{T}+\frac{\lambda_{T}^{2}}{\lambda_{U}+\lambda_{V}}\right)\left[2^{n_{t}}\left(e^{2^{-n_{t}}}-1\right)-(1+2^{-(n_{t}+1)})\right].$
(165)

The complexity of the controlled Hadamards to generate the superpositions over
times is therefore $Kn_{t}$. The same complexity is needed in inverting the
state preparation, giving total cost

$2Kn_{t}.$ (166)

Sorting the times with $n_{t}$ bits gives complexity $2n_{t}$ times the
numbers of comparators given in Table 3. These numbers of comparators are
denoted ${\rm St}\left(K\right)$, and the same cost is needed to invert the
sort, giving a total sorting cost

$4n_{t}{\rm St}\left(K\right).$ (167)

For the times, we need to compute the differences. There is a subtlety here,
in that we need to compute the difference with the maximum time for the final
time difference for a given $k$. That can be achieved by taking all times
where we have _not_ performed the Hadamards to be set to the maximum time (all
ones). One can simply perform $X$ gates after the controlled Hadamards. The
total cost of the differences and absolute values was then given before as
$2(K-1)(n_{t}-1)$.

That accounts for the need to make the exponentiation of the kinetic operator
controlled on the $k$ register, but we also need to do this for the block
encoding of $U+V$. The block encoding is a combination of the $\mathrm{PREP}$
operation, a $\mathrm{SEL}$ operation, and an inverse preparation
$\mathrm{PREP}^{\dagger}$. Because the preparation and inverse preparation
cancel, only the $\mathrm{SEL}$ operation need be made controlled. The cost of
this control may be quantified as follows.

First, we need to control on success of the three equal superposition state
preparations and the preparation of $\nu$, as well as the qubit in the $k$
register. That takes 4 Toffolis to produce a qubit that is the AND of all five
bits. We will apply another 4 Toffolis in erasing this qubit (though we could
also keep the ancillas to erase it without Toffolis).

For the phasing by $-e^{-\mathrm{i}k_{\nu}\cdot R_{\ell}}$, the $R_{\ell}$ can
be made zero in the case where we do not want to perform $\mathrm{SEL}$ by
simply controlling the QROM for outputting of $R_{\ell}$ on that qubit. That
QROM is already controlled on the qubit selecting between $U$ and $V$, and
controlling on another qubit merely needs one more Toffoli. There also needs
to be control of the minus sign on that qubit as well. However, in controlling
the output of the QROM, we are generating an ancilla that can be used to
control that sign flip, and no additional Toffolis are required.

The addition $p+\nu$ is already controlled by the qubit selecting between $U$
and $V$. Again, making it controlled on one more qubit only takes one more
Toffoli. We have costed the subtraction $q-\nu$ as being controlled, but for
the interaction picture that subtraction only needs to be controlled on this
one qubit. We have already given the cost for this subtraction to be
controlled, so there is _no_ additional Toffoli cost over what we have already
presented. As a result, the total additional complexity to make the block
encoding of $U+V$ controlled on a qubit is just _two_ Toffolis, which is
completely trivial compared to the other costs in the circuit.

In the block encoding, we also need to reflect on all qubits used that are not
guaranteed to be returned to zero. The detailed analysis is given in Appendix
C, and yields a Toffoli cost of

$K(n_{\eta\zeta}+2n_{\eta}+4n_{p}+n_{\mathcal{M}}+n_{t}+12)+n_{k}+3.$ (168)

In this costing, we use intermediate checks between the block encodings of
$U+V$ to ensure that the control qubits are rezeroed, and we can reuse them.
Lastly we provide the total cost of a single step of the interaction picture
including inversion in Table 5. There, 10 Toffolis are added to the costs from
Table 4 to account for the complexity of making the operations for $U$ and $V$
controlled.

Procedure |  Toffoli cost
---|---
The preparation of the equal superposition over $\Sigma(0)$ numbers; see (161). |  $2(3n_{k}+2b_{r}-9)$
Inequality tests used for the preparation of the superposition over $k$; see (162). |  $n_{k}+\sum_{k=2}^{K-1}\lceil\log(\Sigma(k))\rceil$
Controlled Hadamards for preparing the superpositions over times; see (166). |  $2Kn_{t}$
The sorting of the time registers; see (167). |  $4n_{t}{\rm St}\left(K\right)$
The differences of the times; see (151). |  $2(K-1)(n_{t}-1)$
The kinetic operator needs to be exponentiated $K+1$ times; see (152). |  $(K+1)[2n_{t}(n_{\eta}+2n_{p})-n_{t}+b_{\rm grad}-2]+2(b_{\rm grad}-2)$
Perform the controlled block encoding of the potential operator $K$ times. |  $K(10+{\rm costs~{}in~{}\hyperref@@ii[tab:intcosts1]{Table~{}\ref*{tab:intcosts1}}})$
The reflection cost for qubitization. |  (168)
Table 5: The costs involved in the block encoding of the interaction picture
representation of the Hamiltonian.

To determine the effective $\lambda$, we need to also account for the
probabilities of success of the amplitude amplification for the various state
preparations we have performed. For the block encoding of $U+V$, these are as
follows.

  1. 1.

The preparation of the state with amplitudes $1/\|\nu\|$, which has
probability of success denoted $p_{\nu}^{\rm amp}$.

  2. 2.

The preparation of the equal superposition state for selecting between $U$ and
$V$, which has probability ${\rm Ps}\left(\eta+2\lambda_{\zeta},b_{r}\right)$.

  3. 3.

Preparing the equal superpositions over $i$ and $j$ have probabilities of
success ${\rm Ps}\left(\eta,b_{r}\right)$ each.

  4. 4.

The probability of obtaining $i\neq j$ is $1-1/\eta$.

The value of $\lambda_{U}+\lambda_{V}$ needs to be divided by all the
probabilities in parts 1 to 3, and just $\lambda_{V}$ needs to be divided by
$1-1/\eta$. It is also possible to perform a joint state preparation over $i$
and $j$, eliminating the need to divide by $1-1/\eta$. We do not explain that
here as the procedure is already quite complicated.

In performing the block encoding of $\sin((U+V)\tau)$ there is a further
preparation of an equal superposition state in order to prepare the
superposition over $k$. That has probability of success ${\rm
Ps}\left(\Sigma(0),b_{r}\right)$. Because we are considering a region where
the sine can be linearized, the effect is the same as multiplying $U+V$ by
this probability, and we should therefore multiply all these probabilities.
Also taking into account the factor $(1+1/2^{2b_{T}+1})$ due to adjusting the
length of time, we get

$P_{\rm eq}=p_{\nu}^{\rm amp}{\rm Ps}\left(\Sigma(0),b_{r}\right){\rm
Ps}\left(\eta+2\lambda_{\zeta},b_{r}\right){\rm
Ps}\left(\eta,b_{r}\right)^{2}/(1+1/2^{2b_{T}+1}).$ (169)

The value of $\lambda_{U}+\lambda_{V}$ then needs to be divided by $P_{\rm
eq}$ to obtain the effective $\lambda$. Therefore, to get the total cost, the
costs for a single step of the interaction picture need to be multiplied by a
number of steps

$\mathcal{N}=\left\lceil\frac{\pi
e(\lambda_{U}+\lambda_{V}/(1-1/\eta))}{2\epsilon_{\rm pha}P_{\rm
eq}}\right\rceil.$ (170)

Now we take account of the nonlinearity in using an operator with eigenvalue
$\arcsin[\sin((E_{j}-E_{0})\tau)/e]$. As discussed above, the error of the
estimation is increased by a relative amount
$\mathcal{O}((\lambda_{U}+\lambda_{V})^{2}\Delta E^{2})$, so we need to
introduce this factor to the number of steps used. Also, the nonlinearity
gives an $\mathcal{O}(1)$ increase to the number of steps needed. If we have
an initial estimate of the energy with precision $\Delta E$, we can choose
$E_{0}$ such that the number of steps needed is

$\mathcal{N}=\frac{\pi e(\lambda_{U}+\lambda_{V}/(1-1/\eta))}{2\epsilon_{\rm
pha}P_{\rm eq}}\left[1+\mathcal{O}((\lambda_{U}+\lambda_{V})^{2}\Delta
E^{2})\right]+\mathcal{O}(1).$ (171)

In this expression, we should also use the values
$\lambda^{\alpha}_{U},\lambda^{\alpha}_{V}$ given by the block encoding to be
precise, though in practice the difference should be negligible.

In computing these costs, we need to ensure that the total error is no larger
than $\epsilon$, where the sources of error are as follows.

  1. 1.

The error due to the phase estimation $\epsilon_{\rm pha}$.

  2. 2.

The error due to the truncation of the Dyson series,
$\epsilon_{K}\approx(\lambda_{U}+\lambda_{V})/(K+1)!$.

  3. 3.

The error $\epsilon_{R}$ due to the finite number of bits used to represent
the nuclear positions, which is upper bounded as in Lemma 2.

  4. 4.

The error $\epsilon_{\mathcal{M}}$ due to the preparation of the state with
amplitudes $1/\|\nu\|$, which is upper bounded as in Lemma 3.

  5. 5.

The error due to the finite number of time steps, $\epsilon_{t}$, which is
upper bounded as in (165).

Because the error $\epsilon_{\rm pha}$ is an RMS value, we combine it with the
other errors as a square root of a sum of squares. We therefore require that

$\epsilon_{\rm
pha}^{2}+(\epsilon_{K}+\epsilon_{R}+\epsilon_{\mathcal{M}}+\epsilon_{t})^{2}\leq\epsilon^{2}.$
(172)

To summarise, the result can be given as in the following theorem.

######  Theorem 5.

Using the interaction picture, assuming an initial estimate with accuracy
$\Delta E$, the eigenvalue of the Hamiltonian given by (8) can be estimated to
within RMS error $\epsilon$, using a number of Toffoli gates

$\displaystyle\mathcal{N}$
$\displaystyle\left(2(3n_{k}+2b_{r}-9)+n_{k}+\sum_{k=2}^{K-1}\lceil\log(\Sigma(k))\rceil+2Kn_{t}+4n_{t}{\rm
St}\left(K\right)\right.$
$\displaystyle\quad+2(K-1)(n_{t}-1)+(K+1)[2n_{t}(n_{\eta}+2n_{p})-n_{t}+b_{\rm
grad}-2]+2(b_{\rm grad}-2)$
$\displaystyle\quad+K\left\\{10+2(4n_{\eta\zeta}+2b_{r}-9)+(14n_{\eta}+8b_{r}-36)+3[3n_{p}^{2}+15n_{p}-7+4n_{\mathcal{M}}(n_{p}+1)]\right.$
$\displaystyle\quad\left.+\lambda_{\zeta}+{\rm
Er}\left(\lambda_{\zeta}\right)+(12\eta
n_{p}+4\eta-4)+24n_{p}+6n_{p}n_{R}+(12n_{p}^{2}+2n_{p}+8n_{\eta})\right\\}$
$\displaystyle\quad\left.\vphantom{\sum_{k=2}^{K-1}}+K(n_{\eta\zeta}+2n_{\eta}+2+4n_{p}+n_{\mathcal{M}}+n_{t}+12)+n_{k}+3+\mathcal{O}(\log(1/\epsilon))\right),$
(173)

where

$\displaystyle\mathcal{N}$ $\displaystyle=\frac{\pi
e(\lambda^{1}_{U}+\lambda^{1}_{V}/(1-1/\eta))}{2\epsilon_{\rm pha}P_{\rm
eq}}\left[1+\mathcal{O}((\lambda_{U}+\lambda_{V})^{2}\Delta
E^{2})\right]+\mathcal{O}(1),$ (174) $\displaystyle P_{\rm eq}$
$\displaystyle=p_{\nu}^{\rm amp}{\rm Ps}\left(\Sigma(0),b_{r}\right){\rm
Ps}\left(\eta+2\lambda_{\zeta},b_{r}\right){\rm
Ps}\left(\eta,b_{r}\right)^{2}/(1+1/2^{2b_{T}+1}),$ (175) $\displaystyle
b_{\rm grad}$
$\displaystyle=b_{T}-\left\lceil\log\left(\frac{\pi}{(\lambda_{U}+\lambda_{V})\Omega^{2/3}}\right)\right\rceil,$
(176)

with $b_{T}\in\mathbb{N}$. For this complexity we require appropriately chosen
phase gradient states to be given as a resource. The quantities $\epsilon_{\rm
pha}>0$ and $K,n_{R},n_{\mathcal{M}},n_{T}\in\mathbb{N}$ are chosen such that

$\displaystyle\epsilon^{2}$ $\displaystyle\geq\epsilon_{\rm
pha}^{2}+(\epsilon_{K}+\epsilon_{R}+\epsilon_{\mathcal{M}}+\epsilon_{t})^{2},$
(177) $\displaystyle\epsilon_{K}$
$\displaystyle=(\lambda_{U}+\lambda_{V})\left(e-\sum_{k=0}^{K}\frac{1}{k!}\right),$
(178) $\displaystyle\epsilon_{\mathcal{M}}$
$\displaystyle=\frac{2\eta}{2^{n_{\mathcal{M}}}\pi\Omega^{1/3}}\left(\eta+2\sum_{\ell=1}^{L}\zeta_{\ell}\right)\left(7\times
2^{n_{p}+1}-9n_{p}-11-3\times 2^{-n_{p}}\right)$ (179)
$\displaystyle\epsilon_{R}$
$\displaystyle=\frac{\eta\lambda_{\zeta}}{2^{n_{R}}\Omega^{2/3}}\sum_{\nu\in
G_{0}}\frac{1}{\left\lVert\nu\right\rVert},$ (180) $\displaystyle\epsilon_{t}$
$\displaystyle=\left(2\lambda_{T}+\frac{\lambda_{T}^{2}}{\lambda_{U}+\lambda_{V}}\right)\left[2^{n_{t}}\left(e^{2^{-n_{t}}}-1\right)-(1+2^{-(n_{t}+1)})\right].$
(181)

###### Proof.

The Toffoli costs given in this theorem are simply the sum of the costs found
above in Table 5. We take $\epsilon_{\rm pha}$, $\epsilon_{K}$,
$\epsilon_{R}$, $\epsilon_{\mathcal{M}}$, and $\epsilon_{t}$ to be the
allowable errors from phase estimation, truncation of the Dyson series at
order $K$, the finite approximation of the positions of the nuclei,
preparation of the state with amplitudes $1/\|\nu\|$, and finite number of
bits for the time, respectively. Because the error in the phase estimate is
independent of all other errors, the sum of the squares of this error and all
other errors needs to be smaller than $\epsilon^{2}$. The expression for
$\epsilon_{K}$ is from (157), the expression for $\epsilon_{R}$ is from Lemma
2, the expression for $\epsilon_{\mathcal{M}}$ is from Lemma 3, and the
expression for $\epsilon_{t}$ is obtained from the upper bound in (165). The
probability of success of the state preparations is given by (169), and the
number of repetitions needed for the state preparation is given by (171). We
have used $\lambda^{\alpha}_{U},\lambda^{\alpha}_{V}$ with $\alpha=1$ so we
can use Lemma 3. The term $\mathcal{O}(\log(1/\epsilon))$ is to account for
complexity of preparing the initial state to use as control for the phase
measurement, and the processing of this control at the end to obtain the phase
estimate. ∎

As discussed previously for the case of qubitization, the general way to use
this theorem is to choose some allowable maxima for the errors, compute the
appropriate values of the parameters required to obtain those errors, then
compute the actual error to obtain $\epsilon_{K}$, $\epsilon_{\mathcal{M}}$,
$\epsilon_{R}$, and $\epsilon_{t}$. Those are then used in (177) to obtain
$\epsilon_{\rm pha}$, which is then used to obtain the total complexity. The
errors $\epsilon_{\mathcal{M}}$ and $\epsilon_{R}$ were discussed above for
Theorem 4. For the case of $\epsilon_{K}$, we can choose $K$ proportional to
$\ln[(\lambda_{U}+\lambda_{V})/\epsilon_{V}]/\ln\ln[(\lambda_{U}+\lambda_{V})/\epsilon_{V}]$.
For $\epsilon_{t}$ we can use the approximation that the expression in the
square brackets is about $2^{-2n_{t}}/6$ to estimate the appropriate value of
$n_{t}$. The value of $b_{T}$ is chosen to be an integer, usually around 8. We
can then adjust the values used for $n_{\mathcal{M}}$, $n_{R}$, $K$, $n_{t}$,
and $b_{T}$ to minimise the complexity.

##  IV Results and discussion

###  IV.1 The resources required to simulate real systems

In this section we numerically compute and compare the costs of the central
algorithms of the preceding two sections to each other, and to other methods
in the literature. The constant factor Toffoli complexity is given for our
qubitization approach by Theorem 4 and given for our interaction picture
approach by Theorem 5. While those expressions are rather complex, the simpler
asymptotic complexities are given in (26) and (33), respectively. We see that
we can characterize the complexity of these approaches in terms of $\eta$,
$N$, $\Omega$ and $\epsilon$. The number of nuclear charges $L$ also appears
in the scaling but only inside of subdominant terms that do not significantly
impact the complexity. Thus, our algorithm has no leading order dependence on
the location, magnitude, or number of nuclear charges. This is a stark
contrast to algorithms based on Gaussian orbitals where the basis itself, and
hence the $\lambda$ values as well as the cost of the quantum walk steps, are
greatly influenced by the particular details of the system being simulated.
Because our costs do not depend strongly on such properties, we will be able
to analyze the scaling of our algorithms much more broadly than is possible
with the molecular orbital based algorithms.

|
---|---
(a) Qubitization Toffoli count | (b) Interaction picture Toffoli count
|
(c) Qubitization logical qubit count | (d) Interaction picture logical qubit count
Figure 1: Toffoli (top) and logical qubit (bottom) counts of the qubitization algorithm (left) and interaction picture algorithm (right) as a function of $\eta$. Here we take $\epsilon=0.0016$ Hartree (corresponding to “chemical accuracy”) and $N=2^{18}\approx 2.6\times 10^{5}$ but we note that outside of the dependence on $\Delta$, algorithm costs depend only logarithmically on $N$. For simplicity these plots correspond to the uniform electron gas (i.e., $U=0$) but block encoding the external potential adds only marginal overhead so we expect these costs to be broadly representative of the resources required to simulate arbitrary molecular systems as a function of $\eta$ and $\Delta$. We emphasize that $\Delta=10^{-3}a_{0}$ is an excessively accurate discretization that would correspond to a grid of one billion points per each Bohr radius cubed of system volume. In most cases, we expect that a resolution of $\Delta=10^{-2}a_{0}$ (which is one million grid points per Bohr radius cubed) would correspond to more accurate energies than large correlation consistent Gaussian basis sets. In the bottom panels we see that the qubit counts of the two algorithms are similar. |
---|---
(a) Ratio of Toffoli counts | (b) Parameter regions of advantage
Figure 2: (a) The ratio of the number of Toffolis for the interaction picture
approach to the number of Toffolis for the qubitization approach. We take
$\epsilon=0.0016$ Hartree (corresponding to “chemical accuracy”) and
$N=2^{18}\approx 2.6\times 10^{5}$. We again note that outside of the
dependence on $\Delta$, algorithm costs depend only logarithmically on $N$.
For simplicity these plots correspond to the uniform electron gas (i.e.,
$U=0$) but block encoding the external potential adds only marginal overhead
so we expect these costs to be broadly representative of the resources
required to simulate arbitrary molecular systems as a function of $\eta$ and
$\Delta$. (b) The parameter region of $\Delta$ versus $\eta$ where there is an
advantage for the Toffoli complexity for the interaction picture approach
versus qubitization. The region where the interaction picture has the
advantage is shown in red, and the region where qubitization has the advantage
is shown in blue. That is, the lower left of (b) is showing the values of
$\Delta$ for which the cost ratio plotted in panel (a) is less than one, as a
function of $\eta$. In this figure we show the crossover for four different
values of $N$, demonstrating that it barely changes as a function of $N$. In
both figures we see that the interaction picture approach only has an
advantage for low particle numbers and/or when simulations target unusually
high precision.

While the logarithmic factors in the complexity sometimes depend on $N$
without $\Omega$ (a result of the $\log N$ sized ancilla and momenta
registers), in the leading order complexity we see that $N$ and $\Omega$
always appear together as a ratio. For example, assuming fixed precision (in
these numerics we will target chemical accuracy, $\epsilon=0.0016$ Hartree) we
can express the leading order algorithmic gate complexities as

${\rm gates}\left({\rm qubitization}\right)=\widetilde{\cal
O}\left(\frac{\eta^{2}}{\Delta^{2}}+\frac{\eta^{3}}{\Delta}\right)\qquad{\rm
gates}\left({\rm interaction\,\,pic}\right)=\widetilde{\cal
O}\left(\frac{\eta^{3}}{\Delta}\right)\qquad\Delta=\left(\frac{\Omega}{N}\right)^{1/3}$
(182)

where we have defined $\Delta$ as the spacing of the plane wave reciprocal
lattice. Intuitively, $\Delta$ can be thought of as the length scale of
resolution within our basis (e.g., wavefunction features smaller than $\Delta$
cannot be resolved but those larger than $\Delta$ can be resolved). It is easy
to see that the $1/\Delta$ dependence arises from the maximum energy of the
Coulomb operator between grid points separated by $\Delta$ whereas the
$1/\Delta^{2}$ dependence arises due to the maximum energy of the kinetic
operator which is proportional to the square of velocities on such a grid. If
one desires to model a non-periodic system using plane waves then the typical
strategy is to choose a sufficiently large box size $\Omega$ such that the
periodic images of the Coulomb operator do not interact555A more efficient
procedure for modeling non-periodic systems using plane waves is to use a
truncated Coulomb operator with a slightly larger cell size Rozzi  _et al._
(2006); Sundararaman and Arias (2013); Ismail-Beigi (2006) or to subtract
terms from the molecular multipole Makov and Payne (1995) or Martyna-Tuckerman
corrections Martyna and Tuckerman (1999). This can likely be incorporated into
our approach in a straightforward fashion but would require minor
modifications to the potential block encodings which we leave for future
work.. At first glace, the larger one makes the box, the cheaper our
algorithms due to the inverse dependence on $\Omega$. However, as $\Omega$
gets larger the resolution $\Delta$ decreases, and thus to maintain reasonable
precision one must also increase $N$. By expressing the complexity as in (182)
we see that if the goal is perform a calculation with fixed spatial resolution
$\Delta$, the exact size of the box will not substantially impact the
complexity.

In Figure 1 we compare the complexities of the two algorithms of this paper by
plotting the total Toffolis required against $\eta$ for several different
values of $\Delta$. While the interaction picture algorithm has slightly
better scaling in terms of $\Delta$, the constant factors are meaningfully
worse. To make the comparison between the complexities more directly, in
Figure 2 we show both the ratio of the Toffoli complexities and the parameter
regimes where each approach is more efficient. It would appear that the
qubitization algorithm has lower cost in almost all regimes of interest. The
exception would be for _extremely_ high resolution simulations where both
$\eta$ and $\Delta$ are on the smaller side of what we would imagine requiring
(there is little dependence on $N$). For example, the interaction picture
algorithm is cheaper when $\eta=20$ and $\Delta=10^{-3}a_{0}$. While such
precisions might be appropriate for performing dynamics, we not anticipate
that those precisions will typically be required for preparing and sampling
observables from molecular eigenstates. In terms of spatial resolution,
$\Delta=10^{-3}a_{0}$ is quite high since it corresponds to a reciprocal
lattice grid with a billion lattice points per Bohr radius cubed. Of course,
at higher values of $\eta$ or $\Delta$ the term with $\eta^{3}/\Delta$ scaling
dominates in the qubitization cost, which gives a scaling that is similar to
that of the interaction picture algorithm but with lower constant factors. A
priori, it was not obvious to us that the qubitization algorithm would be more
efficient than the interaction picture algorithm for modest values of $\eta$.
Thus, the fact that qubitization is usually more practical for chemistry
problems, is an important finding of our work.

system | conventional cell type | $\Omega/a_{0}^{3}$ |  $\eta$ (cell total) |  $r_{s}/a_{0}$ (total) |  $\eta$ (cell valence) |  $r_{s}/a_{0}$ (valence)
---|---|---|---|---|---|---
metallic lithium | body-centered-cubic | 284.94 | 6 | 2.25 | 2 | 3.24
metallic potassium | body-centered-cubic | 961.67 | 38 | 1.82 | 2 | 4.86
diamond | diamond cubic | 307.04 | 48 | 1.15 | 32 | 1.32
crystalline silicon | diamond cubic | 1,080.43 | 112 | 1.32 | 32 | 2.01
iron (II) oxide | NaCl structure | 539.84 | 136 | 0.98 | 52 | 1.35
cobalt oxide | NaCl structure | 522.81 | 140 | 0.96 | 60 | 1.28
aluminium arsenide | zincblende | 1,197.86 | 184 | 1.16 | 32 | 2.08
indium phosphide | zincblende | 1,364.93 | 256 | 1.08 | 32 | 2.17
Table 6:  Values of $\eta$, $\Omega$ and $r_{s}$ for a variety of common materials. Above, $a_{0}$ is the Bohr radius of the hydrogen atom and thus $r_{s}$ values are in units of Bohr and $\Omega$ values are units of ${\rm Bohr}^{3}$. Often one would want to perform a simulation with multiple unit cells in order to resolve correlations in between the units. Doing this would not change $r_{s}$ but it would increase $\eta$ proportional to the number of unit cells included. The valence $\eta$ and $r_{s}$ values correspond to simulations with core electrons removed. Doing this would require the use of pseudopotentials, which is beyond the scope of the current paper but likely an important next step. Note further that the unit cell volumes tabulated here correspond to conventional unit cells, rather than primitive unit cells. Using primitive unit cells would necessitate adapting our approach to reciprocal spaces beyond the cubic Bravais lattices considered here. |
---|---
(a) $N=2^{12}\approx 4.1\times 10^{3}$ | (b) $N=2^{15}\approx 3.3\times 10^{4}$
|
(c) $N=2^{18}\approx 2.6\times 10^{5}$ | (d) $N=2^{21}\approx 2.1\times 10^{6}$
Figure 3: $\mathrm{Toffoli}$ count of qubitization as a function of $\eta$ using the improved techniques of this paper. In contrast to Figure 1, we plot these costs for fixed values of the Wigner-Seitz radius $r_{s}$ to give a sense of the scaling when $\eta$ increases proportionally to $\Omega$ (as would be the case when growing condensed phase simulations towards their thermodynamic limit). Thus, we expect these plots to be representative of costs when simulating periodic materials. We see that the qubitization based algorithm performs better across all regimes studied to compared to the interaction picture based algorithm (Figure 4). |
---|---
(a) $N=2^{12}\approx 4.1\times 10^{3}$ | (b) $N=2^{15}\approx 3.3\times 10^{4}$
|
(c) $N=2^{18}\approx 2.6\times 10^{5}$ | (d) $N=2^{21}\approx 2.1\times 10^{6}$
Figure 4: $\mathrm{Toffoli}$ count of interaction picture algorithm as a
function of $\eta$ using improved techniques. In contrast to Figure 1, we plot
these costs for fixed values of the Wigner-Seitz radius $r_{s}$ to give a
sense of the scaling when $\eta$ increases proportionally to $\Omega$ (as
would be the case when growing condensed phase simulations towards their
thermodynamic limit). Thus, we expect these plots to be representative of
costs when simulating periodic materials. We see that the scaling depends
rather weakly on the number of plane waves, $N$.

Next, we turn to the topic of how our algorithm scales for the important
application of simulating periodic systems such as crystals, metals, surfaces
and polymers. Such systems are ubiquitous in nature and are a major focus of
electronic structure methods deployed to materials science. While it is
important to understand how our algorithm scales as one approaches the
continuum basis set limit where $\Delta\rightarrow 0$ (as discussed in the
context of Figure 1) for condensed phase applications it is also important to
understand how our algorithm scales towards the thermodynamic limit. That is,
the computational cell size should be larger than the correlation length in
order for the simulation to faithfully capture the behavior of the infinite
system. The computational cell will be composed of multiple unit cells (the
simplest repeating unit in the crystalline lattice). For such simulations, one
can imagine adding more unit cells to the computational cell but keeping the
overall density of particles (i.e., $\eta/\Omega$) fixed. Often, behavior in
the thermodynamic limit is extrapolated from a series of progressively larger
calculations. In order to discuss simulations with fixed $\eta/\Omega$ in more
intuitive terms, we will analyze scaling in terms of the Wigner-Seitz radius,
$r_{s}$ (a common measure of a system’s average radius between electrons),
defined through the relation

$\frac{4}{3}\pi r_{s}^{3}=\frac{\Omega}{\eta},\qquad\qquad
r_{s}=\left(\frac{3\Omega}{4\pi\eta}\right)^{1/3}\,.$ (183)

We report Wigner-Seitz radii for several common materials in Table 6 and note
that they are mostly between one and a few Bohr radii. In terms of Wigner-
Seitz radii we can express the overall Toffoli complexities as

${\rm Toffolis}\left({\rm qubitization}\right)\\!=\widetilde{\cal
O}\\!\left(\frac{\eta^{4/3}N^{2/3}}{r_{s}^{2}}+\frac{\eta^{8/3}N^{1/3}}{r_{s}}\right)\qquad{\rm
Toffolis}\left({\rm interaction\,\,pic}\right)\\!=\widetilde{\cal
O}\\!\left(\frac{\eta^{8/3}N^{1/3}}{r_{s}}\right)\,.\qquad$ (184)

In Figure 3 and Figure 4 we plot the estimated resources required to perform
condensed phase simulations for the qubitization and interaction picture
algorithms, respectively. We make these plots between $\eta=20$ and
$\eta=200$, between $N=2^{12}$ and $N=2^{21}$ and between $r_{s}=0.1$ and
$r_{s}=100$. We observe that the qubitization algorithm is most practical due
to lower constant factors across all regimes plotted. As before, in the limit
of large $\eta$, the qubitization algorithm also recovers the better scaling
of the interaction picture algorithm (but with better constant factors). By
comparing Figure 3 with the requirements to simulate various materials in
Table 6 we see that our method can simulate classically intractable real
material Hamiltonians with between $10^{10}$ and $10^{12}$ Toffoli gates.
These numbers of Toffoli gates and qubits are comparable to those (for a
different system) in Lee  _et al._ (2021), where estimates are given of the
resources needed to implement the algorithm on a real device. If we adopt the
same assumptions here, including a surface code cycle time of 1 µs, our
proposed quantum algorithms would require on the order of a few million
physical qubits and a few days of runtime. No prior quantum algorithm for
chemistry has been shown to require anywhere near this few resources to
simulate realistic descriptions of condensed phase molecular systems.

We now address the question of how the more practical of the two algorithms
described in this work (the qubitization one) compares to other methods for
quantum simulating chemistry. To date, the only prior methods that have been
assessed for viability within fault-tolerance have used second quantization.
These papers fall into two categories: those that focus on plane waves (see
Table 1) and those that focus on Gaussian orbitals. The plane wave algorithms
assessed for fault-tolerance are by Babbush _et al._ Babbush  _et al._ (2018a)
and Kivlichan _et al._ Kivlichan  _et al._ (2020). Because the present work
also uses plane waves it is straightforward to make a comparison. As discussed
in Section I.1, for simulating molecules these prior papers using plane waves
are not particularly practical options because at least tens or hundreds of
thousands of plane waves would be required for reasonably accurate simulations
and those works have a space complexity scaling linearly in $N$. Those
algorithms also have significantly worse Toffoli complexity, scaling roughly
like $\widetilde{\cal O}(N^{3})$.

The main motivation for those past papers using plane waves in second
quantization is essentially a different regime than that of realistic
simulations of materials. Those papers emphasize their resource estimates in
terms of the first classically intractable applications of very small fault-
tolerant processors (perhaps only a few hundred thousand physical qubits). In
this sense, the goal is to articulate example problems on which “quantum
supremacy” Preskill (2012); Boixo  _et al._ (2018) might be first demonstrated
for quantum chemistry in an error-corrected context, rather than to highlight
applications corresponding to realistic simulations of nature. For example, at
small $N$ and high filling fraction $\eta/N$, there are a variety of model
systems that are not converged to the continuum limit but provide examples of
electronic structure that would be classically intractable to simulate within
the basis. The value of such simulations for the jellium model (the
homogeneous electron gas that occurs when $U=0$) near $r_{s}=10a_{0}$, is
argued in Section III of Babbush  _et al._ (2018b). As a result, this model
has been used as a benchmark in several subsequent papers. For example, the
largest jellium simulation reported in Table III of Babbush  _et al._ (2018a)
would require $2.1\times 10^{10}$ Toffoli gates666Many of the prior works we
compare to, including Kim _et al._ Kim  _et al._ (2021), Kivlichan _et al._
Kivlichan  _et al._ (2020) and Babbush _et al._ Babbush  _et al._ (2018a),
report the number of T gates required rather than the number of Toffoli gates
required. In order to facilitate a clear comparison we assume that the cost of
a Toffoli gate is equal to the cost of 2 T gates. This equivalence is
justified by the surface code spacetime volume required for a Toffoli gate
compared to a T gate using the catalyzed T gate distillation scheme of Gidney
and Fowler (2019). and $1,\\!136$ logical qubits when $r_{s}=10a_{0}$. This
simulation would involve just $1,\\!024$ plane wave spin-orbitals (equivalent
here to $N=512$). The largest $r_{s}=10a_{0}$ jellium system analyzed in Table
II of Kivlichan _et al._ Kivlichan  _et al._ (2020) involves only $N=125$
plane waves but would require $378$ logical qubits and $2.8\times 10^{9}$
Toffolis. Because those simulations occur at half-filling ($\eta=N$), they are
likely more efficient than the methods of this paper at those system sizes.
However, if we instead target a more physically accurate simulation (with less
basis set discretization error) where $N\geq 10^{3}$ and $\eta/N\leq 1/10$, we
would expect the first quantized approach of this paper to require far fewer
error-corrected resources than these past algorithms utilizing second
quantization.

To compare the the complexities reported here, the smallest value of $N$ we
consider here is for $N=2^{12}$, which is $8$ times larger than in Babbush
_et al._ (2018a). The $N^{3}$ scaling of the second quantization approach
would mean that the complexity is about $10^{13}$ Toffolis, nearly 100 times
larger than the largest complexity we have for $r_{s}=10a_{0}$ for the first
quantized qubitization approach. For the largest value of $N$ considered of
$2^{21}$, the $N^{3}$ scaling would indicate a complexity of over $10^{21}$
Toffolis, around a billion times larger than the complexity of the first
quantized approach here.

Unlike second quantized plane wave methods, second quantized algorithms using
Gaussian basis functions seem promising options for practical quantum
simulations of realistic molecular Hamiltonians. There are now a handful of
papers presenting in depth analysis of new fault-tolerant algorithms Reiher
_et al._ (2017); Berry  _et al._ (2019); von Burg  _et al._ (2020); Lee  _et
al._ (2021) in this representation. Dating back to the first of these papers
by Reiher _et al._ Reiher  _et al._ (2017), these papers have demonstrated
algorithmic improvements by estimating the resources required to simulate a
few benchmark molecules such as the FeMoCo molecule relevant to biological
nitrogen fixation. The work by Lee _et al._ Lee  _et al._ (2021) holds the
record for the Gaussian orbital based algorithms with both the lowest constant
factor and lowest asymptotic Toffoli scaling. We emphasize that all of these
papers have focused on simulating molecular systems in active spaces, and none
of them have developed methods for simulating periodic materials. As a result
it is difficult to directly compare the algorithms developed here to those
prior papers.

Active space Hamiltonians are coarse grained versions of larger molecular
Hamiltonians, designed to capture the physics of a subset of the total
electronic degrees of freedom. Often the active space electrons are the
valence electrons involved in the breaking or formation of chemical bonds. For
example, the FeMoCo molecule (${\rm Fe}_{7}{\rm Mo}{\rm S}_{9}{\rm C}$) is a
large system with 374 electrons; however, the active spaces studied by these
papers require only 54 Reiher  _et al._ (2017) or 113 Li  _et al._ (2019)
electrons. A variety of methods are available for computing active space
Hamiltonians Stein and Reiher (2019); Sayfutyarova  _et al._ (2017); Keller
_et al._ (2015); Bao and Truhlar (2019); Hermes  _et al._ (2020) that are
compatible with molecular orbital based quantum algorithms. Most such
approaches involve removing from the basis molecular orbitals that are
expected to be completely occupied in the ground state and adding the
potential arising from the charge density of those removed electrons to $U$.
Unfortunately, such methods are not immediately compatible with the algorithms
developed in this paper. This is because the molecular orbitals occupied by
the core electrons are not orthogonal to the rest of the plane wave basis and
thus cannot be easily removed from the basis.

However, there is a way to make similar approximations for plane wave methods.
In particular, pseudopotentials Tosoni  _et al._ (2007); Booth  _et al._
(2016) provide a well studied methodology for removing core electrons from
plane wave calculations by modifying the external potential to account for the
core. In addition to reducing the effective system size $\eta$,
pseudopotentials dramatically accelerate the convergence of the basis error.
The reason for this is because the cusps in electronic wavefunction which
occur at the nuclei are smoothed out by pseudpotentials. Because the Fourier
transform of smooth functions is exponentially convergent, the plane wave
basis can then converge dramatically faster. Thus, we suggest that the use of
pseudopotentials Kleinman and Bylander (1982); Dal Corso (2014) is likely to
be essential for the practical deployment of these methods. However, as the
present paper is already rather long, we leave the topic of removing core
electrons via pseudopotentials to future work. Instead, we will compare our
algorithms to Gaussian orbital resource estimates that do not use active
spaces.

reference | system | basis | $\Omega$ | $\eta$ | $N$ | logical qubits | Toffoli count
---|---|---|---|---|---|---|---
Kim _et al._ Kim  _et al._ (2021) | ethylene carbonate | STO-3G | N/A | 46 | 34 | 2,685 | $3.2\times 10^{10}$
Kim _et al._ Kim  _et al._ (2021) | ethylene carbonate | 6-311G | N/A | 46 | 90 | 14,492 | $8.5\times 10^{11}$
Kim _et al._ Kim  _et al._ (2021) | ethylene carbonate | cc-pVDZ | N/A | 46 | 104 | 16,698 | $1.3\times 10^{12}$
Kim _et al._ Kim  _et al._ (2021) | ethylene carbonate | cc-pVTZ | N/A | 46 | 236 | 81,958 | $3.1\times 10^{13}$
This work | ethylene carbonate | plane waves | $10^{5}a_{0}^{3}$ | 46 | $4.2\times 10^{3}$ | 1,395 | $2.5\times 10^{10}$
This work | ethylene carbonate | plane waves | $10^{5}a_{0}^{3}$ | 46 | $3.3\times 10^{4}$ | 1,701 | $6.6\times 10^{10}$
This work | ethylene carbonate | plane waves | $10^{5}a_{0}^{3}$ | 46 | $2.6\times 10^{5}$ | 2,021 | $1.7\times 10^{11}$
This work | ethylene carbonate | plane waves | $10^{5}a_{0}^{3}$ | 46 | $2.1\times 10^{6}$ | 2,355 | $4.2\times 10^{11}$
Kim _et al._ Kim  _et al._ (2021) | ${\rm LiPF}_{6}$ | STO-6G | N/A | 72 | 44 | 3,507 | $1.0\times 10^{11}$
Kim _et al._ Kim  _et al._ (2021) | ${\rm LiPF}_{6}$ | 6-311G | N/A | 72 | 112 | 18,423 | $1.9\times 10^{12}$
Kim _et al._ Kim  _et al._ (2021) | ${\rm LiPF}_{6}$ | cc-pVDZ | N/A | 72 | 116 | 18,600 | $1.7\times 10^{12}$
Kim _et al._ Kim  _et al._ (2021) | ${\rm LiPF}_{6}$ | cc-pVTZ | N/A | 72 | 244 | 84,721 | $3.3\times 10^{13}$
This work | ${\rm LiPF}_{6}$ | plane waves | $10^{5}a_{0}^{3}$ | 72 | $4.1\times 10^{3}$ | 1,758 | $8.0\times 10^{10}$
This work | ${\rm LiPF}_{6}$ | plane waves | $10^{5}a_{0}^{3}$ | 72 | $3.3\times 10^{4}$ | 2,150 | $2.1\times 10^{11}$
This work | ${\rm LiPF}_{6}$ | plane waves | $10^{5}a_{0}^{3}$ | 72 | $2.6\times 10^{5}$ | 2,556 | $5.1\times 10^{11}$
This work | ${\rm LiPF}_{6}$ | plane waves | $10^{5}a_{0}^{3}$ | 72 | $2.1\times 10^{6}$ | 2,976 | $1.3\times 10^{12}$
Table 7:  Table comparing the cost of the qubitization algorithm of this paper
to the cost of a state-of-the-art second quantized Gaussian basis algorithm
recently analyzed by Kim _et al._ Kim  _et al._ (2021). The work of Kim _et
al._ estimates the resources required to simulate systems, which despite being
small single molecules, are presented as relevant to the chemistry of lithium
batteries. Their estimates are a natural point of comparison because they do
not employ active spaces, unlike the focus of most past work using molecular
orbitals Reiher  _et al._ (2017); Berry  _et al._ (2019); von Burg  _et al._
(2020); Lee  _et al._ (2021); Elfving  _et al._ (2020). We focus on two
representative molecules from that work: ethylene carbonate and ${\rm
LiPF}_{6}$, discretized using four different Gaussian orbital basis sets. We
note that the algorithm of Lee _et al._ Lee  _et al._ (2021) has been shown to
scale better than the algorithm of von Burg _et al._ von Burg  _et al._
(2020), on which the work by Kim _et al._ is based, but all resource estimates
using that approach involve active spaces. We expect that in most contexts, if
one uses a number of plane waves that is more than a thousand times the number
of orbitals in a given Gaussian basis set, then one can expect similar or
higher precision from the larger plane wave calculation, even for small non-
periodic systems like the ones analyzed here. This would suggest that for
$N=2.6\times 10^{5}$ plane waves, our algorithms may achieve higher precision
than even the correlation consistent triple zeta (cc-pVTZ) basis set
calculations studied by Kim _et al._ Kim  _et al._ (2021). Thus, the results
in this table may suggest that our algorithms requires significantly less
surface code spacetime volume to reach similar levels of precision compared to
some of the most performant Gaussian based calculations.

A recent paper by Kim _et al._ Kim  _et al._ (2021) provides resource
estimates based on the algorithm of von Burg _et al._ von Burg  _et al._
(2020) (with modifications from Lee _et al._ Lee  _et al._ (2021)) that
happens to not use active spaces, thus enabling a comparison to the present
work which we provide in Table 7. However, it is still difficult to know
precisely how many plane waves would be required to match the accuracy of a
given Gaussian orbital calculation. While both suppress basis error as ${\cal
O}(1/N)$, even if a given Gaussian orbital calculation is more accurate than a
given plane wave calculation at fixed $N$, plane waves tend to converge more
systematically than Gaussian orbitals to the continuum limit. As a result, one
can more reliably predict the complete basis set limit energy from
extrapolations made from a series of calculations in progressively larger
basis sets when using plane waves.

Appendix E of Babbush  _et al._ (2018b) attempts to give an idea of the
relative prefactors between the number of plane waves and Gaussians required
to obtain similar levels of precision. While the precise ratio of plane waves
to Gaussians will depend on system specific details, there it is suggested
that periodic calculations with roughly $20\,N$ plane waves will often give
higher precision than the same calculations with $N$ Gaussian orbitals. It is
more difficult to simulate non-periodic systems like single molecules with
plane waves and reducing this periodicity requires yet more plane waves by
about a factor of $8$ Füsti-Molnar and Pulay (2002); Babbush  _et al._
(2018b). Achieving this convergence with plane waves may require softening
nuclear cusps with pseudopotentials. Moreover, when targeting extremely high
precision, the ratio of the required number of plane waves to Gaussians should
decrease because the challenge becomes resolving the electron-electron cusp.
Since the electron-electron cusp occurs at all points in space, Gaussians
centered on nuclei have no special advantages for resolving it. The fact that
Gaussians derive their advantage from being centered on nuclei also means that
as the number of atoms in a fixed volume increases, the number of Gaussian
orbitals to represent the system to fixed accuracy must also increase roughly
proportionally. On the other hand, plane waves are not biased towards
representing any particular nuclear locations; thus, as the number of atoms
per unit volume is increased, the ratio of the number of plane waves to the
number of Gaussians required for the same accuracy should decrease.

Considering all of these factors, Table 7 suggests that our qubitized first
quantized plane wave algorithm might be more performant than the best Gaussian
orbital basis algorithms applied to even small non-periodic molecules. For
example, when comparing our approach with $N=2.6\times 10^{5}$ plane waves to
the approach of Kim _et al._ Kim  _et al._ (2021) using a triple zeta basis
set for which $N$ is between 236 and 244 for the two systems included in Table
7, see that we are using more than a thousand times the number of plane waves
as Gaussian orbitals but that the first quantized plane wave approach needs
several dozen times fewer logical qubits and about a hundred times fewer
Toffoli gates (amounting to approximately three orders of magnitude less
surface code spacetime volume). Furthermore, it is clear that if the objective
is significantly lower basis set error (e.g., if spectroscopic precision of
$10^{-6}$ Hartree is desired) or if the goal is to simulate periodic systems
such as crystals, metals, polymers, or surfaces, then our first quantized
qubitization approach is dramatically more practical than other fault-tolerant
algorithms described in the literature. Thus, our analysis has shown that
these methods are an important and state-of-the-art tool for the quantum
simulation of chemistry.

###  IV.2 Conclusion and outlook

We have completed the first constant factor analysis of the resources required
to implement first quantized quantum simulations of quantum chemistry.
Unsurprisingly, these methods appear to be the most promising approach for
simulating realistic materials Hamiltonians within quantum error-correction.
More surprisingly, they are also highly competitive with the best second
quantized approaches for simulating single molecules (often requiring fewer
surface code resources by multiple orders of magnitude). We focused on
optimizing and compiling the first quantized algorithms using a plane wave
basis first described by Babbush _et al._ Babbush  _et al._ (2019). Despite
the interaction picture of that work having better scaling, we found that the
qubitization algorithm of that work tends to require fewer resources for the
regimes of interest to us.

Our work goes significantly beyond just translating the prior asymptotic
algorithm descriptions of Babbush _et al._ Babbush  _et al._ (2019) to
concrete circuits, even though just doing that would have still required a
similarly detailed manuscript. We also introduce many improvements to the
implementation of both the qubitized form and the interaction picture
simulation that reduce the complexity. The innovations we introduce that lead
to the most significant reduction in algorithmic complexity are summarized
below.

  * •

In the case of qubitizing the Hamiltonian, we entirely eliminate the cost of
the arithmetic for the kinetic energy term $T$ (see Section II.2). Instead of
computing the sum of squares of momenta, control registers are used to select
the momenta as well as the individual bits to multiply. The only arithmetic is
the multiplication of bits, which is performed with a single Toffoli. This is
a significant improvement, because multiplication would otherwise be the main
cost.

  * •

For the interaction picture simulations we greatly reduce the cost of
computing the total kinetic energy (see Section III.2). Instead of computing
it at each step, we compute the sum of squares of momenta at the beginning.
Then, each time we block encode the potential operator, which involves
updating just two momenta, we also update the register with the total momentum
by adding the change in momentum.

  * •

In both cases it is necessary to prepare a state with amplitudes $\propto
1/\|\nu\|$, where $\nu$ is a 3-dimensional vector of integers. This
preparation normally involves amplitude amplification, which triples the cost
of state preparation Brassard  _et al._ (2002). For qubitizing the
Hamiltonian, rather than using amplitude amplification, we can use the insight
that this state preparation is only needed for the potential term, but not for
the kinetic term. As a result, instead of using amplitude amplification, we
can implement the kinetic term in branches of the quantum state that would
otherwise correspond to failure (see Appendix D). That reduces the cost of
this term by a factor of 3.

  * •

In the case of the interaction picture one would normally require amplitude
amplification on the Dyson series in order to implement a step of the time
evolution deterministically, which again triples the cost. Instead, we can
qubitize a single step of the Dyson series (see Section III.1). This does not
deterministically implement time evolution, but gives a step with eigenvalues
related to those of the Hamiltonian. That eliminates one part where amplitude
amplification was required, but in that case we still need to use amplitude
amplification on the preparation of the state with amplitudes $\propto
1/\|\nu\|$.

  * •

Rather than performing arithmetic in place on the momentum registers, we swap
the correct momentum registers into ancilla registers controlled on ancillae
that select the registers (see Section II.4 and Section III.3). That means the
arithmetic only needs to be performed once in those ancillae, rather than for
every momentum register.

  * •

For the discretization of the integrals, we use times in the centers of the
intervals, which results in the error being higher order (see Appendix H).
With a careful analysis of the error due to the time discretization, the
number of bits needed to represent the time is about half of that using the
methods from prior work.

As a result of these improvements, in cases where there are a large number of
electrons the majority of the cost in each step is just selecting the momentum
registers. There is an overhead in the interaction picture approach
corresponding to the order of the Dyson series. That means that to provide an
improvement via the interaction picture, there should be a ratio between the
kinetic and potential energy that is at least as large as the order of the
Dyson series. That ratio is large in cases where the number of electrons is
small, but in those cases the cost of preparing the state with amplitudes
$1/\|\nu\|$ tends to be dominant. In those cases the direct qubitization
approach has the additional advantage that one can avoid the amplitude
amplification. That means the ratio needs to be even larger before there is an
advantage due to using the interaction picture. We find that there is a minor
improvement via the interaction picture with only $\eta=20$ electrons and
$N=2^{21}\approx 2\times 10^{6}$ plane waves, but as the number of electrons
grows this advantage diminishes.

We have demonstrated that our first quantized qubitization algorithm compares
favorably to all prior fault-tolerant algorithms for chemistry for both
periodic, and non-periodic systems. However, there are several obvious next
steps that should be pursued in order to further refine and generalize our
algorithms. While we have alluded to all of these ideas in the main text, we
also summarize them below.

  * •

We should modify the algorithm for block encoding $U$ to include various types
of accurate pseudopotentials. Using pseudopotentials will allow us to remove
core electrons from the simulation in order to only treat electronic degrees
of freedom that are involved in chemical reactivity and bonding, thus
effectively reducing $\eta$ while increasing $r_{s}$. Furthermore, the
screened potentials resulting from removing core electrons (even if only 1s
electrons are removed) are smooth functions. With sharp nuclear cusps removed
from the Coulomb potential (and thus, the wavefunction), the plane wave basis
will more rapidly converge to the continuum limit. It will be straightforward
to make this modification for spherically symmetric pseudopotentials such as
those that remove just the 1s electrons, but it will be more involved to adapt
the block encodings to more precise pseudopotentials for higher angular
momenta electrons (e.g., using non-local pseudopotentials). Including
pseudopotentials will also facilitate a more straightforward comparison to
Gaussian basis approaches using active spaces.

  * •

In (6) of this work we define the reciprocal lattice associated with our plane
wave mesh on a cubic lattice. However, for systems that lack a cubic
conventional unit cell (e.g., those with hexagonal symmetry), such lattices
cannot be used. Instead, we would need to define the reciprocal lattice using
non-orthogonal Bravais vectors. This would change the definition of $k_{p}$
from what is given in (6) to

$k_{p}=\frac{2\pi}{\Omega^{1/3}}\left(p_{1}g^{(1)}+p_{2}g^{(2)}+p_{3}g^{(3)}\right)\qquad\qquad
p=\left(p_{1},p_{2},p_{3}\right)\in G\,,$ (185)

where $G$ is defined in (6) and $g^{(1)}$, $g^{(2)}$ and $g^{(3)}$ are the
three dimensional non-orthogonal unit vectors of the Bravais lattice defined
by the _primitive_ unit cell geometry. Even for systems with cubic
conventional unit cells, there are likely convergence advantages associated
with being able to represent non-cubic primitive unit cells. However, using
this more general Bravais lattice will also come at a cost since it will make
computing $\|k_{p}\|^{2}$ (a bottleneck of our approach) more costly. Further
work should explore optimizations of these dot products in order to generalize
our methods while retaining high efficiency.

  * •

Related to the above point, there are different ways to converge the system
physics towards its thermodynamic limit. The approach explored here involves
growing the system in real space (i.e., by adding more unit cells to the
computational cell). However, one can also grow the system towards the
thermodynamic limit by evaluating an integral over the Brillouin zone in a
fashion referred to as “k-point sampling” Martin (2004). These schemes make
use of generalized Bravais lattices such as the Monkhorst-Pack mesh, which
allows one to represent systems of arbitrary symmetry while distributing the
reciprocal lattice vector points homogeneously in the Brillouin zone Martin
(2004). Thus far we have only described how to perform so-called
“$\Gamma$-point” calculations. Future work should investigate the benefits of
k-point sampling in the context of these algorithms and, if advantageous to do
so, make the relatively simple modification of the algorithms required to
support sampling outside of the $\Gamma$-point.

  * •

We have discussed that the methods of this paper can be used to simulate non-
periodic systems despite the underlying periodicity of plane waves. In
particular, this can be accomplish by using a sufficiently large cell volume
$\Omega$ so that periodic images do not interact. However, there are more
efficient ways of preventing these periodic images from interacting. For
example, one can use a range truncated Coulomb operator Rozzi  _et al._
(2006); Sundararaman and Arias (2013); Ismail-Beigi (2006). For an isolated
molecule, the total electronic density decays exponentially quickly away from
its center, and thus the molecule can be represented in a box of volume
$\Omega=D^{3}$ with only exponentially small parts of its true wavefunction
density outside of the box. By using a Coulomb kernel truncated at distance
$D$ Füsti-Molnar and Pulay (2002), such that

$V\left(r,r^{\prime}\right)=\begin{cases}\frac{1}{\left|r-r^{\prime}\right|}&\left|r-r^{\prime}\right|\leq
D\\\ 0&\left|r-r^{\prime}\right|>D,\end{cases}$ (186)

and by carrying out the simulation in a box of size $8\,\Omega=(2\,D)^{3}$, we
ensure that there is no Coulomb interaction at all between the repeated images
of the molecule, up to exponentially small terms in $\Omega$ arising from the
density of the molecule outside of the box. While the Fourier amplitudes of
the normal Coulomb kernel are $4\pi/k^{2}$, the Fourier amplitudes of the
truncated Coulomb interaction shown above are $4\pi(1-\cos(|k|D))/k^{2}$. This
would modify the potential operator coefficients of (7) to be

$U_{pq}\\!=\sum_{\ell=1}^{L}\frac{4\pi\left(1\\!-\\!\cos\left(\left\|k_{p-q}\right\|D\right)\right)\zeta_{\ell}e^{ik_{q-p}\cdot
R_{\ell}}}{\Omega\left\lVert k_{p-q}\right\rVert^{2}}\qquad
V_{pqrs}^{(\alpha,\beta)}\\!=\delta_{p-q,r-s}\frac{4\pi\left(1\\!-\\!\cos\left(\left\|k_{p-q}\right\|D\right)\right)\zeta_{\alpha}\zeta_{\beta}}{\Omega\left\|k_{p-q}\right\|^{2}}\,.$
(187)

Thus, future work focused on simulating systems of reduced periodicity should
adapt the algorithms presented here to block encode operators with the above
modification to the coefficients.

  * •

This paper has provided first quantized algorithms for preparing molecular
eigenstates given an initial state with sufficient overlap on an eigenstate of
interest. However, methods for state preparation in first quantization are
under-developed compared to their second quantized counterparts. The work of
Berry  _et al._ (2018) provided an algorithm for symmetrizing or
antisymmetrizing initial states with complexity ${\cal O}(\eta\log
N\log\eta)$, but another common step in state preparation would be to
initialize an arbitrary Slater determinant. Such states are described by
${\cal O}(\eta N)$ coefficients, and in second quantization can be prepared
with gate complexity of $\widetilde{\cal O}(\eta N)$ Kivlichan  _et al._
(2018). We would expect that this is also a lower bound on the gate complexity
required in first quantization. It would not be ideal for such state
preparations to need to scale linearly in $N$. However, one can choose to
involve fewer plane waves in the initial state preparation and so long as the
number of plane waves taken is less than $M=\widetilde{\cal
O}(\eta^{1/3}N^{2/3}/\epsilon+\eta^{5/3}N^{1/3}/\epsilon)$ then state
preparation will require a number of gates that is less than the ${\cal
O}(\eta M)$ of our qubitization algorithm, and thus be a negligible additive
cost (we also expect lower constant factors). But the most straightforward
translation of this second quantized algorithm into first quantization gives
an approach scaling as $\widetilde{\cal O}(\eta^{2}M)$, which would require
taking $M=\widetilde{\cal
O}(N^{2/3}/(\eta^{2/3}\epsilon)+\eta^{2/3}N^{1/3}/\epsilon)$ in order to
remain an additive negligible cost. This is likely fine for some applications
but future work should investigate whether this scaling can be improved.
Likewise, the topic of state preparation beyond single Slater determinants
(even in second quantization) has often been ignored by papers providing
error-corrected cost estimates for chemistry. However, when complex initial
states are required (e.g., for challenging systems like FeMoCo), preparing
these states may introduce a non-negligible cost that should be taken into
account more carefully.

  * •

The qubitization algorithm of this work can be used to perform time-evolution,
as opposed to used for preparing eigenstates for phase estimation, with the
very minimal additional overhead (see (20)) of quantum signal processing Low
and Chuang (2019). We expect that time-evolution of the non-Born-Oppenheimer
Hamiltonian would be especially useful for computing chemical reaction rates
Lidar and Wang (1999); Kassal  _et al._ (2008). There would also need to be
minor modifications to the block encoding algorithms of this paper to simulate
the nuclear kinetic operators and explicit nuclei interaction terms appearing
in (9). Future work should analyze the constant factors of such modifications
in the context of a specific chemical dynamics problem and estimate the
resources required to solve that problem within fault-tolerance. This would
constitute a fundamentally new type of application assessed for viability on
an error-corrected quantum processor. Furthermore, we note that our
interaction picture technique is trivially compatible with state preparation
methods based on adiabatic state preparation Babbush  _et al._ (2014) as we
can easily use that approach to perform time-dependent Hamiltonian evolutions.
Such methods might prove important for preparing complex initial states when
single Slater determinants fail to have sufficient overlap with the
eigenstates of interest. However, more work is required to figure out how to
best use the qubitization based approach in the context of adiabatic state
preparation.

  * •

Finally, in Appendix K we asymptotically describe how one can instead simulate
a real space version of these algorithms using the grid-like representation
employed by Kassal _et al._ Kassal  _et al._ (2008), but employing more modern
quantum simulation methods. We introduce algorithms for block encoding $U$ and
$V$ in this representation that lead to qubitization and interaction picture
algorithms with asymptotically the same complexity as the main algorithms
explored in this work. However, the constant factors are different and the
non-periodic nature of that representation may have advantages when simulating
non-periodic systems. Thus, we propose that future work should more thoroughly
explore this simulation strategy and compare the constant factors associated
with it to those of the techniques in this paper.

In conclusion, there are still many aspects of these algorithms (and their
applications) to explore. Nonetheless, this first study of fault-tolerant
overheads of first quantized quantum simulations of chemistry suggests many
potential advantages over their second quantized counterparts. While Gaussian
orbital based algorithms in second quantization will continue to compete with
first quantized algorithms when simulating small non-periodic molecules, we
believe that the longer term future of quantum computing for chemistry will
focus much more on first quantized simulations. We expect this is the case due
to the ease of removing the Born-Oppenheimer approximation, the fact that one
can more efficiently refine the accuracy of simulations in first quantization,
the fact that one does not require complicated and costly classical
precomputation of active spaces and molecular orbitals, and the fact that, on
a whole, these methods scale better than their second quantized counterparts,
while still having reasonable constant factors (as our work has shown).

### Acknowledgements

The authors thank Garnet Kin-Lic Chan, Matthias Degroote, Craig Gidney, Cody
Jones, Joonho Lee, Jarrod McClean, Yuval Sanders, Norm Tubman and the quantum
computing team at BASF for helpful discussions. DWB worked on this project
under a sponsored research agreement with Google Quantum AI. DWB is also
supported by Australian Research Council Discovery Projects DP190102633 and
DP210101367. NW was funded by a grant from Google Quantum AI as well as the US
Department of Energy, Office of Science, National Quantum Information Science
Research Centers, Co-Design Center for Quantum Advantage under contract number
DE-SC0012704 which supported his work on the asymptotic analysis of the
algorithms. YS was partly supported by the National Science Foundation RAISE-
TAQS 1839204. The Institute for Quantum Information and Matter is an NSF
Physics Frontiers Center PHY-1733907.

## References

  * Babbush  _et al._ (2019) Ryan Babbush, Dominic W Berry, Jarrod R McClean, and Hartmut Neven, “Quantum Simulation of Chemistry with Sublinear Scaling in Basis Size,” npj Quantum Information 5, 92 (2019).
  * Low and Chuang (2019) Guang Hao Low and Isaac L Chuang, “Hamiltonian Simulation by Qubitization,” Quantum 3, 163 (2019).
  * Low and Wiebe (2018) Guang Hao Low and Nathan Wiebe, “Hamiltonian Simulation in the Interaction Picture,” arXiv:1805.00675 (2018).
  * Feynman (1982) Richard P Feynman, “Simulating physics with computers,” International Journal of Theoretical Physics 21, 467–488 (1982).
  * Lloyd (1996) Seth Lloyd, “Universal Quantum Simulators,” Science 273, 1073–1078 (1996).
  * Aspuru-Guzik  _et al._ (2005) Alan Aspuru-Guzik, Anthony D Dutoi, Peter J Love, and Martin Head-Gordon, “Simulated Quantum Computation of Molecular Energies,” Science 309, 1704 (2005).
  * Berry  _et al._ (2015) Dominic W Berry, Andrew M Childs, Richard Cleve, Robin Kothari, and Rolando D Somma, “Simulating Hamiltonian Dynamics with a Truncated Taylor Series,” Physical Review Letters 114, 90502 (2015).
  * Wiebe and Granade (2015) Nathan Wiebe and Christopher E Granade, “Efficient Bayesian Phase Estimation,” arXiv:1508.00869 (2015).
  * Poulin  _et al._ (2017) David Poulin, Alexei Y Kitaev, Damian Steiger, Matthew Hastings, and Matthias Troyer, “Fast Quantum Algorithm for Spectral Properties,” Physical Review Letters 121, 010501 (2017).
  * Haah  _et al._ (2018) Jeongwan Haah, Matthew B. Hastings, Robin Kothari, and Guang Hao Low, “Quantum algorithm for simulating real time evolution of lattice Hamiltonians,” SIAM Journal on Computing , FOCS18–250–FOCS18–284 (2018).
  * Campbell (2019) Earl Campbell, “Random Compiler for Fast Hamiltonian Simulation,” Physical Review Letters 123, 070503 (2019).
  * Childs  _et al._ (2021) Andrew M Childs, Yuan Su, Minh C Tran, Nathan Wiebe, and Shuchen Zhu, “Theory of Trotter Error with Commutator Scaling,” Physical Review X 11, 011020 (2021).
  * Motta  _et al._ (2020a) Mario Motta, Chong Sun, Adrian T K Tan, Matthew J O’Rourke, Erika Ye, Austin J Minnich, Fernando G S L Brandao, and Garnet Kin-Lic Chan, “Determining eigenstates and thermal states on a quantum computer using quantum imaginary time evolution,” Nature Physics 16, 205–210 (2020a).
  * Şahinoğlu and Somma (2020) Burak Şahinoğlu and Rolando D Somma, “Hamiltonian simulation in the low energy subspace,” arXiv:2006.02660 (2020).
  * Tran  _et al._ (2021) Minh C Tran, Yuan Su, Daniel Carney, and Jacob M Taylor, “Faster digital quantum simulation by symmetry protection,” PRX Quantum 2, 010323 (2021).
  * Whitfield  _et al._ (2011) James D Whitfield, Jacob Biamonte, and Alan Aspuru-Guzik, “Simulation of electronic structure Hamiltonians using quantum computers,” Molecular Physics 109, 735–750 (2011).
  * Wecker  _et al._ (2014) David Wecker, Bela Bauer, Bryan K Clark, Matthew B Hastings, and Matthias Troyer, “Gate-count estimates for performing quantum chemistry on small quantum computers,” Physical Review A 90, 022305 (2014).
  * Hastings  _et al._ (2015) Matthew B Hastings, Dave Wecker, Bela Bauer, and Matthias Troyer, “Improving Quantum Algorithms for Quantum Chemistry,” Quantum Information & Computation 15, 1–21 (2015).
  * Poulin  _et al._ (2015) David Poulin, M B Hastings, Dave Wecker, Nathan Wiebe, Andrew C Doherty, and Matthias Troyer, “The Trotter Step Size Required for Accurate Quantum Simulation of Quantum Chemistry,” Quantum Information & Computation 15, 361–384 (2015).
  * Babbush  _et al._ (2016) Ryan Babbush, Dominic W Berry, Ian D Kivlichan, Annie Y Wei, Peter J Love, and Alan Aspuru-Guzik, “Exponentially More Precise Quantum Simulation of Fermions in Second Quantization,” New Journal of Physics 18, 33032 (2016).
  * Babbush  _et al._ (2018a) Ryan Babbush, Craig Gidney, Dominic W Berry, Nathan Wiebe, Jarrod McClean, Alexandru Paler, Austin Fowler, and Hartmut Neven, “Encoding electronic spectra in quantum circuits with linear T complexity,” Physical Review X 8, 041015 (2018a).
  * Berry  _et al._ (2018) Dominic W Berry, Maria Kieferová, Artur Scherer, Yuval R Sanders, Guang Hao Low, Nathan Wiebe, Craig Gidney, and Ryan Babbush, “Improved Techniques for Preparing Eigenstates of Fermionic Hamiltonians,” npj Quantum Information 4, 22 (2018).
  * Berry  _et al._ (2019) Dominic W Berry, Craig Gidney, Mario Motta, Jarrod McClean, and Ryan Babbush, “Qubitization of Arbitrary Basis Quantum Chemistry Leveraging Sparsity and Low Rank Factorization,” Quantum 3, 208 (2019).
  * Kivlichan  _et al._ (2018) Ian D Kivlichan, Jarrod McClean, Nathan Wiebe, Craig Gidney, Alan Aspuru-Guzik, Garnet Kin-Lic Chan, and Ryan Babbush, “Quantum Simulation of Electronic Structure with Linear Depth and Connectivity,” Physical Review Letters 120, 110501 (2018).
  * Kivlichan  _et al._ (2020) Ian D Kivlichan, Craig Gidney, Dominic W Berry, Nathan Wiebe, Jarrod McClean, Wei Sun, Zhang Jiang, Nicholas Rubin, Austin Fowler, Alán Aspuru-Guzik, Hartmut Neven, and Ryan Babbush, “Improved Fault-Tolerant Quantum Simulation of Condensed-Phase Correlated Electrons via Trotterization,” Quantum 4, 296 (2020).