CROSSREFERENCE TO RELATED APPLICATIONS
This application claims the benefit of U.S. Provisional Application No. 62/656,517 entitled “HAMILTONIAN SIMULATION IN THE INTERACTION PICTURE” and filed on Apr. 12, 2018, which is hereby incorporated herein in its entirety.
SUMMARY
In this disclosure, quantum algorithms are presented for simulating Hamiltonian timeevolution e^{−i(A+B)t }in the interaction picture of quantum mechanics on a quantum computer. The interaction picture is a known analytical tool for separating dynamical effects due to trivial freeevolution A from those due to interactions B. This is especially useful when the energyscale of the trivial component is dominant, but of little interest. Whereas stateofart simulation algorithms scale with the energy ∥A∥+∥B∥ of the full Hamiltonian, embodiments of the disclosed approach generally scale with only the lowenergy component ∥B∥. Applied to simulating a periodic Nsite Hubbard model with arbitrary longranged densitydensity interactions, the disclosed gate complexity (N^{2}t) is a quadratic improvement in N over prior art for electronic structure simulation in the planewave basis. More generally in the abstract query model, diagonallydominant Hamiltonians with sparsity d can be simulated with cost independent of the diagonal component. The various contributions disclosed herein include a rigorous analysis of a lowspace overhead simulation algorithm based on the Dyson series for general timedependent Hamiltonians, which also enables a quadratic improvement in sparsity for timedependent simulation with (td) queries.
In some embodiments of the disclosed technology, a quantum computer is configured to simulate a quantum system, and a Hamiltonian in the simulation is represented in the interaction picture. The simulation of the quantum system is then performed using the quantum computer. In certain implementations, the simulation of the quantum system is a subroutine that is repeated two or more times. In some implementations, the simulation is performed using linear combinations of unitaries. For instance, in some cases, the simulation uses linear combinations of unitaries performed on a diagonally dominant matrix (e.g., using linear combinations of unitaries performed on the diagonally dominant components of the diagonally dominant matrix). In certain implementations, the quantum system is modelled by a Hubbard model. In particular implementations, the quantum system describes a physical chemical system or molecule. In some implementations, the Hamiltonian is sparse and the simulation uses the a state of an auxiliary qubit to encode the matrix elements of the Hamiltonian instead of using graph decomposition techniques. In certain implementations, the simulation is performed by compressing ancillas for quantum simulation of a timedependent Hamiltonian.
In further embodiments, a quantum algorithm is implemented on a quantum computer for simulating a general sparse timedependent quantum system. In this embodiment, the quantum algorithm does not use graph decomposition techniques. In certain implementations, a quantum Hamiltonian used in the simulation is represented in the interaction picture. In some implementations, the simulation uses linear combinations of unitaries performed on a diagonally dominant matrix. In particular implementations, an interaction Hamiltonian is chosen to be a diagonal matrix. In certain implementations, the quantum system is modelled by a Hubbard model, a physical chemical system, or a molecule. In further implementations, the simulation includes compressing ancillas use to index the time for the simulation.
The foregoing and other objects, features, and advantages of the disclosed technology will become more apparent from the following detailed description, which proceeds with reference to the accompanying figures.
BRIEF DESCRIPTION OF THE DRAWINGS
FIG. 1 is a schematic block diagram showing a circuit that implements HAM with normalization constant α=Σ_{j=1}^{l }α_{j }using example unitary oracles.
FIG. 2 shows a quantum circuit representations of the components implementing the truncated Dyson series algorithm TDS.
FIG. 3 is a quantum circuit representations of the gadget V for applying a sequence of probabilistic operators H_{k }. . . H_{2}H_{1}, encoded in (0_{α}⊗I_{s})U_{k}(0_{α}⊗I_{s})=H_{k}, controlled on number state k_{b}, k∈{0,1, . . . , K}.
FIG. 4 depicts the quantum circuit representation of an example implementation of HAM_{K }from Eq. (75) using K queries to controlledHAM, a single step of the truncated Taylor series algorithm before oblivious amplitude amplification, and a single step of timeevolution by the truncated Taylor series algorithm from Eq. (78).
FIG. 5 is a schematic block diagram showing quantum circuit representations of the unitary DYS_{K }that encodes timeordered products of Hamiltonians, using fewer ancilla than the construction in FIG. 2 by applying the compression gadget of FIG. 4.
FIG. 6 is an example method for a timedependent simulation algorithm as disclosed herein.
FIG. 7 is an example method for an example interaction picture simulation method as disclosed herein.
FIG. 8 is an example method simulating chemistry and Hubbard models as disclosed herein.
FIG. 9 is an example method of a compression method as disclosed herein.
FIG. 10 illustrates a generalized example of a suitable classical computing environment in which aspects of the described embodiments can be implemented.
FIG. 11 illustrates an example of a possible network topology (e.g., a clientserver network) for implementing a system according to the disclosed technology.
FIG. 12 illustrates another example of a possible network topology (e.g., a distributed computing environment) for implementing a system according to the disclosed technology.
FIG. 13 illustrates an exemplary system for implementing the disclosed technology.
FIG. 14 illustrates an example method for performing a quantum simulation in accordance with an embodiment of the disclosed technology.
FIG. 15 illustrates a further example method for performing a quantum simulation in accordance with an embodiment of the disclosed technology.
DETAILED DESCRIPTION
I. General Considerations
As used in this application, the singular forms “a,” “an,” and “the” include the plural forms unless the context clearly dictates otherwise. Additionally, the term “includes” means “comprises.” Further, the term “coupled” does not exclude the presence of intermediate elements between the coupled items. Further, as used herein, the term “and/or” means any one item or combination of any items in the phrase.
Although the operations of some of the disclosed methods are described in a particular, sequential order for convenient presentation, it should be understood that this manner of description encompasses rearrangement, unless a particular ordering is required by specific language set forth below. For example, operations described sequentially may in some cases be rearranged or performed concurrently. Moreover, for the sake of simplicity, the attached figures may not show the various ways in which the disclosed systems, methods, and apparatus can be used in conjunction with other systems, methods, and apparatus. Additionally, the description sometimes uses terms like “produce” and “provide” to describe the disclosed methods. These terms are highlevel abstractions of the actual operations that are performed. The actual operations that correspond to these terms will vary depending on the particular implementation and are readily discernible by one of ordinary skill in the art.
II. Introduction
Manybody dynamics are often described by a Hamiltonian H=A+B with two noncommuting parts: a wellunderstood but trivial freeparticle theory A, and a more interesting component B describing interactions. One often finds a large separation of energy scales ∥A∥»∥B∥ in nature. For instance, the ultraviolet cutoff of a free scalar field is usually assumed to be significantly higher than the energy scale of its interactions. In practical applications, the energy contribution of electronelectron correlations in molecules and materials are typically larger than chemical accuracy, at 4 kJ mol^{−1}, but still much smaller than the HartreeFock component.
Dynamical effects arising from the interaction are often the subject of interest. Whereas the quantum statevector φ_{I}(t) evolves under the full Hamiltonian by the Schrödinger equation i∂_{t}φ(t)=(A+B)φ(t), the interaction picture moves to the rotating frame of the trivial dynamics. This frame allows one to focus on effects of the interaction, and is particularly fruitful when interactions are perturbative corrections to the freetheory. For instance, an analytic evaluation of the timeordered propagator of the interaction picture Hamiltonian H_{I}(t)=e^{iAt}Be^{−iAt }is possible by perturbative expansions based on Green'"'"'s functions and Feynman diagrams. However, such techniques fail for interactions that are nonperturbative, even if they are weak relative to the free theory. In such situations, one resorts to numerical simulations.
Numerical simulations of quantum dynamics scale exponentially with particle number, so the interaction picture provides little advantage here. One then relies on nonperturbative numerical approximations such as densityfunctional theory or Monte Carlo pathintegral evaluation. Though highly successful in elucidating qualitative behavior, obtaining accurate quantitative predictions to high precision is difficult. For instance, achieving chemical accuracy in abinitio quantum chemistry calculations of reaction rates often only possible by exact diagonalization of the Hamiltonian. This has an exponentially growing cost ^{(n) }with the system size n. (The standard big notation defines f(n)=(g(n)) for positive functions f(n),g(n)>0 as the existence of absolute constants α>0,b>0 such that for any n>α, f(n)≤bg(n). One can also use f(n)=(g(n)) when f(n)≤bg(n) polylog(g(n)), and f(n)=Ω(g(n)) when f(n)≥bg(n), and f(n)=⊖(g(n)) when both f(n)=(g(n)) and f(n)=Ω(g(n)) are true.)
One approach that is tractable and nonperturbative is simulating quantum dynamics on a quantum computer. Given a (poly(n))sized description of a timeindependent Hamiltonian H, one devises a sequence of (poly(n)) primitive quantum gates that approximates the unitary timeevolution operator e^{−iHt }with error ϵ. (Primitive gate costs are defined as the number N of arbitrary single and twoqubit gates. In faulttolerant architectures, one multiples by an additional overhead (log (N/ϵ)) for approximating N such gates to total error ϵ using a universal discrete gate set, e.g. CLIFFORD+T.) This is then applied to the system state φ(0) encoded in (poly(n)) logical quantum bits to obtain the timeevolved state φ(t)=e^{−iHt}φ(0). Algorithms for Hamiltonian simulation have progressed rapidly culminating recently in schemes with optimal query complexity (nt+log (1/ϵ)) with respect to all parameters for generic sparse Hamiltonians, and algorithms for geometricallylocal exponentiallydecaying interactions, with gate complexity (nt polylog(nt/ϵ)) that is optimal up to polylogarithmic factors.
It would be desirable if the cost of simulating timeevolution scaled with energy of the smaller interaction term, say (t∥B∥). As the interaction picture Hamiltonian has norm ∥H_{I}(t)∥=∥B∥, this possibility is highly suggestive. Nevertheless, stateofart performs simulation in the Schrödinger picture, where the gate cost of simulating e^{−iHt }generally scales at least like Ω(t(∥A∥+∥B∥)log(1/ϵ)), even if the trivial dynamics e^{−iAt }of term A may be evaluated in closedform. Even with specialized techniques that scale with a geometric mean (^{(k)}t∥A∥(t∥B∥/ϵ)^{1/k}), ∀k∈2^{+}. the uninteresting component (t∥A∥) is the dominant prefactor.
Unfortunately, simulation in the interaction picture is difficult. There, evolution by a timeindependent Hamiltonian H is transformed in the rotating frame φ_{I}(t)=e^{iAt}φ(t) to evolution by a timedependent Hamiltonian H_{I}(t)=e^{iAt}Be^{−iAt}. This follows from an elementary manipulation of the Schrödinger equation:
i∂_{t}φ(t)=(A+B)φ(t)→i∂_{t}φ_{I}(t)=e^{iAt}Be^{−iAt}φ_{I}(t). (1)
Implementing the timeordered propagator [exp(−i∫_{0}^{t}H(s)ds)] that solves Eq. (1) on a quantum computer requires timedependent simulation algorithms. These are generally more complicated than timeindependent algorithms, and exhibit different cost tradeoffs that do not appear favorable. For instance, an orderk timedependent TrotterSuzuki product formula has a cost that scales with the rate of change of H(t) like (^{(k)}(t∧)^{1+1/(2k)}), where ∧=max_{s}∥{dot over (H)}(s)∥^{1/2}=(∥[A,B]∥^{1/2}). More advanced techniques based on compressed fractional queries appear to scale better like ˜
$t\ue89e\uf605B\uf606\ue89e\frac{\mathrm{log}\ue8a0\left(\Lambda \ue89e\phantom{\rule{0.3em}{0.3ex}}\ue89et/\u03f5\right)}{\mathrm{log}\ue89e\phantom{\rule{0.3em}{0.3ex}}\ue89e\mathrm{log}\ue8a0\left(\Lambda \ue89e\phantom{\rule{0.3em}{0.3ex}}\ue89et/\u03f5\right)}$
but in terms of queries to a unitary oracle that obscures the gate complexity as it expresses Hamiltonian matrix elements at different times in binary, and may be difficult to implement in practice. One proposed technique directly implements a truncated Dyson series of a the timeordered propagator and argues, though without proof, a similar scaling in terms of queries to a different type of oracle.
Here, it is shown that simulation in the interaction picture can substantially improve the efficiency of simulation. In Section III, the general timedependent simulation algorithm is completed by a truncated Dyson series by providing a rigorous analysis of the approximation and explicit circuit constructions, with improvements in gate and space complexity over previously expected costs. In Section IV, situations are identified where the gate complexity of implementing these queries scale with the interaction strength (∥B∥), and not the larger uninteresting component (∥A∥). Such are the cases where simulation in the interaction picture is advantageous. In Section V, the potential of interactionpicture simulation by an electronic structure application in the planewave basis is discussed. Further, the cost of simulating the timeevolution of N spinorbitals subject to longrange electronelectron interactions to (N^{2}t) gates is discussed, which is close to a quadratic improvement over prior art of (N^{11/3}t). In Section VI, a complexity theoretic perspective of the disclosure is provided by considering the abstract problem of simulating timedependent sparse Hamiltonians in the standard query model. A quadratic improvement in sparsity scaling is shown, and optimized algorithms for simulating diagonally dominant Hamiltonians are provided. A more detailed summary of main results in each section is as follows.
In Section III (“Timedependent Hamiltonian simulation by a truncated Dyson series”), an example algorithm for a general timedependent simulation algorithm Theorem 1 is provided. The example algorithm is based on approximating a truncation and discretization of the Dyson series for general timedependent Hamiltonians H(t) characterized by spectralnorm α≥max_{t}∥H(t)∥ and average rateofchange
$\u3008\uf605\stackrel{.}{H}\uf606\u3009=\frac{1}{t}\ue89e{\int}_{0}^{t}\ue89e\uf605\stackrel{.}{H}\ue8a0\left(s\right)\uf606\ue89e\mathrm{ds}.$
A rigorous analysis of its performance is provided in together with explicit circuit constructions. Bounds on the approximation error is evaluated in Section III A, which is used to obtain the cost of simulating the timeordered evolution operator. This cost is found to be
$\ue52e\ue8a0\left(\alpha \ue89e\phantom{\rule{0.3em}{0.3ex}}\ue89et\ue89e\phantom{\rule{0.3em}{0.3ex}}\ue89e\frac{\mathrm{log}\ue89e\phantom{\rule{0.3em}{0.3ex}}\ue89e\left(\alpha \ue89e\phantom{\rule{0.3em}{0.3ex}}\ue89et/\u03f5\right)}{\mathrm{log}\ue89e\phantom{\rule{0.3em}{0.3ex}}\ue89e\mathrm{log}\ue8a0\left(\alpha \ue89e\phantom{\rule{0.3em}{0.3ex}}\ue89et/\varepsilon \right)}\right)$
queries to a certain unitary oracle HAMT that encodes H(t) evaluated
$M=\ue52e\left(\frac{1}{\u03f5}\ue89e\frac{\u3008\uf605\stackrel{.}{H}\uf606\u3009}{{\alpha}^{3}}\right)$
points. Example implementations as disclosed in Section III B, characterized by Theorem 2, completes the scheme proposed by D. W. Berry, A. M, Childs, and R. Kothari, “Hamiltonian simulation with nearly optimal dependence on all parameters,” in Proceedings of the 56th Annual IEEE Symposium on Foundations of Computer Science (FOCS), pages 792809 (October 2015), with some improvements in gate complexity. One obtains scaling with the average rateofchange (log(∥H∥)), instead of the worst case (log(max_{t}∥H(t)∥)). This allows one to simulate Hamiltonians with arbitrary elementwise slewrates, so long as ∥H(t)∥ remains bounded. One can then devise a compression gadget in Section III C that reduces the qubit overhead complexity of this implementation to (log(M)), compared to
$\ue52e\ue8a0\left(\mathrm{log}\ue8a0\left(M\right)\ue89e\frac{\mathrm{log}\ue8a0\left(\alpha \ue89e\phantom{\rule{0.3em}{0.3ex}}\ue89et/\u03f5\right)}{\mathrm{log}\ue89e\phantom{\rule{0.3em}{0.3ex}}\ue89e\mathrm{log}\ue89e\phantom{\rule{0.3em}{0.3ex}}\ue89e\left(\alpha \ue89e\phantom{\rule{0.3em}{0.3ex}}\ue89et/\u03f5\right)}\right)$
of the original. This gadget is of independent interest and also applicable to the timeindependent truncated Taylor series algorithm, discussed in Section XA for completeness.
In Section IV (“Accelerated interaction picture simulation”), one can apply this truncated Dyson series algorithm to simulate timeevolution by e^{−i(A+B)t }in the interaction picture. One can evaluate the gate complexity in Section IV A for one possible construction of these queries HAMT for the interaction picture Hamiltonian H_{I}(t)=e^{iAt}Be^{−iAt}. This leads to the simulation algorithm Theorem 3 for timeevolution by e^{−i(A+B)t }with gate complexity
((C_{B}+C_{e}_{iA/α}_{B})α_{B}tpolylog(∥A∥,α_{B},t,ϵ)), (2)
where C_{B }is the gate complexity a unitary oracle O_{B }such that
$(\u30080\ue89e{}_{a}\ue89e\otimes {I}_{s})\ue89e{O}_{B}\ue8a0\left({\uf6030\u3009}_{a}\otimes {I}_{s}\right)=\frac{B}{{\alpha}_{B}},$
and C_{e}_{iA/α}_{B }is the cost of simulating timeevolution by A for time (α_{B}^{−1}). This may be compared to stateofart Schrödinger picture simulation algorithms for timeindependent Hamiltonians, which have gate complexity
((C_{B}+C_{A})(α_{B}+α_{A})tpolylog(∥A∥,∥B∥,t,ϵ)), (3)
where C_{A }is analogous to C_{B }but for the Hamiltonian term A. The disclosed result Eq. (2) is then advantageous in cases where ∥A∥»∥B∥ and C_{e}_{iA/α}=(C_{B}). There, the gate complexity scales with the interaction strength ∥B∥, and not the larger uninteresting component ∥A∥, up to polylogarithmic factors.
In Section V (“Application to the Hubbard model with longranged interactions”),
One can demonstrate a concrete advantage for Hamiltonian simulation in the interaction picture over timeindependent simulation algorithms. One can consider a general Hubbard periodic model on N lattice sites in an arbitrary number of dimensions, which has Hamiltonian
$\begin{array}{cc}H=\sum _{\overrightarrow{x},\overrightarrow{y},\sigma}\ue89e\phantom{\rule{0.3em}{0.3ex}}\ue89eT\ue8a0\left(\overrightarrow{x}\overrightarrow{y}\right)\ue89e{a}_{\overrightarrow{x}\ue89e\sigma}^{\u2020}\ue89e{a}_{\overrightarrow{y},\sigma}+\sum _{\overrightarrow{x},\sigma}\ue89eU\ue8a0\left(\overrightarrow{x},\sigma \right)\ue89e{n}_{\overrightarrow{x}\ue89e\sigma}+\sum _{\left(\overrightarrow{x},\sigma \right)\ne \left(\overrightarrow{y},{\sigma}^{\prime}\right)}\ue89eV\ue8a0\left(\overrightarrow{x}\overrightarrow{y}\right)\ue89e{n}_{\overrightarrow{x}\ue89e\sigma}\ue89e{n}_{\overrightarrow{y}\ue89e{\sigma}^{\prime}},& \left(4\right)\end{array}$
where α_{{right arrow over (x)}σ}^{†}, α_{{right arrow over (x)}σ }are Fermionic creation and annihilation operators indexed by lattice position {right arrow over (x)} and spin σ, and n_{{right arrow over (x)}σ }is a number operator. Note that one can allow for arbitrary translationally invariant kinetic T(⋅) hopping terms and longrange densitydensity V(⋅) interactions, subject to arbitrary local disorder U(⋅). Provided that the kinetic term is extensive, embodiments of the disclosed interactionpicture algorithm have gate complexity (N^{2}t), independent of U and V up to polylogarithmic factors. This model generalizes electronic structure simulations in the planewave basis, considered in Section V A, where one obtains almost a quadratic improvement over the prior art of (N^{11/3}t) gates for this special case. This is complementary to recent work by Jeongwan Haah, Matthew B Hastings, Robin Kothari, and Guang Hao Low, Quan tum algorithm for simulating real time evolution of lattice Hamiltonians. arXiv preprint arXiv:1801.03922, 2018, which achieves (Nt) scaling, but under the assumption of shortrange exponentially decaying coefficients.
In Section VI (“Application to sparse Hamiltonian simulation”), the application to the Hubbard model with longranged interactions is explained.
Consider next a complexitytheoretic generalization of the technology that has been developed for timedependent simulation and simulation in the interaction picture. The standard query model for blackbox dsparse Hamiltonian simulation assumes access to a unitary oracle that provides the positions and values of nonzero entries of the Hamiltonian. Each row has at most d nonzero entries, and the maximum absolute value of any entry is ∥H∥_{max}. This infortmation is also provided as a function of a time index for timedependent Hamiltonians, In Section VI A, one considers this timedependent case and show with Theorem 4 that the timeordered evolution operator may be simulated using
$\ue52e\ue8a0\left(\mathrm{td}\ue89e{\uf605H\uf606}_{\mathrm{max}}\ue89e\frac{\mathrm{log}\ue8a0\left(\mathrm{td}\ue89e{\uf605H\uf606}_{\mathrm{max}}/\u03f5\right)}{\mathrm{log}\ue89e\phantom{\rule{0.3em}{0.3ex}}\ue89e\mathrm{log}\ue8a0\left(\mathrm{td}\ue89e{\uf605H\uf606}_{\mathrm{max}}/\u03f5\right)}\right)$
queries. Though linear scaling with respect to d is wellknown in the timeindependent case, this is a quadratic improvement in sparsity scaling over prior art for the timedependent case. An analogous treatment for simulating sparse timeindependent Hamiltonian in the interaction picture in Section VI B, Theorem 5 has query complexity that scales only with ∥H_{od}∥_{max}, which is the maximum absolute value of any offdiagonal entry. This improved cost
$\ue52e\ue8a0\left(\mathrm{td}\ue89e{\uf605{H}_{\mathrm{od}}\uf606}_{\mathrm{max}}\ue89e\frac{\mathrm{log}\ue8a0\left(\mathrm{td}\ue89e{\uf605{H}_{\mathrm{od}}\uf606}_{\mathrm{max}}/\u03f5\right)}{\mathrm{log}\ue89e\phantom{\rule{0.3em}{0.3ex}}\ue89e\mathrm{log}\ue8a0\left(\mathrm{td}\ue89e{\uf605{H}_{\mathrm{od}}\uf606}_{\mathrm{max}}/\u03f5\right)}\right),$
is particularly advantageous for the common situation of diagonally dominant Hamiltonians.
III TimeDependent Hamiltonian Simulation by A Truncated Dyson Series
The timedependent simulation technique applied throughout this work is motivated by the Dyson series which solves the timedependent Schrödinger equation. The Dyson series is a powerseries expansion of the unitary propagator that involves integrals of products of H(t) at different times. The technique implements a truncation of this powerseries, where integrals in each term is discretized in time.
This technology is presented with a rigorous analysis of the truncation and discretization error in Section III A, followed by an explicit circuit construction that is in the spirit of earlier work in Section III B. Two improvements are of note. First, the gate complexity of the algorithm scales with the average rateofchange (log(∥H∥)), instead of the worst case (log (max_{t}∥H(t)∥). This allows one to simulate Hamiltonians with arbitrary elementwise slewrates, so long as ∥H(t)∥ remains bounded. Second, in Section III C, an alternate implementation of the technique is introduced that reduces the space overhead of the original proposal substantially, and has the same query and gate complexity.
The complexity of any simulation algorithm depends on assumptions made about how information on the timedependent Hamiltonian H(t) is accessed by the quantum computer. This section begins by defining an input model for the unitary quantum oracle HAMT that encodes H(t). The complexity of timedependent simulation is then expressed in terms of queries to HAMT (It is assumed that the query complexity to a controlledunitary blackbox is the same as that to the original blackbox. In general, this will not affect the gate complexity. Though there are often cleverer ways to implement an arbitrary controlledunitary, in the worstcase, all quantum gates may be replaced by their controlled versions with constant overhead.), any additional primitive gates required beyond HAMT, and the total number of qubits required including that for HAMT. In later sections where the timedependent simulation algorithm is applied to the interaction picture, this oracle is “opened up” and and the gate complexity of its implementation is discussed.
A common assumption in the timeindependent situation is that H=Σ_{j=1}^{l}α_{j}H_{j }is a sum of l local Hermitian terms with positive coefficients α_{j}, and that each term e^{−iH}^{j}^{t }may be implemented with (l) quantum gates for any t. In abstract models for sparse Hamiltonians, a binary representation of matrix entries H_{jk }and positions are returned in superposition by querying a certain unitary quantum oracle that is a natural generalization of classical Boolean circuits, Here, Hamiltonians H is considered encoded in a socalled standardform. For any timeindependent Hamiltonian H, it is assumed that there exists a unitary oracle HAM such that
(0_{a}⊗I_{s})HAM(0_{a}⊗I_{s})=H/α, (5)
where HAM acts jointly on registers α, s consisting of n_{α}, n_{s }qubits respectively, and α≥∥H∥ is a normalization constant that should be made small for best performance. In other words, postselected on the n_{α}qubit measurement outcome 0_{α }the Hamiltonian H is applied on any arbitrary state in the system register s. Use of this is justified as it generalizes a variety of different input models. Thus any simulation algorithm that queries Eq. (5) for this general case readily specializes to more structured inputs. As an example, H=Σ_{j=1}^{2}^{n}^{α}α_{j}U_{j }could be a linear combination of l=2^{n}^{α }unitaries. Then the circuit depicted in FIG. 1 implements HAM with normalization constant α=Σ_{j=1}^{l }α_{j }using the unitary oracles
$\begin{array}{cc}\mathrm{HAM}=\left({\mathrm{PREP}}^{\u2020}\otimes {I}_{s}\right)\xb7\mathrm{SEL}\xb7\left(\mathrm{PREP}\otimes {I}_{s}\right),\text{}\ue89e\mathrm{PREP}\ue89e\phantom{\rule{0.8em}{0.8ex}}\ue89e{\uf6030\u3009}_{a}=\sum _{j=1}^{i}\ue89e\phantom{\rule{0.3em}{0.3ex}}\ue89e\sqrt{\frac{{a}_{j}}{\alpha}}\ue89e{\uf603j\u3009}_{a},\text{}\ue89e\mathrm{SEL}=\sum _{j=1}^{i}\ue89e\phantom{\rule{0.3em}{0.3ex}}\ue89e\uf603j\u3009\ue89e{\u3008j\uf604}_{a}\otimes {U}_{j}.& \left(6\right)\end{array}$
These unitaries each cost (l) gates—PREP is implemented by arbitrary quantum state preparation on l dimensions, and SEL is implemented by binarytree control logic for l different inputs. Typically, the number of terms in the Hamiltonian scales like l=(poly(n_{s})), so n_{α}=(log(l)). However, the number of ancilla qubits can also be polynomial in n_{s}, such as in the simulation of sparse Hamiltonians where n_{α}≥n_{s}+2, or where a spacetime complexity tradeoff is possible.
FIG. 1 shows a quantum circuit representation of (left) an oracle HAM from Eq. (5) encoding a timeindependent Hamiltonian, (center) an oracle HAMT from Eq. (7) encoding a timedependent Hamiltonian, and (right) an example implementation of HAM from with a linearcombination of unitaries from Eq. (6). Bold horizontal lines depict registers that in general comprise of multiple qubits. Vertical lines connecting boxes depict unitaries that act jointly on all registers covered by the boxes. A small square box marked by T indicates control by a time index.
A number of techniques may be applied to approximate the unitary timeevolution operator e^{−iHt }to error ϵ by querying HAM. The combined Qubitization and quantum signal processing approach is one such algorithm with query complexity (αt+log(1/ϵ)) that is optimal for sparse Hamiltonians. Alternatively, one could apply a truncated Taylor series of e^{−iHt}, which has a slightly worse query and space complexity, but benefits from being conceptually simpler. As the Dyson series algorithm is timedependent generalization of the truncated Taylor series algorithm, a review and some useful background is provided before in Section X A for completeness.
A natural timedependent generalization of Eq. (6) is the unitary oracle HAMT encoding the Hamiltonian H(s) defined over time s ∈[0,t], with t>0. The continuous parameter s is discretized into an integer number of M>0 time bins which are indexed by m ∈ {0, 1, . . . , M−1}. Here, it is assumed that H(s) is encoded in an analogous format
$\begin{array}{cc}\left({\u30080\uf604}_{a}\otimes {I}_{s}\right)\ue89e\mathrm{HAM}\ue89e\text{}\ue89eT\ue8a0\left({\uf6030\u3009}_{a}\otimes {I}_{s}\right)=\sum _{m=0}^{M1}\ue89e\phantom{\rule{0.3em}{0.3ex}}\ue89e\uf603m\u3009\ue89e{\u3008m\uf604}_{d}\otimes \frac{H\ue8a0\left(m\ue89e\phantom{\rule{0.3em}{0.3ex}}\ue89e\Delta \right)}{\alpha}.& \left(7\right)\end{array}$
Given blackbox access to HAMT, details of a quantum algorithm that approximates the Dyson series Eq. (8) are provided. The overall algorithm is efficient so long as HAMT has an efficient implementation in terms of primitive quantum gates. More formally, the performance of the algorithm is captured by the following theorem, proven in the remainder of this section.
Theorem 1 (Hamiltonian simulation by a truncated Dyson series). Let a timedependent Hamiltonian H(s) be characterized by spectral norm α≥max_{s }∥H(s)∥, and average rateofchange
$\u3008\uf605\stackrel{.}{H}\uf606\u3009=\frac{1}{t}\ue89e{\int}_{o}^{t}\ue89e\uf605\frac{\mathrm{dH}\ue8a0\left(s\right)}{\mathrm{ds}}\ue89e\phantom{\rule{0.2em}{0.2ex}}\uf606$
ds. This Hamiltonian is defined for s∈[0, t], where αt≤½, and is encoded in an oracle HAMT acting on n_{s}+n_{α}+n_{d }qubits, as per Eq. (7), evaluated at
$M=\ue52e\left(\frac{{t}^{2}}{\u03f5}\ue89e\left(\frac{\u3008\uf605\stackrel{.}{H}\uf606\u3009}{\alpha}+\frac{{\mathrm{max}}_{s}\ue89e{\uf605H\ue8a0\left(s\right)\uf606}^{2}}{{\alpha}^{2}}\right)\right)$
points, where n_{d}=(log (M)). Then the timeordered evolution operator [exp(−i∫_{0}^{t}H(s)ds)] may be approximated to error ϵ with the following cost.
$1.\ue89e\phantom{\rule{0.8em}{0.8ex}}\ue89e\mathrm{Queries}\ue89e\phantom{\rule{0.8em}{0.8ex}}\ue89e\mathrm{to}\ue89e\phantom{\rule{0.8em}{0.8ex}}\ue89e\mathrm{HAM}\ue89e\text{}\ue89eT\ue89e\text{:}\ue89e\phantom{\rule{0.8em}{0.8ex}}\ue89e\ue52e\ue8a0\left(\frac{\mathrm{log}\ue8a0\left(t/\u03f5\right)}{\mathrm{log}\ue89e\phantom{\rule{0.3em}{0.3ex}}\ue89e\mathrm{log}\ue8a0\left(t/\u03f5\right)}\right).\text{}\ue89e2.\ue89e\phantom{\rule{0.8em}{0.8ex}}\ue89e\mathrm{Qubits}\ue89e\text{:}\ue89e\phantom{\rule{0.8em}{0.8ex}}\ue89e{n}_{s}+\ue52e\ue8a0\left({n}_{a}+\mathrm{log}\ue8a0\left(M\right)\right).\text{}\ue89e3.\ue89e\phantom{\rule{0.8em}{0.8ex}}\ue89e\mathrm{Primitive}\ue89e\phantom{\rule{0.8em}{0.8ex}}\ue89e\mathrm{gates}\ue89e\text{:}\ue89e\phantom{\rule{0.8em}{0.8ex}}\ue89e\ue52e\ue8a0\left(\left({n}_{a}+\mathrm{log}\ue8a0\left(M\right)\right)\ue89e\frac{\mathrm{log}\ue8a0\left(1/\u03f5\right)}{\mathrm{log}\ue89e\phantom{\rule{0.3em}{0.3ex}}\ue89e\mathrm{log}\ue8a0\left(1/\u03f5\right)}\right).$
Note that, this algorithm performs timeevolution for a relatively short duration of t≤α^{−1}. Timeevolution for longer durations is simulated by breaking the timeordered evolution operator [exp(−i∫_{0}^{t}H(s)ds)] into L=(αt) timesteps 0=t_{0}<t_{1 }. . . <t_{L}=t of size (α^{−1}), and each step is implemented by the quantum circuit underlying Theorem 1 that queries a different HAMT now defined for a Hamiltonian H(s) for s∈[t_{j},t_{j−1}]. As the final error is Lϵ by a simple triangle inequality, one can rescale ϵ→ϵ/L. Thus simulation for the full duration to error ϵ has query complexity
$\ue52e\ue8a0\left(\alpha \ue89e\phantom{\rule{0.3em}{0.3ex}}\ue89et\ue89e\frac{\mathrm{log}\ue8a0\left(\alpha \ue89e\phantom{\rule{0.3em}{0.3ex}}\ue89et/\u03f5\right)}{\mathrm{log}\ue89e\phantom{\rule{0.3em}{0.3ex}}\ue89e\mathrm{log}\ue8a0\left(\alpha \ue89e\phantom{\rule{0.3em}{0.3ex}}\ue89et/\u03f5\right)}\right).$
A. Dyson Series Truncation Error
The choice of parameters in Theorem 1 depends on the scaling of error from a finite and discrete approximation of the Dyson series. Consider a unitary propagator U(t) defined to solve the timedependent Schrödinger equation i∂_{t}φ(t)=H(t)φ(t) by providing the state φ(t)=U(t)φ(0) given the initial state φ(0) at time t=0. For any t>0 and bounded ∥H(t)∥, U(t) has an absolutely convergent infinite expansion known as the Dyson series
U(t)=I−i∫_{0}^{t}H(t_{1})dt_{1}−∫_{t}_{2}^{t}∫_{0}^{t}^{2}H(t_{2})H(t_{1})dt_{1}dt_{2}+i∫_{t}_{3}^{t}∫_{t}_{2}^{t}^{3}∫_{0}^{t}^{2}H(t_{3})H(t_{2})H(t_{1})dt_{1}dt_{2}dt_{3}+ . . . . (8)
This may be compactly represented using the timeordering operator z,99 which sorts any sequence of k operators according to the times t_{j }of their evaluation, that is, [H(t_{k}) . . . H(t_{2})H(t_{1})]=H(t_{σ(k)}) . . . H(t_{σ(2)})H(t_{σ(1)}), where σ is a permutation such that t_{σ(1)}≤t_{σ(2)}≤ . . . ≤t_{σ(k)}. For instance, [H(t_{2})H(t_{1})]=θ(t_{2}−t_{1})H(t_{2})H(t_{1})+θ(t_{1}−t_{2})H(t_{1})H(t_{2}) using the Heaviside step function θ. With this notation, the propagator is commonly expressed as a timeordered evolution operator U(t)=[e^{−i∫}^{0}^{t}^{H(s)ds}], and Eq. (8) is written as
$\begin{array}{cc}\ue533\ue8a0\left[{e}^{i\ue89e{\int}_{0}^{t}\ue89eH\ue8a0\left(s\right)\ue89e\mathrm{ds}}\right]=\sum _{k=0}^{\infty}\ue89e{\left(i\right)}^{k}\ue89e{D}_{k},\text{}\ue89e{D}_{k}=\frac{1}{k!}\ue89e{\int}_{0}^{t}\ue89e\dots \ue89e{\int}_{0}^{t}\ue89e\ue533\ue8a0\left[H\ue8a0\left({t}_{k}\right)\ue89e\phantom{\rule{0.8em}{0.8ex}}\ue89e\dots \ue89e\phantom{\rule{0.8em}{0.8ex}}\ue89eH\ue8a0\left({t}_{1}\right)\right]\ue89e{d}^{k}\ue89et.& \left(9\right)\end{array}$
The Dyson series of Eq. (9) may be approximated on a computer by truncating the infinite expansion at some finite order k=K≥0, and evaluating H(t_{j}) at some finite number of M timesteps of size Δ=t/M. Thus, one has the approximation
$\begin{array}{cc}\ue533\ue8a0\left[{e}^{i\ue89e{\int}_{0}^{t}\ue89eH\ue8a0\left(s\right)\ue89e\mathrm{ds}}\right]\approx \sum _{k=0}^{K}\ue89e{\left(i\right)}^{k}\ue89e{D}_{k}\approx \sum _{k=0}^{\infty}\ue89e\frac{{\left(\mathrm{it}\right)}^{k}}{k!\ue89e{M}^{k}}\ue89e{\stackrel{~}{B}}_{k},\text{}\ue89e{\stackrel{~}{B}}_{k}=\sum _{{m}_{1},\dots ,{m}_{k}=0}^{M1}\ue89e\ue533\ue8a0\left[H\ue8a0\left({m}_{k}\ue89e\Delta \right)\ue89e\phantom{\rule{0.8em}{0.8ex}}\ue89e\dots \ue89e\phantom{\rule{0.8em}{0.8ex}}\ue89eH\ue8a0\left({m}_{1}\ue89e\Delta \right)\right],& \left(10\right)\end{array}$
which converges to U(t) in the limit K, M→∞ if H(t) is Riemann integrable. One can obtain some intuition on how the error of this approximation depends on choices of K, M, H(t). As ∥H(t)∥ is bounded,
$\uf605{D}_{k}\uf606\le \frac{{t}^{k}}{k!}\ue89e{\mathrm{max}}_{s}\ue89e{\uf605H\ue8a0\left(s\right)\uf606}^{k}.$
Using Taylor'"'"'s remainder theorem for the exponential function, the error of truncation for t max_{s }∥H(s)∥=(1) is at most ϵ=(1/K!). Thus the order of truncation is at most
$K=\ue52e\ue8a0\left(\frac{\mathrm{log}\ue8a0\left(1/\u03f5\right)}{\mathrm{log}\ue89e\phantom{\rule{0.3em}{0.3ex}}\ue89e\mathrm{log}\ue8a0\left(1/\u03f5\right)}\right).$
As {tilde over (B)}_{k }is a left Riemann sum of D_{k}, one expects that the error of discretization should scale like
$\u03f5=\ue52e\ue8a0\left(\frac{1}{M}\right).$
These two errors may then be added by a triangle inequality to obtain the overall error from approximating [e^{−i∫}^{0}^{t}^{H(s)ds}] with
$\sum _{k=0}^{K}\ue89e\frac{{\left(\mathrm{it}\right)}^{k}}{{M}^{k}}\ue89e{B}_{k}.$
The explicit timeordering in {tilde over (B)}_{k }may be removed with a slightly different approximation
$\begin{array}{cc}{\stackrel{~}{B}}_{k}=k!\ue89e{B}_{k}+{C}_{k},{B}_{k}=\sum _{0\le {m}_{1}<\dots <{m}_{k}<M}\ue89eH\ue8a0\left({m}_{k}\ue89e\Delta \right)\ue89e\phantom{\rule{0.8em}{0.8ex}}\ue89e\dots \ue89e\phantom{\rule{0.8em}{0.8ex}}\ue89eH\ue8a0\left({m}_{1}\ue89e\Delta \right),& \left(11\right)\end{array}$
where C_{k }captures terms where at least one pair of indices m_{j}=m_{k }collide for j≠k. A rigorous statement of this intuition is given by the following Lemma 1, with a detailed proof in Section XI.
Lemma 1 (Error from truncating and discretizing the Dyson series). Let H(s): →^{N×N }be differentiable on the domain [0, t], and let
$\u3008\uf605\stackrel{.}{H}\uf606\u3009\ue89e\text{:}=\frac{1}{t}\ue89e{\int}_{0}^{t}\ue89e\uf605\frac{\mathrm{dH}\ue8a0\left(s\right)}{\mathrm{ds}}\uf606\ue89e\mathrm{ds}.$
For any ϵ∈[0, 2^{1e}], an approximation to the time ordered operator exponential of −iH(s) can be constructed such that
$\uf605\ue533\ue8a0\left[{e}^{i\ue89e{\int}_{0}^{t}\ue89eH\ue8a0\left(s\right)\ue89e\mathrm{ds}}\right]\sum _{k=0}^{K}\ue89e{\left(\mathrm{it}/M\right)}^{k}\ue89e{B}_{k}\uf606\le \u03f5$
if one takes all the following are true.
$1.\ue89e\phantom{\rule{0.8em}{0.8ex}}\ue89e{\mathrm{max}}_{s}\ue89e\uf605H\ue8a0\left(s\right)\uf606\ue89et\le \mathrm{ln}\ue89e\phantom{\rule{0.3em}{0.3ex}}\ue89e2.$
$2.\ue89e\phantom{\rule{0.8em}{0.8ex}}\ue89eK=\lceil 1+\frac{2\ue89e\phantom{\rule{0.3em}{0.3ex}}\ue89e\mathrm{ln}\ue8a0\left(2/\u03f5\right)}{\mathrm{ln}\ue89e\phantom{\rule{0.3em}{0.3ex}}\ue89e\mathrm{ln}\ue8a0\left(2/\u03f5\right)+1}\rceil .\text{}\ue89e3.\ue89e\phantom{\rule{0.8em}{0.8ex}}\ue89eM\ge \mathrm{max}\ue89e\left\{\frac{16\ue89e\phantom{\rule{0.3em}{0.3ex}}\ue89e{t}^{2}}{\u03f5}\ue89e\left(\u3008\uf605\stackrel{.}{H}\uf606\u3009+{\mathrm{max}}_{s}\ue89e{\uf605H\ue8a0\left(s\right)\uf606}^{2}\right),{K}^{2}\right\}.$
On a quantum computer, it becomes possible to compute the B_{k }exactly and efficiently even it they sum over exponentially many points M. In contrast, computing these Riemann sums on classical computer would be prohibitive, even by approximate MonteCarlo sampling, which is exacerbated by the sign problem. However, this efficient quantum computation crucially assumes that information on terms of the Hamiltonian are made accessible in a certain coherent manner. In the case discussed here, this is information is accessed by querying the blackbox unitary oracle HAMT in Eq. (7).
B. Algorithm by Duplicating Control Registers
A quantum algorithm is presented that applies the ϵapproximation
$\stackrel{~}{U}=\sum _{k=0}^{\infty}\ue89e\frac{{\left(\mathrm{it}\right)}^{k}}{{M}^{k}}\ue89e{B}_{k}$
to the timeordered evolution operator [e^{−i∫}^{0}^{t}^{H(s)ds}], where the truncation order K and the number of discretization points M are given by Lemma 1. This version is the prelude to Theorem 1 and has worse space complexity. The algorithm proceeds in two steps: First, one synthesizes a unitary quantum circuit DYS_{K }that applies Ũ with success probability ½ on any input state to the s register. As B_{K }contains K a product of K Hamiltonian, this step requires (K) queries to HAMT. Second, since Ũ is ϵclose to unitary, one applies one round of oblivious amplitude amplification to boost this probability to 1−(ϵ). Among the contributions here are rigorous bounds on K and M, and the implementation of a step not discussed previously—the efficient preparation of a particular quantuni state that correctly selects a desired linear combination of timeordered products of Hamiltonians. This step is nonobvious as the state has (M!) different amplitudes, and in the worstcase would take (M!) gates to create by arbitrary state preparation techniques. The cost of this implementation is captured by the following theorem.
Theorem 2 (Hamiltonian simulation by a truncated Dyson series with duplicated registers). Let a timedependent Hamiltonian H(s) be characterized by spectral norm α≥max_{s }∥H(s)∥, and average rateofchange
$\u3008\uf605\stackrel{.}{H}\uf606\u3009=\frac{1}{t}\ue89e{\int}_{0}^{t}\ue89e\uf605\frac{\mathrm{dH}\ue8a0\left(s\right)}{\mathrm{ds}}\uf606\ue89e\mathrm{ds}.$
This Hamiltonian is defined for s ∈ [0,t], where αt≤½, and is encoded in an oracle HAMT acting on n_{s}+n_{α}+n_{d }qubits, as per Eq. (7), evaluated at
$M=\ue52e\left(\frac{{t}^{2}}{\u03f5}\ue89e\left(\frac{\u3008\uf605\stackrel{.}{H}\uf606\u3009}{\alpha}+\frac{{\mathrm{max}}_{s}\ue89e{\uf605H\ue8a0\left(s\right)\uf606}^{2}}{{\alpha}^{2}}\right)\right)$
points, where n_{d}=(log (M)). Then the timeordered evolution operator [exp (−i∫_{0}^{t}H(s)ds)] may be approximated to error ϵ with the following cost.
$1.\ue89e\phantom{\rule{0.8em}{0.8ex}}\ue89e\mathrm{Queries}\ue89e\phantom{\rule{0.8em}{0.8ex}}\ue89e\mathrm{to}\ue89e\phantom{\rule{0.8em}{0.8ex}}\ue89e\mathrm{HAM}\ue89e\text{}\ue89eT\ue89e\text{:}\ue89e\phantom{\rule{0.8em}{0.8ex}}\ue89e\ue52e\ue8a0\left(\frac{\mathrm{log}\ue8a0\left(t/\u03f5\right)}{\mathrm{log}\ue89e\phantom{\rule{0.3em}{0.3ex}}\ue89e\mathrm{log}\ue8a0\left(t/\u03f5\right)}\right).\text{}\ue89e2.\ue89e\phantom{\rule{0.8em}{0.8ex}}\ue89e\mathrm{Qubits}\ue89e\text{:}\ue89e\phantom{\rule{0.8em}{0.8ex}}\ue89e{n}_{s}+\ue52e\ue8a0\left({n}_{a}+\mathrm{log}\ue8a0\left(M\right)\ue89e\frac{\mathrm{log}\ue8a0\left(t/\u03f5\right)}{\mathrm{log}\ue89e\phantom{\rule{0.3em}{0.3ex}}\ue89e\mathrm{log}\ue8a0\left(t/\u03f5\right)}\right).\text{}\ue89e3.\ue89e\phantom{\rule{0.8em}{0.8ex}}\ue89e\mathrm{Primitive}\ue89e\phantom{\rule{0.8em}{0.8ex}}\ue89e\mathrm{gates}\ue89e\text{:}\ue89e\phantom{\rule{0.8em}{0.8ex}}\ue89e\ue52e\ue8a0\left(\left({n}_{a}+\mathrm{log}\ue8a0\left(M\right)\right)\ue89e\frac{\mathrm{log}\ue8a0\left(1/\u03f5\right)}{\mathrm{log}\ue89e\phantom{\rule{0.3em}{0.3ex}}\ue89e\mathrm{log}\ue8a0\left(1/\u03f5\right)}\right).$
Proof Let HAMT_{K }be a unitary that acts jointly on registers s, {right arrow over (α)}, {right arrow over (b)}, c, {right arrow over (d)}. This unitary is depicted in FIG. 2 and is defined to apply products of Hamiltonians
$\begin{array}{cc}\left(\u30080\ue89e{}_{\stackrel{>}{a}\ue89ee}\ue89e\otimes {I}_{s})\ue89e\mathrm{HAM}\ue89e\text{}\ue89e{{T}_{K}(0\u3009}_{\stackrel{>}{a}\ue89ee}\otimes {I}_{s}\right)\ue89e\text{:}=\hspace{1em}(\sum _{k=0}^{K}\ue89e\uf603k\u3009\ue89e\phantom{\rule{0.8em}{0.8ex}}\ue89e{\u3008k\uf604}_{\stackrel{>}{b}}\otimes \left(\sum _{\stackrel{>}{m}\in {\left[M\right]}^{k}}\ue89e\uf603\stackrel{>}{m}\u3009\ue89e\phantom{\rule{0.8em}{0.8ex}}\ue89e{\u3008\stackrel{>}{m}\uf604}_{{d}_{1}\ue89e\dots \ue89e\phantom{\rule{0.3em}{0.3ex}}\ue89e{d}_{k}}\otimes {I}_{{d}_{k+1}\ue89e\dots \ue89e\phantom{\rule{0.3em}{0.3ex}}\ue89e{d}_{K}}\otimes \left(\prod _{j=1}^{k}\ue89eH\ue8a0\left({m}_{j}\ue89e\Delta \right)\right)\right)+\dots \ue89e\phantom{\rule{0.8em}{0.8ex}})\otimes {\mathrm{SWAP}}_{c},& \left(12\right)\end{array}$
were SWAP_{c }swaps the two qubits of register c. Note that the action of HAMT_{K }is only defined for input states to register {right arrow over (b)} that are spanned by basis states of the unary encoding k_{{right arrow over (b)}}=0^{⊗k}1^{⊗Kk}, which determines the number of terms in the product. As seen in the figure, HAMT_{K }makes K queries to HAMT and copies the α, b, and d registers K times.
FIG. 2 shows a quantum circuit representations of the components implementing the truncated Dyson series algorithm TDS.
Now, consider a state s_{k}_{{right arrow over (d)}}, which is a uniform superposition of timeindex states under the dimensionk unit simplex
$\begin{array}{cc}{{{s}_{k}\u3009}_{\overrightarrow{d}}\ue89e\text{:}=\sqrt{\frac{k!\ue89e\left(Mk\right)!}{M!}}\ue89e{(\sum _{0\le {m}_{1}<{m}_{2}<\dots <{m}_{k}<M}\ue89e\overrightarrow{m}\u3009}_{{d}_{1}\ue89e\dots \ue89e\phantom{\rule{0.3em}{0.3ex}}\ue89e{d}_{k}})0\u3009}_{{d}_{k+1}\ue89e\dots \ue89e\phantom{\rule{0.3em}{0.3ex}}\ue89e{d}_{K}}.& \left(13\right)\end{array}$
This state is easy to prepare when k=1—there, it is simply a uniform superposition over M number states, and costs (log M) gates. Otherwise, naive methods based on rejection sampling have some success probability γk^{2 }that decreases exponentially with large k. Let PREP_{K }be one such unitary that prepares s_{k}_{{right arrow over (d)} }on measurement outcome 00_{c}.
PREP_{K}k_{{right arrow over (b)}}0_{c{right arrow over (d)}}:=k_{{right arrow over (b)}}(γk00_{c}s_{k}_{{right arrow over (d)}}+√{square root over (1−γk^{2})}01_{c }. . . ). (14)
For each order k, the Riemann sum B_{k }may be implemented by DYS_{K}:=(PREP_{K}^{†}⊗I_{αs}). HAMT_{K}·(PREP_{K}⊗I_{αs}), as depicted in FIG. 2. The unitary DYS_{K }encodes precisely terms B_{k }of the Dyson series as follows
$\begin{array}{cc}\left(\u30080\ue89e{}_{\overrightarrow{a}\ue89ec\ue89e\phantom{\rule{0.3em}{0.3ex}}\ue89e\overrightarrow{d}}\ue89e\otimes {I}_{\overrightarrow{b}\ue89es})\ue89e{{\mathrm{DYS}}_{K}(0\u3009}_{\overrightarrow{a}\ue89ec\ue89e\phantom{\rule{0.3em}{0.3ex}}\ue89e\overrightarrow{d}}\otimes {I}_{\overrightarrow{b}\ue89es}\right)=\sum _{k=0}^{K}\ue89e\uf603k\u3009\ue89e\phantom{\rule{0.8em}{0.8ex}}\ue89e{\u3008k\uf604}_{\overrightarrow{b}}\otimes \frac{{\uf603{\gamma}_{k}\uf604}^{2}\ue89ek!\ue89e\left(Mk\right)!}{M!}\ue89e{B}_{k}.& \left(15\right)\end{array}$
Now, a linear combination of Dyson series terms is implemented by preparing a state with the appropriate amplitudes in the basis k_{{right arrow over (b)}}. The required state preparation operators are
$\begin{array}{cc}{{{{\phantom{\rule{4.4em}{4.4ex}}\ue89e{\mathrm{COEF}}_{K}0\u3009}_{\overrightarrow{b}}\ue89e\text{:}=\frac{1}{\sqrt{\beta}}\ue89e\sum _{k=0}^{K}\ue89e\sqrt{\frac{M!\ue89e{\left(\mathrm{it}\right)}^{k}}{{M}^{k}\ue89e{\uf603{\gamma}_{k}\uf604}^{2}\ue89ek!\ue89e\left(Mk\right)!}}k\u3009}_{\stackrel{>}{b}},\text{}\ue89e\phantom{\rule{1.1em}{1.1ex}}\ue89e\beta =\sum _{k=0}^{K}\ue89e\frac{M!\ue89e{t}^{k}}{{M}^{k}\ue89e{\uf603{\gamma}_{k}\uf604}^{2}\ue89ek!\ue89e\left(Mk\right)!},\text{}\ue89e{\mathrm{COEF}}_{K}^{\prime}0\u3009}_{\stackrel{>}{b}}\ue89e\text{:}=\frac{1}{\sqrt{\beta}}\ue89e\sum _{k=0}^{K}\ue89e\sqrt{\frac{M!\ue89e{t}^{k}}{{M}^{k}\ue89e{\uf603{\gamma}_{k}\uf604}^{2}\ue89ek!\ue89e\left(Mk\right)!}}k\u3009}_{\stackrel{>}{b}},& \left(16\right)\end{array}$
and may be implemented using (K) primitive gates. Up to a proportionality factor β, one can obtain the desired linear combination for simulating timeevohnion.
$\begin{array}{cc}{\phantom{\rule{4.4em}{4.4ex}}\ue89e{\mathrm{TDS}}_{\beta}:=\left({\mathrm{COEF}}_{K}^{\mathrm{\prime \u2020}}\otimes {I}_{\overrightarrow{a}\ue89e\phantom{\rule{0.3em}{0.3ex}}\ue89ec\ue89e\phantom{\rule{0.3em}{0.3ex}}\ue89e\overrightarrow{d}\ue89e\phantom{\rule{0.3em}{0.3ex}}\ue89es}\right)\xb7{\mathrm{DYS}}_{K}\xb7\left({\mathrm{COEF}}_{K}\otimes {I}_{\overrightarrow{a}\ue89e\phantom{\rule{0.3em}{0.3ex}}\ue89ec\ue89e\phantom{\rule{0.3em}{0.3ex}}\ue89e\overrightarrow{d}\ue89e\phantom{\rule{0.3em}{0.3ex}}\ue89es}\right)\ue89e\text{}\ue89e\left({\u30080\uf604}_{\overrightarrow{a}\ue89e\phantom{\rule{0.3em}{0.3ex}}\ue89e\overrightarrow{b}\ue89e\phantom{\rule{0.3em}{0.3ex}}\ue89ec\ue89e\phantom{\rule{0.3em}{0.3ex}}\ue89e\overrightarrow{d}}\otimes {I}_{s}\right)\ue89e{\mathrm{TDS}}_{\beta}(\uf604\ue89e0\u3009}_{\overrightarrow{a}\ue89e\phantom{\rule{0.3em}{0.3ex}}\ue89e\overrightarrow{b}\ue89e\phantom{\rule{0.3em}{0.3ex}}\ue89ec\ue89e\phantom{\rule{0.3em}{0.3ex}}\ue89e\overrightarrow{d}}\otimes {I}_{s})=\frac{{\sum}_{k=0}^{K}\ue89e\phantom{\rule{0.3em}{0.3ex}}\ue89e{\left(\mathrm{it}\right)}^{k}\ue89e{B}_{k}}{{M}^{k}\ue89e\beta}\approx \frac{\ue533\ue89e\phantom{\rule{0.3em}{0.3ex}}\ue89e{e}^{i\ue89e{\int}_{0}^{t}\ue89eH\ue8a0\left(s\right)\ue89e\phantom{\rule{0.2em}{0.2ex}}\ue89e\mathrm{ds}}}{\beta}.& \left(17\right)\end{array}$
If one chooses t to be sufficiently small such that β=2, one can then obtain a single timestep of the truncated Dyson series algorithm TDS in FIG. 2. All that remains is to find an implementation of PREP_{K }that prepares s_{k}_{{right arrow over (d)} }with an amplitude that γk that is sufficiently large so that t=⊖(1).
The state s_{k}_{{right arrow over (b)}c{right arrow over (d)} }can be prepared in a number of ways. The most straightforward approach creates a uniform superposition of states over the dimensionk hypercube using n_{d}×k Hadamard gates HAD, then uses k reversible adders to flag states {right arrow over (m)}_{d}_{1 }_{. . . d}_{k }with the correct ordering. This circuit produces s_{k}_{{right arrow over (d)} }with amplitude
${\gamma}_{k}=\sqrt{\frac{M!}{{M}^{k}\ue89ek!\ue89e\left(Mk\right)!}}.\phantom{\rule{0.6em}{0.6ex}}\ue89e{\mathrm{PREP}}_{K}$
is then obtained by controlling on input state k_{{right arrow over (b)}}. Thus
$\begin{array}{cc}\beta =\sum _{k=0}^{K}\ue89e\phantom{\rule{0.3em}{0.3ex}}\ue89e{t}^{k}\le \sum _{k=0}^{\infty}\ue89e\phantom{\rule{0.3em}{0.3ex}}\ue89e{t}^{k}=\frac{1}{1t}.& \left(18\right)\end{array}$
Thus by choosing t=⊖(1)≈½, s=2 and a single round of oblivious amplitude amplification suffices. Notably, even though the success probability of naive state preparation γk^{2 }decays rapidly, this only amounts to a constant factor slowdown compared to more sophisticated techniques that effectively prepare s_{k}_{{right arrow over (d)} }with success probability≈1. For example, rather than rejection sampling, one may perform a reversible sort on on uniform superposition of states
$\begin{array}{c}\frac{1}{\sqrt{{M}^{k}}}\ue89e{\sum}_{\overrightarrow{m}}\ue89e{\uf603\overrightarrow{m}\u3009}_{{d}_{1}\ue89e\dots \ue89e\phantom{\rule{0.3em}{0.3ex}}\ue89e{d}_{k}}\ue89e{\uf6030\u3009}_{\mathrm{garbage}}\to \hspace{1em}\hspace{1em}\ue89e\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\ue89e\hspace{1em}\frac{1}{\sqrt{{M}^{k}}}\ue89e{\sum}_{\overrightarrow{m}}\ue89e{\uf603\ue533\ue8a0\left[{m}_{{d}_{1}}\ue89e\phantom{\rule{0.8em}{0.8ex}}\ue89e\dots \ue89e\phantom{\rule{0.8em}{0.8ex}}\ue89e{m}_{{d}_{k}}\right]\u3009}_{{d}_{1}\ue89e{\mathrm{\dots d}}_{k}}\uf604\\ {\overrightarrow{m}\u3009}_{\mathrm{garbage}},\end{array}$
such as with the quanthm bitonic sorting network. This effectively increases γ_{k}^{2 }by a factor of k!, and uses significantly more ancilla qubits, but ultimately allows us to implement time steps t≈ln2≈0.693 larger by a constant factor. □
C. Compression of Control Registers
The dominant contribution to the space overhead of the truncated Dyson series algorithm in Section III B is the Kfold duplication of registers {right arrow over (α)}, {right arrow over (d)} in the oracle DYS_{K}. A general technique is now presented that when applied to DYS_{K }avoids this duplication and completes the proof of Theorem 1. Suppose one has a sequence of K unitary oracles U_{1}, U_{2}, . . . , U_{K}, defined as
(0_{α}⊗I_{s})U_{k}(0_{α}⊗I_{s})=H_{k}, (19)
where H_{j }is some arbitrary matrix with bounded spectral norm ∥H_{j}∥≤1. In the general problem, one can construct a quantum circuit V that when controlled on index k in register b, applies the sequence U_{k }. . . U_{2}U_{1}, that is
$\begin{array}{cc}{\left({\u30080\uf604}_{\mathrm{ac}}\otimes {I}_{s}\right)\ue89eV(\uf604\ue89e0\u3009}_{\mathrm{ac}}\otimes {I}_{s})=\uf604\ue89e0\u3009\ue89e{\u30080\uf604}_{b}\otimes {I}_{s}=\sum _{k=1}^{K}\ue89e\phantom{\rule{0.3em}{0.3ex}}\ue89e\uf603k\u3009\ue89e{\u3008k\uf604}_{b}\otimes \left(\prod _{j=1}^{k}\ue89e\phantom{\rule{0.3em}{0.3ex}}\ue89e{H}_{j}\right).& \left(20\right)\end{array}$
Though binary control logic for this sequence is trivial when H_{k }is unitary, the complication here is that H_{k }is in general nonunitary and so the probability of successfully measuring 0_{α }is less than one. Any other measurement outcome corresponds failure as it applies on register α an operator that is not H_{k}. This complication is overcome by introducing two more registers b, c of size (log(K)) gubits that count the number of successful measurements.
Lemma 2 (Compression gadget). Let {U_{k},k∈[K]} be a set of K unitaries that encode matrices H_{k }as defined in Eq. (19). Then there exists a quantum circuit V satisfying Eq. (20) such that the number of qubits n_{b}=(log (K)), n_{c}=(log (K)). The cost of V is one query to each of controlledcontrolledU_{k}, and (K (n_{α}+log(K))) additional primitive quantum gates.
Proof. Let b, c be counter registers. Using notation where commas in the subscript indicates respective assignments, these represent an n_{b,c}bit integer l_{b,c}=Σ_{r=0}^{n}^{b,c}^{−1}2^{r}q_{r }in the number state l_{b,c}_{b,c}:=q_{0}q_{1 }. . . q_{n}_{b,c }. . . 1_{b,c}, where q_{r }∈{0, 1}. The size of these integers are determined by n_{b}=n_{c}+1=┌log_{2}(K+1)┐+1. The unitaries U_{j }will be applied conditional on the trailing bit q_{0}=0 in the c register, and the leading bit q_{n}_{b}_{−1}=0 in the b register, that is
CC−U_{k}:=I^{⊗n}^{b}^{n}^{c}^{−2}⊗(00_{b}_{nb}_{−1}⊗00_{c}_{0}⊗U_{k}+ . . . ). (21)
Consider the circuit in FIG. 3. There, one can apply CCU_{k}, then increment k by one, increment l_{c }by one conditional on the α register not being in the 0_{α }state, and decrement l_{b }by one conditional on the α register being in the 0_{α }state. This is accomplished by multiplycontrolled modular addition
$\begin{array}{cc}{\mathrm{ADD}}_{\mathrm{ca}}={\mathrm{ADD}}_{b}^{\u2020}\otimes {I}_{c}\otimes \uf6030\u3009\ue89e{\u30080\uf604}_{a}+{I}_{b}\otimes {\mathrm{ADD}}_{c}\otimes \sum _{l=1}^{{2}^{{n}_{a}}1}\ue89e\phantom{\rule{0.3em}{0.3ex}}\ue89e\uf603l\u3009\ue89e{\u3008l\uf604}_{a},\text{}\ue89e{\mathrm{ADD}}_{b,c}=\sum _{l=0}^{{2}^{{n}_{b,c}}1}\ue89e\phantom{\rule{0.3em}{0.3ex}}\ue89e\uf603l+1\ue89e\phantom{\rule{0.8em}{0.8ex}}\ue89e\mathrm{mod}\ue89e\phantom{\rule{0.8em}{0.8ex}}\ue89e{2}^{{n}_{b,c}}\u3009\ue89e{\u3008l\uf604}_{b,c}.& \left(22\right)\end{array}$
As one adds integers of size (K), each application of modular addition costs (log(K)) primitive gates and requires (log(K)) qubits. Implementing the multiple controls costs (n_{α}) primitive gates and up to n_{α }extra qubits. FIG. 3 is a quantum circuit representations of the gadget V for applying a sequence of probabilistic operators H_{k }. . . H_{2}H_{1}, encoded in (0_{α}⊗I_{s})U_{k}(0_{α}⊗I_{s})=H_{k}, controlled on number state k_{b}, k∈{0,1, . . . , K}.
Restricted to input states 0_{α}l_{b}0_{c}, where l ∈{2^{n}^{b}−1,0,1,2,3, . . . , K−1}, this implements V. For example, consider the evolution of an input state 0_{α}2_{b}0_{c}φ_{s }for K=3.
$\begin{array}{cc}{{{{\uf6030\u3009}_{a}\uf604\ue89e1\u3009}_{b}\ue89e{\uf6030\u3009}_{c}\uf604\ue89e\psi \u3009}_{s}\ue89e\underset{\mathrm{CC}{U}_{1}}{\to}\ue89e{\uf6030\u3009}_{a}\ue89e{\uf6031\u3009}_{b}\ue85c0\u3009}_{c}\ue89e{H}_{1}\ue89e{\uf603\psi \u3009}_{s}+{\uf603{0}^{\perp \mathrm{.1}}\u3009}_{a}\ue89e{\uf6032\u3009}_{b}\ue89e{\uf6030\u3009}_{c}\ue89e\phantom{\rule{0.8em}{0.8ex}}\ue89e\dots \ue89e\phantom{\rule{0.8em}{0.8ex}}\ue89e\underset{{\mathrm{ADD}}_{\mathrm{ca}}}{\to}\ue89e{\hspace{1em}{\uf6030\u3009}_{a}\ue89e{\uf6030\u3009}_{b}0\u3009}_{c}\ue89e{H}_{1}\ue89e{\uf603\psi \u3009}_{s}+{\uf603{0}^{\perp \mathrm{.1}}\u3009}_{a}\ue89e{\uf6032\u3009}_{b}\ue89e{\uf6031\u3009}_{c}\ue89e\phantom{\rule{0.8em}{0.8ex}}\ue89e\dots \ue89e\phantom{\rule{0.6em}{0.6ex}}\ue89e\underset{\mathrm{CC}{U}_{2}}{\to}\ue89e\phantom{\rule{0.3em}{0.3ex}}{\uf6030\u3009}_{a}\ue89e{\uf6030\u3009}_{b}\ue89e{\uf6030\u3009}_{c}\ue89e{H}_{2}\ue89e{H}_{1}\ue89e{\uf603\psi \u3009}_{s}+{\uf603{0}^{\perp \mathrm{.2}}\u3009}_{a}\ue89e{\uf6031\u3009}_{b}\ue89e{\uf6030\u3009}_{c}\ue89e\phantom{\rule{0.8em}{0.8ex}}\ue89e\dots +{\uf603{0}^{\perp \mathrm{.1}}\u3009}_{a}\ue89e{\uf6032\u3009}_{b}\ue89e{\uf6031\u3009}_{c}\ue89e\phantom{\rule{0.8em}{0.8ex}}\ue89e\dots \ue89e\phantom{\rule{0.8em}{0.8ex}}\ue89e\underset{{\mathrm{ADD}}_{\mathrm{ca}}}{\to}\ue89e{\uf6030\u3009}_{a}\ue89e{\uf6031\u3009}_{b}\ue89e{\uf6030\u3009}_{c}\ue89e{H}_{2}\ue89e{H}_{1}\ue89e{\uf603\psi \u3009}_{s}+{\uf603{0}^{\perp \mathrm{.2}}\u3009}_{a}\ue89e{\uf6031\u3009}_{b}\ue89e{\uf6031\u3009}_{c}\ue89e\phantom{\rule{0.8em}{0.8ex}}\ue89e\dots +{\uf603{0}^{\perp \mathrm{.1}}\u3009}_{a}\ue89e\uf603{2}_{b}\ue89e{\uf6032\u3009}_{c}\ue89e\phantom{\rule{0.8em}{0.8ex}}\ue89e\dots \ue89e\phantom{\rule{0.8em}{0.8ex}}\ue89e\underset{\mathrm{CC}{U}_{3}}{\to}\ue89e{\uf6030\u3009}_{a}\ue89e{\uf6031\u3009}_{b}\ue89e{\uf6030\u3009}_{c}\ue89e{H}_{2}\ue89e{H}_{1}\ue89e{\uf603\psi \u3009}_{s}+{\uf603{0}^{\perp \mathrm{.2}}\u3009}_{a}\ue89e{\uf6031\u3009}_{b}\ue89e{\uf6031\u3009}_{c}\ue89e\phantom{\rule{0.8em}{0.8ex}}\ue89e\dots +{\uf603{0}^{\perp \mathrm{.1}}\u3009}_{a}\ue89e{\uf6032\u3009}_{b}\ue89e{\uf6032\u3009}_{c}\ue89e\phantom{\rule{0.8em}{0.8ex}}\ue89e\dots \ue89e\phantom{\rule{0.8em}{0.8ex}}\ue89e\underset{{\mathrm{ADD}}_{\mathrm{ca}}}{\to}\ue89e{\uf6030\u3009}_{a}\ue89e{\uf6032\u3009}_{b}\ue89e{\uf6030\u3009}_{c}\ue89e{H}_{2}\ue89e{H}_{1}\ue89e{\uf603\psi \u3009}_{s}+{\uf603{0}^{\perp \mathrm{.2}}\u3009}_{a}\ue89e{\uf6031\u3009}_{b}\ue89e{\uf6032\u3009}_{c}\ue89e\phantom{\rule{0.8em}{0.8ex}}\ue89e\dots +{\uf603{0}^{\perp \mathrm{.1}}\u3009}_{a}\ue89e{\uf6032\u3009}_{b}\ue89e{\uf6033\u3009}_{c}\ue89e\phantom{\rule{0.8em}{0.8ex}}\ue89e\dots & \left(23\right)\end{array}$
In the above, negative integers correspond to the trailing bit q_{n}_{b}_{−1}=1 so U_{k }in Eq. (21) is not applied. Note that Eq. (20) applies H_{k }. . . H_{1 }controlled on k_{b}, so one can simply relabel k=l_{b}+1 mod 2^{n}^{b}.
Proof of Theorem 1. Thus DYS_{K }may be implemented with reduced ancilla overhead through Lemma 2 provided that one finds a sequence {U_{k}} such that H_{k }. . . H_{2}H_{1}∝B_{k}—in other words,
$\begin{array}{cc}{\left({\u30080\uf604}_{\mathrm{ac},\mathrm{others}}\otimes {I}_{\mathrm{bs}}\right)\ue89e{\mathrm{DYS}}_{K}(\uf604\ue89e0\u3009}_{\mathrm{ac},\mathrm{others}}\otimes {I}_{\mathrm{bs}})=\sum _{k=0}^{K}\ue89e\phantom{\rule{0.3em}{0.3ex}}\ue89e\uf603k\u3009\ue89e{\u3008k\uf604}_{b}\otimes {\gamma}_{k}\ue89e{B}_{k},& \left(24\right)\end{array}$
where ‘others’ represent registers with size independent of K, and γ_{k }is a scaling factor depends on the choice of U_{k}. This sequence is obtained by combining three matrices. First, a unitary matrix U that prepares a uniform superposition
${{U0\u3009}_{d}={\sum}_{m=0}^{M1}\ue89e\phantom{\rule{0.3em}{0.3ex}}\ue89e\frac{1}{\sqrt{M}}\uf604\ue89em\u3009}_{d}.$
Second, the blockdiagonal matrix D=Σ_{m=0}^{M−1}mm⊗H(Δm) implemented by HAMT. Third, a strictly uppertriangular matrix G ∈^{M×M }with elements
$\begin{array}{cc}{G}_{\mathrm{ij}}=\{\begin{array}{cc}\frac{1}{M},& i<j,\\ 0& \mathrm{otherwise},\end{array}\ue89eG=\frac{1}{M}\ue89e\sum _{i=0}^{M1}\ue89e\phantom{\rule{0.3em}{0.3ex}}\ue89e\sum _{j=i+1}^{M1}\ue89e\phantom{\rule{0.3em}{0.3ex}}\uf604\ue89ei\u3009\ue89e\u3008j\uf604.& \left(25\right)\end{array}$
The nonunitary triangular operator G is implemented by using an integer comparator COMP acting on registers d, e, f consisting of n_{d}=n_{e }and n_{f}=1 qubits. For any input number state index j_{d}, let one compare j with a uniform superposition state
${{\sum}_{i=0}^{M1}\ue89e\phantom{\rule{0.3em}{0.3ex}}\ue89e\frac{1}{\sqrt{M}}\uf604\ue89ei\u3009}_{e}.$
Conditional on i≥j, perform a NOT gate on register f. One can then swap registers d,e, and unprepare the uniform superposition. On input j_{d}0_{e}0_{f}, this implements the sequence.
$\begin{array}{cc}{{{{{\uf604j\u3009}_{d}\uf604\ue89e0\u3009}_{e}\uf604\ue89e0\u3009}_{f}\ue89e\underset{U\ue89e\phantom{\rule{0.6em}{0.6ex}}\ue89e\mathrm{on}\ue89e\phantom{\rule{0.8em}{0.8ex}}\ue89ee}{\to}\ue89e\frac{{\uf604j\u3009}_{d}}{\sqrt{M}}\ue89e\sum _{i=0}^{M1}\ue89e{\uf603i\u3009}_{e}\uf604\ue89e0\u3009}_{f}\ue89e\underset{\mathrm{COMP}}{\to}\ue89e\frac{{j\u3009}_{d}}{\sqrt{M}}\ue89e\sum _{i=0}^{M1}i\u3009}_{e}i\ge {\hspace{1em}j\u3009}_{f}\ue89e\underset{{\mathrm{SWAP}}_{d\in}}{\to}& \left(26\right)\\ {{{{{\frac{{j\u3009}_{e}}{\sqrt{M}}\ue89e\sum _{i=0}^{M1}i\u3009}_{d}i\ge j\u3009}_{f}\ue89e\underset{{{U}^{\u2020}}_{\mathrm{on}\ue89e\phantom{\rule{0.6em}{0.6ex}}\ue89ee}}{\to}\ue89e\frac{1}{M}\ue89e\sum _{i=0}^{M1}i\u3009}_{d}0\u3009}_{e}i\ge j\u3009}_{f}+\dots \ue89e\phantom{\rule{0.8em}{0.8ex}},& \phantom{\rule{0.3em}{0.3ex}}\end{array}$
where i≥j_{f}=1_{f }if i≥j and is 0_{f }if i<j. This defines the following circuit that encodes G.
$\begin{array}{cc}\mathrm{LT}=\left({I}_{f}\otimes U\otimes {I}_{d}\right)\xb7\mathrm{COMP}\xb7\left({I}_{f}\otimes {\mathrm{SWAP}}_{\mathrm{de}}\right)\xb7\left({I}_{f}\otimes {U}^{\u2020}\otimes {I}_{d}\right),{\hspace{1em}\left({\u30080\uf604}_{\mathrm{ef}}\otimes {I}_{d}\right)\ue89e\mathrm{LT}(\uf604\ue89e0\u3009}_{\mathrm{ef}}\otimes {I}_{d})=\frac{1}{M}\ue89e\sum _{i=0}^{M1}\ue89e\phantom{\rule{0.3em}{0.3ex}}\ue89e\sum _{j=i+1}^{M1}\ue89e\phantom{\rule{0.3em}{0.3ex}}\ue89e\uf603i\u3009\ue89e{\u3008j\uf604}_{d}.& \left(27\right)\end{array}$
The gate complexity of G is dominated by the comparator COMP, which costs (n_{d}) primitive quantum gates.
FIG. 5 is a schematic block diagram showing quantmn circuit representations of the unitary DYS_{K }that encodes timeordered products of Hamiltonians, using fewer ancilla than the construction in FIG. 2 by applying the compression gadget in FIG. 4.
One may then verify that the terms B_{k }are generated by the following sequence
$\begin{array}{cc}{\u30080\ue89e{\uf603}_{d}\ue89e{U}^{\u2020}\xb7D\xb7U\uf604\ue89e0\u3009}_{d}=\frac{{B}_{1}}{M}=\frac{1}{M}\ue89e\sum _{{m}_{1}=0}^{M1}\ue89e\phantom{\rule{0.3em}{0.3ex}}\ue89eH\ue89e\left(\Delta \ue89e\phantom{\rule{0.3em}{0.3ex}}\ue89em\right),& \left(28\right)\\ {\u30080\ue89e{\uf603}_{d}\ue89e{U}^{\u2020}\xb7\left(D\xb7G\right)\xb7D\xb7U\uf604\ue89e0\u3009}_{d}=\frac{{B}_{2}}{{M}^{2}}=\frac{1}{{M}^{2}}\ue89e\sum _{0\le {m}_{1}<{m}_{2}<M}^{\phantom{\rule{0.3em}{0.3ex}}}\ue89e\phantom{\rule{0.3em}{0.3ex}}\ue89eH\ue8a0\left(\Delta \ue89e\phantom{\rule{0.3em}{0.3ex}}\ue89e{m}_{2}\right)\ue89eH\ue89e\left(\Delta \ue89e\phantom{\rule{0.3em}{0.3ex}}\ue89e{m}_{1}\right),& \phantom{\rule{0.3em}{0.3ex}}\\ \vdots & \phantom{\rule{0.3em}{0.3ex}}\\ {\u30080\ue89e{\uf603}_{d}\ue89e{U}^{\u2020}\xb7{\left(D\xb7G\right)}^{k1}\xb7D\xb7U\uf604\ue89e0\u3009}_{d}=\frac{{B}_{k}}{{M}^{k}}=& \phantom{\rule{0.3em}{0.3ex}}\\ \frac{1}{{M}^{k}}\ue89e\sum _{0\le {m}_{1}<{m}_{2}<\dots \ue89e\phantom{\rule{0.3em}{0.3ex}}\ue89e{m}_{j}<M}^{\phantom{\rule{0.3em}{0.3ex}}}\ue89e\phantom{\rule{0.3em}{0.3ex}}\ue89eH\ue8a0\left(\Delta \ue89e\phantom{\rule{0.3em}{0.3ex}}\ue89e{m}_{j}\right)\ue89e\dots \ue89e\phantom{\rule{0.3em}{0.3ex}}\ue89eH\ue8a0\left(\Delta \ue89e\phantom{\rule{0.3em}{0.3ex}}\ue89e{m}_{2}\right)\ue89e\dots \ue89e\phantom{\rule{0.3em}{0.3ex}}\ue89eH\ue8a0\left(\Delta \ue89e\phantom{\rule{0.3em}{0.3ex}}\ue89e{m}_{1}\right).& \phantom{\rule{0.3em}{0.3ex}}\end{array}$
Thus one can make the choice
$\begin{array}{cc}{U}_{k}:=\{\begin{array}{cc}\left({U}^{\u2020}\otimes {I}_{\mathrm{aefs}}\right)\xb7\left(\mathrm{HAM}\ue89e\text{}\ue89eT\otimes {I}_{\mathrm{ef}}\right)\xb7\left(U\otimes {I}_{\mathrm{aefs}}\right),& k=1,\\ \left({U}^{\u2020}\otimes {I}_{\mathrm{aefs}}\right)\xb7\left(\mathrm{HAM}\ue89e\text{}\ue89eT\otimes {I}_{\mathrm{ef}}\right)\xb7\left(\mathrm{LT}\otimes {I}_{\mathrm{as}}\right)\xb7\left(U\otimes {I}_{\mathrm{aefs}}\right),& k>1.\end{array}& \left(29\right)\end{array}$
Combined with Lemma 2, this leads to the circuit of FIG. 5, which implements DYS_{K }in Eq. (24) by identifying ‘others’ with the e and f registers, and scaling factor
${\gamma}_{k}=\frac{1}{{M}^{k}}.$
$\begin{array}{cc}\left({\u30080\uf604}_{\mathrm{acef}}\otimes {I}_{\mathrm{bs}}\right)\ue89e{\mathrm{DYS}}_{K}\ue8a0\left({\uf6030\u3009}_{\mathrm{acef}}\otimes {I}_{\mathrm{bs}}\right)=\sum _{k=0}^{K}\ue89e\phantom{\rule{0.3em}{0.3ex}}\ue89e\uf603k\u3009\ue89e{\u3008k\uf604}_{b}\otimes \frac{{B}_{k}}{{M}^{k}},& \left(30\right)\end{array}$
One can then select the desired linear combination of different orders in the Dyson series with the state
$\begin{array}{cc}{\mathrm{COEF}}_{K}\ue89e{\uf6030\u3009}_{b}=\frac{1}{\sqrt{\beta}}\ue89e\sum _{k=0}^{K}\ue89e\sqrt{{\left(\mathrm{it}\right)}^{k}}\ue89e{\uf603k\u3009}_{b},& \left(31\right)\\ {\mathrm{COEF}}_{K}^{\prime}\ue89e{\uf6030\u3009}_{b}=\frac{1}{\sqrt{\beta}}\ue89e\sum _{k=0}^{K}\ue89e\sqrt{{t}^{k}}\ue89e{\uf603k\u3009}_{b},\beta =\sum _{j=0}^{K}\ue89e{t}^{k}\le \sum _{k=0}^{\infty}\ue89e{t}^{k}=\frac{1}{1t}.& \left(32\right)\end{array}$
where one assumes that t<1 for convergence. Thus by choosing t=⊖(1)≈½, the success probability of approximating e^{−1∫}^{0}^{t}^{H(s)ds }may be boosted from β^{2}=4to 1−(ϵ) using a single round of oblivious amplitude amplification. □
IV. Accelerated Interaction Picture Simulation
Timeindependent Hamiltonians H become timedependent H_{I}(t) in the interaction picture. This requires the use of timedependent Hamiltonian simulation algorithm, which scale with parameters of H_{I}(t) that differ from those for the timeindependent case. For certain broad classes of Hamiltonian identified in Section IV A, these different dependencies allow us to improve the gate complexity of approximating the timeevolution operator e^{−iHt }by instead performing simulation in the interaction picture using the truncated Dyson series algorithm Theorem 1.
The interaction picture can be viewed as an intermediate between the Schrödinger and Heisenberg pictures wherein some of the dynamics is absorbed into the state and the remainder is absorbed into the dynamics of the operators. If the Hamiltonian in the Schrödinger picture is H=A+B and φ(t)=e^{−iHt}φ(0) then the Hamiltonian in the interaction picture is H_{I}(t)=e^{iAt}Be^{−iAt }and i∂_{t}φ_{I}(t)=H_{I}(t)φ_{I}(t) with φ_{I}(t)=e^{iAt}φ(t) for all t. These relations can easily be seen by substituting into the Schrödinger equation:
$\begin{array}{cc}i\ue89e{\partial}_{t}\ue89e\uf603{\psi}_{I}\ue8a0\left(t\right)\u3009=i\ue89e{\partial}_{t}\ue89e\left({e}^{\mathrm{iAt}}\ue89e\uf603\psi \ue8a0\left(t\right)\u3009\right)={e}^{\mathrm{iAt}}\ue8a0\left(A+H\right)\ue89e\uf603\psi \ue8a0\left(t\right)\u3009={e}^{\mathrm{iAt}}\ue89e{\mathrm{Be}}^{\mathrm{iAt}}\ue89e{e}^{\mathrm{iAt}}\ue89e\uf603\psi \ue8a0\left(t\right)\u3009={H}_{I}\ue8a0\left(t\right)\ue89e\uf603{\psi}_{I}\ue8a0\left(t\right)\u3009.& \left(33\right)\end{array}$
Note that if one started with timedependent B(t), that is H(t)=A+B(t), the interaction picture Hamiltonian is H_{I}(t)=e^{iAt}B(t)e^{−iAt}. The following results generalize easily to this situation, and so one can consider timeindependent B for simplicity.
The advantage of this representation is a Hanmiltonian with a smaller norm ∥H(t)∥=∥B∥≤∥A∥+∥B∥, but at the price of introducing timedependence. In general, one cannot write a closed form expression for the timeevolution operator. The following notation is commonly used to express the time evolution operator [e^{−i∫}^{0}^{t}^{H(s)ds}]=lim_{r→∞}II_{j=1}^{r}e^{−iH(jt/r)t/r }where this product is implicitly defined to be time ordered. Given an initial state φ(0), the state after evolution for t>0 may thus be written as
$\begin{array}{cc}\uf603\psi \ue8a0\left(t\right)\u3009={e}^{\mathrm{iAt}}\ue89e\uf603{\psi}_{I}\ue8a0\left(t\right)\u3009={e}^{\mathrm{iAt}}\ue89e\ue533\ue8a0\left[{e}^{i\ue89e{\int}_{0}^{t}\ue89e{H}_{I}\ue8a0\left(s\right)\ue89e\mathrm{ds}\ue89e\phantom{\rule{0.2em}{0.2ex}}}\right]\ue89e\uf603{\psi}_{I}\ue8a0\left(0\right)\u3009={e}^{\mathrm{iAt}}\ue89e\ue533\ue8a0\left[{e}^{i\ue89e{\int}_{0}^{t}\ue89e{H}_{I}\ue8a0\left(s\right)\ue89e\mathrm{ds}\ue89e\phantom{\rule{0.2em}{0.2ex}}}\right]\ue89e\uf603\psi \ue8a0\left(0\right)\u3009={\left({e}^{\mathrm{iAt}}\ue89e\ue533\ue8a0\left[{e}^{i\ue89e{\int}_{0}^{t}\ue89e{H}_{I}\ue8a0\left(s\right)\ue89e\mathrm{ds}\ue89e\phantom{\rule{0.2em}{0.2ex}}}\right]\right)}^{L}\ue89e\uf603\psi \ue8a0\left(0\right)\u3009,& \left(34\right)\end{array}$
and evolution by the full duration is decomposed into evolution by L shorter segments of duration .
Using the simulation algorithm Theorem 1 to simulate each segment in Eq (34) leads to the following result
Lemma 3 (Query complexity of Hamiltonian simulation in the interaction picture). For any Hamiltonian H=A+B, let HAMT be a unitary oracle that encodes H_{I}(t)=e^{iAt}Be^{−iAt }at
$M=\ue52e\ue8a0\left(\frac{t}{\u03f5}\ue89e\frac{\uf605A\uf606}{{\alpha}_{B}}\right)$
uniformly spaced values of t ∈[0,=(α_{B}^{−1})].
$\begin{array}{cc}\left({\u30080\uf604}_{a}\otimes {I}_{s}\right)\ue89e\mathrm{HAM}\ue89e\text{}\ue89eT\ue8a0\left({\uf6030\u3009}_{a}\otimes {I}_{s}\right)=\sum _{m=0}^{M1}\ue89e\phantom{\rule{0.3em}{0.3ex}}\ue89e\uf603m\u3009\ue89e{\u3008m\uf604}_{b}\otimes \frac{{e}^{\mathrm{iA\tau m}/M}\ue89e{\mathrm{Be}}^{\mathrm{iA\tau m}/M}}{{\alpha}_{B}},\text{}\ue89e\phantom{\rule{4.4em}{4.4ex}}\ue89e{\alpha}_{B}\ge \uf605B\uf606.& \left(35\right)\end{array}$
Then e^{−iHt }may be approximated to error ϵ using
$1.\ue89e\phantom{\rule{0.8em}{0.8ex}}\ue89e\mathrm{Queries}\ue89e\phantom{\rule{0.8em}{0.8ex}}\ue89e\mathrm{to}\ue89e\phantom{\rule{0.8em}{0.8ex}}\ue89e\mathrm{HAM}\ue89e\text{}\ue89eT\ue89e\text{:}\ue89e\phantom{\rule{1.1em}{1.1ex}}\ue89e\ue52e\ue8a0\left({\alpha}_{B}\ue89et\ue89e\frac{\mathrm{log}\ue8a0\left({\alpha}_{B}\ue89et/\u03f5\right)}{\mathrm{log}\ue89e\phantom{\rule{0.3em}{0.3ex}}\ue89e\mathrm{log}\ue8a0\left({\alpha}_{B}\ue89et/\u03f5\right)}\right).\text{}\ue89e2.\ue89e\phantom{\rule{0.8em}{0.8ex}}\ue89e\mathrm{Queries}\ue89e\phantom{\rule{0.8em}{0.8ex}}\ue89e\mathrm{to}\ue89e\phantom{\rule{0.8em}{0.8ex}}\ue89e{e}^{\mathrm{iA}\ue89e\phantom{\rule{0.3em}{0.3ex}}\ue89e\tau}\ue89e\text{:}\ue89e\phantom{\rule{0.8em}{0.8ex}}\ue89e\ue52e\ue8a0\left({\alpha}_{B}\ue89et\right).\text{}\ue89e3.\ue89e\phantom{\rule{0.8em}{0.8ex}}\ue89e\mathrm{Qubits}\ue89e\text{:}\ue89e\phantom{\rule{0.8em}{0.8ex}}\ue89e{n}_{s}+\ue52e\ue8a0\left({n}_{a}+\mathrm{log}\ue8a0\left(M\right)\right).\text{}\ue89e4.\ue89e\phantom{\rule{0.8em}{0.8ex}}\ue89e\mathrm{Primitive}\ue89e\phantom{\rule{0.8em}{0.8ex}}\ue89e\mathrm{gates}\ue89e\text{:}\ue89e\phantom{\rule{0.8em}{0.8ex}}\ue89e\ue52e\ue8a0\left(\left({n}_{a}+\mathrm{log}\ue8a0\left(M\right)\right)\ue89e{\alpha}_{B}\ue89et\ue89e\frac{\mathrm{log}\ue8a0\left({\alpha}_{B}\ue89et/\u03f5\right)}{\mathrm{log}\ue89e\phantom{\rule{0.3em}{0.3ex}}\ue89e\mathrm{log}\ue8a0\left({\alpha}_{B}\ue89et/\u03f5\right)}\right).$
Proof. The number of segments L=(α_{B}t) in Eq. (34) is determined by a normalization constant α_{B}≥max_{t}∥H_{I}(t)=∥B∥, and so each segment is of size =(α_{B}^{−1}). After rescaling ϵ→ϵ/L by the number of segments, the total query complexity is obtained from the truncated Dyson series algorithm Theorem 1. Using the facts max_{s }∥H_{I}(s)∥≤∥B∥, ∥H∥=∥[A,B]∥≤2∥A∥∥B∥, and =(α_{B}^{−1}), it suffices to choose
$\begin{array}{cc}M=\ue52e\ue8a0\left(\frac{{\alpha}_{B}\ue89et\ue89e\phantom{\rule{0.3em}{0.3ex}}\ue89e{\tau}^{2}}{\u03f5}\ue89e\left(\frac{\u3008\uf605\stackrel{.}{H}\uf606\u3009}{{\alpha}_{B}}+\frac{{\mathrm{max}}_{s}\ue89e{\uf605H\ue8a0\left(s\right)\uf606}^{2}}{{\alpha}_{B}^{2}}\right)\right)=\ue52e\ue8a0\left(\frac{t}{\u03f5}\ue89e\frac{\uf605A\uf606}{{\alpha}_{B}}\right).& \left(36\right)\end{array}$
Each segment is simulated
$\ue52e\ue8a0\left(\frac{\mathrm{log}\ue8a0\left({\alpha}_{B}\ue89et/\u03f5\right)}{\mathrm{log}\ue89e\phantom{\rule{0.3em}{0.3ex}}\ue89e\mathrm{log}\ue8a0\left({\alpha}_{B}\ue89et/\u03f5\right)}\right)$
queries to HAMT, thus simulation for the full duration has query complexity
$\phantom{\rule{0.3em}{0.3ex}}\ue89e\ue52e\ue89e\left({\alpha}_{B}\ue89et\ue89e\frac{\mathrm{log}\ue8a0\left({\alpha}_{B}\ue89et/\u03f5\right)}{\mathrm{log}\ue89e\phantom{\rule{0.3em}{0.3ex}}\ue89e\mathrm{log}\ue8a0\left({\alpha}_{B}\ue89et/\u03f5\right)}\right).$
A. Comparison with Simulation of TimeIndependent Hamiltonians in the Schrödinger Picture
The cost of simulation in the interaction picture using the truncated Dyson series is now compared with stateofart simulation in the Schrödinger picture with timeindependent Hamiltonians using the truncated Taylor series approach outlined in section X A. Up to logarithmic factors, this comparison is valid as the truncated Taylor series algorithm cost differs from optimal algorithms by only logarithmic factors.
For any Hamiltonian H=A+B, let us assume access to the oracles
$\begin{array}{cc}\left({\u30080\uf604}_{a}\otimes {I}_{s}\right)\ue89e{O}_{A}\ue8a0\left({\uf6030\u3009}_{a}\otimes {I}_{s}\right)=\frac{A}{{\alpha}_{A}},\text{}\ue89e\left({\u30080\uf604}_{a}\otimes {I}_{s}\right)\ue89e{O}_{B}\ue8a0\left({\uf6030\u3009}_{a}\otimes {I}_{s}\right)=\frac{B}{{\alpha}_{B}},& \left(37\right)\end{array}$
which have gate complexity C_{A},C_{B }respectively. The gate complexity of timeindependent simulation e^{−i(A+B)t }is then
$\begin{array}{cc}\begin{array}{c}{C}_{\mathrm{TTS}}=\ue89e\ue52e\ue8a0\left(\left({C}_{A}+{C}_{B}+{n}_{a}\right)\ue89e\left({\alpha}_{A}+{\alpha}_{B}\right)\ue89et\ue89e\frac{\mathrm{log}\ue8a0\left(\left({\alpha}_{A}+{\alpha}_{B}\right)\ue89et/\u03f5\right)}{\mathrm{log}\ue89e\phantom{\rule{0.3em}{0.3ex}}\ue89e\mathrm{log}\ue8a0\left(\left({\alpha}_{A}+{\alpha}_{B}\right)\ue89et/\u03f5\right)}\right)\\ =\ue89e\ue52e\ue8a0\left(\left({C}_{A}+{C}_{B}+{n}_{a}\right)\ue89e\left({\alpha}_{A}+{\alpha}_{B}\right)\ue89e\mathrm{tpoly}\ue89e\phantom{\rule{0.3em}{0.3ex}}\ue89e\mathrm{log}\ue8a0\left({\alpha}_{A},{\alpha}_{B},t,\u03f5\right)\right).\end{array}& \left(38\right)\end{array}$
In the interactionpicture, one can prove the following theorem
Theorem 3 (Gate complexity of Hamiltonian simulation in the interaction picture). For any Hamiltonian H=A+B, let
 1. C_{B }be the gate cost of implementing the oracle
$\left({\u30080\uf604}_{a}\otimes {I}_{s}\right)\ue89e{O}_{B}\ue8a0\left({\uf6030\u3009}_{a}\otimes {I}_{s}\right)=\frac{B}{{\alpha}_{B}}.$

 2. C_{e}_{−iAt[ϵ]}=(tpolylog(t,ϵ)) be the gate cost of approximating e^{−iAt }to error ϵ.
Then the total gate complexity of simulating e^{−iHt }is
C_{TDS}=((C_{B}+C_{e}_{iAt}[ϵ]+n_{α})α_{B}tpolylog(∥A∥,α_{B},t,ϵ)). (39)
Proof. let C_{HAMT}ϵ be the gate complexity of approximating HAMT in Eq. (35) to error ϵ. Then from Lemma 3, the total gate complexity of simulating e^{−iHt}=(e^{−iAt}[e^{−i∫}^{0}^{t}^{H}^{I}^{(s)ds}])^{t/T}, T=(α_{B}^{−1}), is
$\begin{array}{cc}\ue52e\ue8a0\left(\left({C}_{\mathrm{HAM}\ue89e\text{}\ue89eT}\ue8a0\left[\frac{{\alpha}_{B}\ue89et}{\u03f5}\right]+{n}_{a}+\mathrm{log}\ue8a0\left(\frac{\uf605A\uf606\ue89et}{{\mathrm{\u03f5\alpha}}_{B}}\right)\right)\ue89e\left({\alpha}_{B}\ue89et\ue89e\phantom{\rule{0.3em}{0.3ex}}\ue89e\frac{\mathrm{log}\ue8a0\left({\alpha}_{B}\ue89et/\u03f5\right)}{\mathrm{log}\ue89e\phantom{\rule{0.3em}{0.3ex}}\ue89e\mathrm{log}\ue8a0\left({\alpha}_{B}\ue89et/\u03f5\right)}\right)\right).& \left(40\right)\end{array}$
One possible decomposition of HAMT is
$\begin{array}{cc}\phantom{\rule{4.4em}{4.4ex}}\ue89e\mathrm{HAM}\ue89e\text{}\ue89eT=& \left(41\right)\\ \left(\sum _{m=0}^{M1}\ue89e\uf603m\u3009\ue89e{\u3008m\uf604}_{d}\otimes {I}_{a}\otimes {e}^{\mathrm{iA}\ue89e\phantom{\rule{0.3em}{0.3ex}}\ue89e\tau \ue89e\phantom{\rule{0.3em}{0.3ex}}\ue89em/M}\right)\xb7\left({I}_{d}\otimes {O}_{B}\right)\ue89e\phantom{\rule{0.em}{0.ex}}\xb7\left(\sum _{m=0}^{M1}\ue89e\uf603m\u3009\ue89e{\u3008m\uf604}_{d}\otimes {I}_{a}\otimes {e}^{\mathrm{iA}\ue89e\phantom{\rule{0.3em}{0.3ex}}\ue89e\tau \ue89e\phantom{\rule{0.3em}{0.3ex}}\ue89em/M}\right).& \left(42\right)\end{array}$
Using C_{e}_{−iAt}[ϵ]=(t polylog(t, ϵ)), synthesis of Σ_{m=0}^{M−1}mm_{d}⊗e^{iArm/M }by exponentiated power of controllede^{iAt/M},e^{iA2t/M},e^{iA4t/M}, . . . , e^{iA2}^{└log}^{2}^{(M)┘}^{T/M}, has cost C_{e}_{iAt2j/M}[ϵ/log (M)]=(C_{e}_{iAt}[ϵ/log (M)]). Thus C_{HAMT}[ϵ]=(C_{B}+C_{e}_{iAt}[ϵ/log (M))]). □
From comparing Eqs. (38) and (39), one may iminediately state sufficient criteria for when simulation in the interaction picture is advantageous over simulation in the Schrödinger picture.
 1. The upper bound on the spectral norms α_{A}≥∥A∥, α_{B}≥∥A∥ of the encoding in Eq. (37) satisfy α_{A}»α_{B}. Generally speaking, this is correlated with term A representing fast dynamics ∥A∥»∥B∥.
 2. The gate complexity of timeevolution by A alone for time =(α_{B}^{−1}) is comparable to that of synthesizing the oracle O_{B}, that is C_{e}_{iAα}_{B}_{−1}=(C_{B}).
Note that satisfying condition (2) depends strongly on the structure of A, B. For instance, a simulation of A for time =(α_{B}^{−1}) using generic timeindependent techniques has gate complexity (C_{A}_{α}_{A}/α_{B}). As the interest is in the case ∥A∥/∥B∥»1, this quantity could be large and scale poorly with the problem size. One sufficient possibility is the very strong assumption that e^{−iAt }is cheap and can be fastforwarded., such that the gate complexity C_{e}_{−iAt}[ϵ]=(polylog(t,ϵ)) is constant up to logarithmic factor. This turns out to be a reasonable assumption in the application to next consider.
V. Application to the Hubbard Model with LongRanged Interactions
One can now apply the technology developed in Section III and Section IV for Hamiltonian simulation in the interaction picture to physical problems of practical interest. Here, the focus is on the Hubbbard model in ddimensions with N lattice sites subject to local disorder and translationalinvariant twobody couplings that may be longranged in general. One can perform a gate complexity comparison with simulation by timeindependent techniques, and later in Section V A, this model is specialized to that of quantum chemistry simulations in the planewave and dual basis.
The Hubbard Hamiltonian considered has the form H=T+U+V, where T is the kinetic energy hopping operator, U is the local singlesite potential, and V is a symmetric translationallyinvariant twobody density coupling term between opposite spins. In the dual basis, H is expressed in terms of singlesite Fermionic creation and annihilation operators {α_{{right arrow over (x)}σ}, α_{{right arrow over (y)}σ′}}={α_{{right arrow over (x)}σ}^{554 }, α_{{right arrow over (y)}σ′}^{554 }}=0, {α_{{right arrow over (x)}σ}, α_{{right arrow over (y)}σ′}^{554 }}=δ_{{right arrow over (x)}{right arrow over (y)}}δ_{σσ′}, and the number operator n_{{right arrow over (x)}σ}=α_{{right arrow over (x)}σ}^{554 }α_{{right arrow over (x)}σ}. The subscript {right arrow over (x)} ∈ [−N^{1/d},N^{1/d}]^{d }indexes one of N lattice sites in d dimensions, and σ ∈ {−1,1} is a spin½ index. Explicitly,
$\begin{array}{cc}H=\sum _{\overrightarrow{x},\overrightarrow{y},\sigma}\ue89e\phantom{\rule{0.3em}{0.3ex}}\ue89eT\ue8a0\left(\overrightarrow{x}\overrightarrow{y}\right)\ue89e{a}_{\overrightarrow{x}\ue89e\sigma}^{\u2020}\ue89e{a}_{\overrightarrow{y}\ue89e\sigma}+\sum _{\overrightarrow{x},\sigma}\ue89eU\ue8a0\left(\overrightarrow{x},\sigma \right)\ue89e{n}_{\overrightarrow{x}\ue89e\sigma}+\sum _{\left(\overrightarrow{x},\sigma \right)\ne \left(\overrightarrow{y},{\sigma}^{\prime}\right)}\ue89eV\ue8a0\left(\overrightarrow{x}\overrightarrow{y}\right)\ue89e{n}_{\overrightarrow{x}\ue89e\sigma}\ue89e{n}_{\overrightarrow{y}\ue89e{\sigma}^{\prime}},& \left(43\right)\end{array}$
where the coefficients T({right arrow over (s)}), U({right arrow over (s)},σ), V({right arrow over (s)}) are real functions of the lattice index {right arrow over (s)} ∈ [−N^{1/d},N^{1/d}]^{d}.
Further simplification of Eq. (43) is possible as the kinetic energy operator is diagonal in the planewave basis. This basis related to the dual basis by a unitary rotation FFFT, an acronym for ‘FastFermionicFourierTransform’ that implements a Fourier transform over the lattice site indices, resulting in Fermionic creation and annihilation operators c_{{right arrow over (p)}σ}^{554 }, c_{{right arrow over (p)}σ}.
$\begin{array}{cc}{c}_{\overrightarrow{p}\ue89e\sigma}=\frac{1}{\sqrt{N}}\ue89e\sum _{\overrightarrow{x}}\ue89e{a}_{\overrightarrow{x}\ue89e\sigma}\ue89e{e}^{i\ue89e\phantom{\rule{0.3em}{0.3ex}}\ue89e2\ue89e\pi \ue89e\overrightarrow{p}\xb7\overrightarrow{x}/{N}^{1/d}}={\mathrm{FFFT}}^{\u2020}\ue89e{a}_{\overrightarrow{p}\ue89e\sigma}\ue89e\mathrm{FFFT},& \left(44\right)\\ {c}_{\overrightarrow{p}\ue89e\sigma}^{\u2020}=\frac{1}{\sqrt{N}}\ue89e\sum _{\overrightarrow{x}}\ue89e{a}_{\overrightarrow{x}\ue89e\sigma}^{\u2020}\ue89e{e}^{i\ue89e\phantom{\rule{0.3em}{0.3ex}}\ue89e2\ue89e\pi \ue89e\overrightarrow{p}\xb7\overrightarrow{x}/{N}^{1/d}}={\mathrm{FFFT}}^{\u2020}\ue89e{a}_{\overrightarrow{p}\ue89e\sigma}^{\u2020}\ue89e\mathrm{FFFT},& \left(45\right)\end{array}$
By substituting the Fourier transform of the kinetic term
$T\ue8a0\left(\overrightarrow{s}\right)=\frac{1}{N}\ue89e{\sum}_{\overrightarrow{p}}\ue89e\stackrel{~}{T}\ue8a0\left(\overrightarrow{p}\right)\ue89e{e}^{i\ue89e\phantom{\rule{0.3em}{0.3ex}}\ue89e2\ue89e\pi \ue89e\overrightarrow{p}\xb7\overrightarrow{s}/{N}^{1/d}},$
an equivalent expression for the Hubbard Hamiltonian is
$\begin{array}{cc}H={\mathrm{FFFT}}^{\u2020}\xb7\left(\sum _{\overrightarrow{x},\sigma}\ue89e\phantom{\rule{0.3em}{0.3ex}}\ue89e\stackrel{~}{T}\ue8a0\left(\overrightarrow{x}\right)\ue89e{n}_{\overrightarrow{x}\ue89e\sigma}\right)\xb7\mathrm{FFFT}+\sum _{\overrightarrow{x},\sigma}\ue89eU\ue8a0\left(\overrightarrow{x},\sigma \right)\ue89e{n}_{\overrightarrow{x}\ue89e\sigma}+\sum _{\left(\overrightarrow{x},\sigma \right)\ne \left(\overrightarrow{y},{\sigma}^{\prime}\right)}\ue89eV\ue8a0\left(\overrightarrow{x}\overrightarrow{y}\right)\ue89e{n}_{\overrightarrow{x}\ue89e\sigma}\ue89e{n}_{\overrightarrow{y}\ue89e{\sigma}^{\prime}},& \left(46\right)\end{array}$
where each term is now diagonal in their respective bases.
A simulation of this Hamiltonian on a qubit quantum computer requires a map from its Fermionic operators to spin operators. One possibility is the JordanWigner transformation, which requires some map from Fermionic indices to spin indices, such as
$f\ue8a0\left(\overrightarrow{x},\sigma \right)=N\ue89e\frac{1\sigma}{2}+\left({\sum}_{j=0}^{d1}\ue89e\phantom{\rule{0.3em}{0.3ex}}\ue89e{\overrightarrow{x}}_{j}\ue89e{N}^{j/d}\right).$
Subsequently, one can replace
$\begin{array}{cc}{a}_{\overrightarrow{x}\ue89e\sigma}^{\u2020}>\frac{1}{2}\ue89e\left({X}_{f\ue8a0\left(\overrightarrow{x},\sigma \right)}{\mathrm{iY}}_{f\ue8a0\left(\overrightarrow{x},\sigma \right)}\right)\ue89e\underset{j=0}{\stackrel{f\ue8a0\left(\overrightarrow{x},\sigma \right)1}{\otimes}}\ue89e{Z}_{j},\text{}\ue89e{a}_{\overrightarrow{x}\ue89e\sigma}>\frac{1}{2}\ue89e\left({X}_{f\ue8a0\left(\overrightarrow{x},\sigma \right)}+{\mathrm{iY}}_{f\ue8a0\left(\overrightarrow{x},\sigma \right)}\right)\ue89e\underset{j=0}{\stackrel{f\ue8a0\left(\overrightarrow{x},\sigma \right)1}{\otimes}}\ue89e{Z}_{j}.& \left(47\right)\end{array}$
Note the very useful property where number operators map to singlesite spin operators n_{{right arrow over (x)}σ}→½(I−Z_{∫({right arrow over (x)},σ)}) (under this encoding. Moreover, FFFT, can be implemented using (N log (N)) primitive quantum gates in the JordanWigner representation.
Now evaluate the worstcase gatecomplexity for timeevolution by H in the Schrödinger picture. As an example of stateofart using the truncated Taylor series approach in Section X A, e^{−i(T+U+V)t }may be simulated using ((α_{T}+α_{U}+α_{V})t log (1/ϵ)) queries to oracles that encode T, U, and V.
$\begin{array}{cc}{(\u30080\uf604}_{a}\otimes {I}_{s})\ue89e{O}_{T}\ue8a0\left({\uf6030\u3009}_{a}\otimes {I}_{s}\right)=\frac{T}{{\alpha}_{T}},\text{}\ue89e{(\u30080\uf604}_{a}\otimes {I}_{s})\ue89e{O}_{U}\ue8a0\left({\uf6030\u3009}_{a}\otimes {I}_{s}\right)=\frac{U}{{\alpha}_{U}},\text{}\ue89e{(\u30080\uf604}_{a}\otimes {I}_{s})\ue89e{O}_{V}\ue8a0\left({\uf6030\u3009}_{a}\otimes {I}_{s}\right)=\frac{V}{{\alpha}_{V}}.& \left(48\right)\end{array}$
The cost of simulation depends strongly on the coefficients {tilde over (T)} ({right arrow over (p)}), U ({right arrow over (x)}), V ({right arrow over (x)}). The most straightforward approach synthesizes these oracles using the linearcombination of unitaries outline in Eq. (6). For instance, O_{T}=(PREP_{T}^{†}⊗FFFT^{†})·SEL_{T}·(PREP_{T}⊗FFFT), where
$\begin{array}{cc}{\mathrm{PREP}}_{T}\ue89e{\uf6030\u3009}_{a}=\sum _{\overrightarrow{p},\sigma}\ue89e\phantom{\rule{0.3em}{0.3ex}}\ue89e\sqrt{\frac{\phantom{\rule{0.3em}{0.3ex}}\ue89e\stackrel{~}{T}\ue8a0\left(\overrightarrow{p}\right)}{{\alpha}_{T}}}\ue89e{\uf603\overrightarrow{p},\sigma \u3009}_{a},\text{}\ue89e{\mathrm{SEL}}_{T}=\sum _{\overrightarrow{p},\sigma}\ue89e\uf603\overrightarrow{p},\sigma \u3009\ue89e{\u3008\overrightarrow{p},\sigma \uf604}_{a}\otimes {n}_{\overrightarrow{p}\ue89e\sigma},\text{}\ue89e{\alpha}_{T}=\sum _{\overrightarrow{p},\sigma}\ue89e\uf603\stackrel{~}{T}\ue8a0\left(\overrightarrow{p}\right)\uf604.& \left(49\right)\end{array}$
and similarly for U and V. As there are (N) distinct coefficients in the worstcase, each of PREP_{T,U,V }costs (N). As V has (N^{2}) terms, SEL_{V }has the largest cost of (N^{2}). Thus overall gate complexity is (N^{2}(α_{T}+α_{U}+α_{V})t log (1/ϵ)). As there are (N^{2}) coefficients, max{α_{T}, α_{U}, α_{V}}=(α_{T}+N^{2}), and so the cost of simulation is
(N^{2}(α_{T}+N^{2})t log (1/ϵ)). (50)
The worstcase gatecomplexity may be substantially improved by instead simulating H in the interaction picture using the truncated Dyson series algorithm in Section III. The key idea is to simulate in the rotating frame of the interactions e^{−i(U+V)t}, where the Hamiltonian becomes a timedependent H_{I}(t)=e^{i(U+V)t}_{Te}^{−i(U+V)t}. Using the same oracle O_{T }in Eq. (48) for the kinetic term, the cost of timeevolution e^{−i(T+U+V)t }by this technique is given by Eq. (39):
$\begin{array}{cc}\begin{array}{c}{C}_{\mathrm{TDS}}=\ue89e\ue52e(\left({C}_{T}+{C}_{{e}^{i\ue8a0\left(U+V\right)/{\alpha}_{T+{n}_{a}}}}\ue8a0\left[\u03f5\right]\right)\\ \ue89e{\alpha}_{T}\ue89e\mathrm{tpoly}\ue89e\phantom{\rule{0.3em}{0.3ex}}\ue89e\mathrm{log}\ue8a0\left(\uf605U+V\uf606,{\alpha}_{T},t,\u03f5\right))\\ =\ue89e\ue52e(\left(N+{C}_{{e}^{i\ue8a0\left(U+V\right)/{\alpha}_{T}}}\ue8a0\left[\u03f5\right]\right)\\ \ue89e{\alpha}_{T}\ue89e\mathrm{tpoly}\ue89e\phantom{\rule{0.3em}{0.3ex}}\ue89e\mathrm{log}\ue8a0\left(N,\uf605U\uf606,\uf605V\uf606,{\alpha}_{T},t,\u03f5\right)).\end{array}& \left(51\right)\end{array}$
All that remains is to bound the cost of timeevolution by the term C_{e}_{i(U+V)/α}_{T}[ϵ]. Using the fact that this is diagonal in the Pauli Z basis, the Hamiltonian may be fastforwarded and so has cost that is independent of the evolution time. Thus the most straightforward approach decomposes
$\begin{array}{cc}{e}^{i\ue8a0\left(U+V\right)\ue89et}=\left(\prod _{\overrightarrow{x},\sigma}\ue89e\phantom{\rule{0.3em}{0.3ex}}\ue89e{e}^{\mathrm{iU}\ue8a0\left(\overrightarrow{x},\sigma \right)\ue89e{n}_{\overrightarrow{x}\ue89e\sigma}\ue89et}\right)\ue89e\left(\prod _{\left(\overrightarrow{x},\sigma \right)\ne \left(\overrightarrow{y},\sigma \right)}\ue89e{e}^{\mathrm{iV}\ue8a0\left(\overrightarrow{x}\overrightarrow{y}\right)\ue89e{n}_{\overrightarrow{x}\ue89e\sigma}\ue89e{n}_{\overrightarrow{y}\ue89e{\sigma}^{\prime}}\ue89et}\right).& \left(52\right)\end{array}$
There are (N^{2}) exponentials, and so C_{e}_{i(U+V)t}[ϵ]=(N^{2}). Compared to Eq. (50), one can already see an improvement by a factor (N) in cases where the kinetic energy is extensive, that is α_{T}=(N), so C_{TDS}=(N^{2}α_{T}tpolylog(N, α_{T}, t, ϵ)).
A further improvement to
C_{TDS}=(Nα_{T}tpolylog(N, α_{T}, t, ϵ)) (53)
is possible by a more creative evaluation of the gate complexity of e^{i(U+V)t }to reduce its cost from (N^{2}) to (N log (N)). Clearly, C_{e}_{iUt}=(N) with N commuting terms poses no problem. The difficulty lies in constructing timeevolution by the twobody term e^{iVt }such that C_{e}_{iVt}=(N log N). As V is a sum of (N^{2}) commuting terms, a gate cost (N^{2}) appear unavoidable. However, this may be reduced by exploiting the translation symmetry of its coefficients with a discrete Fourier transform. As V({right arrow over (x)})=V(−{right arrow over (x)}) is real and symmetric, its discrete Fourier transform {tilde over (V)}({right arrow over (k)})=Σ_{{right arrow over (x)}}V({right arrow over (x)})e^{i2π{right arrow over (x)}·{right arrow over (k)}/N}^{1/d }only has real coefficients. Let one rewrite V from Eq. (43) as
$\begin{array}{cc}V=\sum _{\left(\overrightarrow{x},\sigma \right)\ne \left(\overrightarrow{y},\sigma \right)}\ue89eV\ue8a0\left(\overrightarrow{x}\overrightarrow{y}\right)\ue89e{n}_{\overrightarrow{x}\ue89e\sigma}\ue89e{n}_{\overrightarrow{y}\ue89e{\sigma}^{\prime}}=\sum _{\left(\overrightarrow{x},\sigma \right)\ne \left(\overrightarrow{y},\sigma \right)}\ue89e\frac{1}{N}\ue89e\sum _{\overrightarrow{k}}\ue89e\stackrel{~}{V}\ue8a0\left(\overrightarrow{k}\right)\ue89e{e}^{{\phantom{\rule{0.3em}{0.3ex}}}^{i\ue89e\phantom{\rule{0.3em}{0.3ex}}\ue89e2\ue89e\pi \ue8a0\left(\overrightarrow{x}\overrightarrow{y}\right)\xb7\overrightarrow{k}/N}}\ue89e{n}_{\overrightarrow{x}\ue89e\sigma}\ue89e{n}_{\overrightarrow{y}\ue89e{\sigma}^{\prime}}=\sum _{\overrightarrow{k}}\ue89e\frac{\stackrel{~}{V}\ue8a0\left(\overrightarrow{k}\right)}{N}\ue89e\left(\sum _{\left(\overrightarrow{x},\sigma \right),\left(\overrightarrow{y},{\sigma}^{\prime}\right)}\ue89e{e}^{{\phantom{\rule{0.3em}{0.3ex}}}^{i\ue89e\phantom{\rule{0.3em}{0.3ex}}\ue89e2\ue89e\pi \ue8a0\left(\overrightarrow{x}\overrightarrow{y}\right)\xb7\overrightarrow{k}/N}}\ue89e{n}_{\overrightarrow{x},\sigma}\ue89e{n}_{\overrightarrow{y},{\sigma}^{\prime}}\sum _{\overrightarrow{x},\sigma}\ue89e{n}_{\overrightarrow{x}\ue89e\sigma}\right)=\sum _{\overrightarrow{k}}\ue89e\stackrel{~}{V}\ue8a0\left(\overrightarrow{k}\right)\ue89e\underset{{\stackrel{~}{\chi}}_{k}}{\underset{\uf613}{\left(\frac{1}{\sqrt{N}}\ue89e\sum _{\overrightarrow{x}}\ue89e{e}^{{\phantom{\rule{0.3em}{0.3ex}}}^{i\ue89e\phantom{\rule{0.3em}{0.3ex}}\ue89e2\ue89e\pi \ue89e\overrightarrow{x}\xb7\overrightarrow{k}/N}}\ue89e\sum _{\sigma}\ue89e\phantom{\rule{0.3em}{0.3ex}}\ue89e{n}_{\overrightarrow{x},\sigma}\right)}}\ue89e\underset{{\stackrel{~}{\chi}}_{k}^{\u2020}}{\underset{\uf613}{\left(\frac{1}{\sqrt{N}}\ue89e\sum _{\overrightarrow{y}}\ue89e{e}^{{\phantom{\rule{0.3em}{0.3ex}}}^{i\ue89e\phantom{\rule{0.3em}{0.3ex}}\ue89e2\ue89e\pi \ue89e\overrightarrow{y}\xb7\overrightarrow{k}/N}}\ue89e\sum _{{\sigma}^{\prime}}\ue89e\phantom{\rule{0.3em}{0.3ex}}\ue89e{n}_{\overrightarrow{y},{\sigma}^{\prime}}\right)}}\sum _{\overrightarrow{p},\sigma}\ue89e\left(\sum _{\overrightarrow{k}}\ue89e\stackrel{~}{V}\ue8a0\left(\stackrel{~}{k}\right)\right)\ue89e{n}_{\overrightarrow{p},\sigma}.& \left(54\right)\end{array}$
The strategy for implementing e^{−iVt }is based on the following observation: Suppose one had a unitary oracle O_{{right arrow over (A)}}j0_{o}0_{garb}=jA_{j}_{o}g(j)_{garbage }that on input j ∈^{dim[{right arrow over (A)}]}, outputs on the l quoit o register, the value of the j^{th }element of some complex vector {right arrow over (A)}, together with some garbage state g(j)_{garb }of lesser interest required to make the operation reversible. One may then perform a phase rotation that depends on A_{j }as follows:
$\begin{array}{cc}\uf603j\u3009\ue89e{\uf6030\u3009}_{o}\ue89e{\uf6030\u3009}_{\mathrm{garb}}\ue89e\uf6030\u3009\ue89e\underset{{O}_{\overrightarrow{A}}}{\to}\ue89e\uf603j\u3009\ue89e{\uf603{A}_{j}\u3009}_{o}\ue89e{\uf603g\ue8a0\left(j\right)\u3009}_{\mathrm{garb}}\ue89e\uf6030\u3009\ue89e\underset{\mathrm{PHASE}}{\to}\ue89e{e}^{{\mathrm{iA}}_{j}\ue89et}\ue89e\uf603j\u3009\ue89e{\uf603{A}_{j}\u3009}_{o}\ue89e{\uf603g\ue8a0\left(j\right)\u3009}_{\mathrm{garb}}\ue89e\uf6030\u3009\ue89e\underset{{O}_{\overrightarrow{A}}^{\u2020}}{\to}\ue89e{e}^{{\mathrm{iA}}_{j}\ue89et}\ue89e\uf603j\u3009\ue89e{\uf6030\u3009}_{o}\ue89e{\uf6030\u3009}_{\mathrm{garb}}\ue89e\uf6030\u3009.& \left(55\right)\end{array}$
If A_{j }were represented in binary, say, A_{j}=Σ_{k=0}^{l−1}q_{k}2^{−k}, PHASE could be implemented using (l) controlledphase 00⊗I+11⊗e^{−it2}^{−k}^{Z }rotations.
Thus, one can construct a unitary O_{V,binary }with the property that
$\begin{array}{cc}{O}_{V,\mathrm{binary}}\ue8a0\left(\underset{\overrightarrow{x},\sigma}{\otimes}\ue89e\uf603{n}_{\overrightarrow{x},\sigma}\u3009\right)\ue89e\uf6030\u3009\ue89e{\uf6030\u3009}_{\mathrm{garb}}=\left(\underset{\overrightarrow{x},\sigma}{\otimes}\ue89e\uf603{n}_{\overrightarrow{x},\sigma}\u3009\right)\ue89e\uf603f\ue8a0\left(\overrightarrow{n}\right)\u3009\ue89e{\uf603g\ue8a0\left(\overrightarrow{n}\right)\u3009}_{\mathrm{garb}},\text{}\ue89ef\ue8a0\left(\overrightarrow{n}\right)=\sum _{\left(\overrightarrow{x},\sigma \right)\ne \left(\overrightarrow{y},{\sigma}^{\prime}\right)}\ue89eV\ue89e\left(\overrightarrow{x}\overrightarrow{y}\right)\ue89e{n}_{\overrightarrow{p}\ue89e\sigma}\ue89e{n}_{\overrightarrow{q}\ue89e{\sigma}^{\prime}},& \left(56\right)\end{array}$
where the value ƒ({right arrow over (n)}) is encoded in l=(log(1/ϵ)) bits. This is implemented by the following sequence, where one has omitted the garbage register for clarity.
$\begin{array}{cc}\left(\underset{\overrightarrow{x},\sigma}{\otimes}\ue89e\uf603{n}_{\overrightarrow{x},\sigma}\u3009\right)\ue89e\uf6030\u3009\ue89e\underset{\mathrm{ADD}}{\to}\ue89e\underset{\overrightarrow{x}}{\otimes}\ue89e\uf603\sum _{\sigma}\ue89e{n}_{\overrightarrow{x},\sigma}\u3009\ue89e\underset{\mathrm{FFT}}{\to}\ue89e\underset{\overrightarrow{k}}{\otimes}\ue89e\uf603{\stackrel{~}{\chi}}_{\overrightarrow{k}}\u3009\ue89e\underset{{\uf603\xb7\uf604}^{2}}{\to}\ue89e\underset{\overrightarrow{k}}{\otimes}\ue89e\uf603{\uf603{\stackrel{~}{\chi}}_{\overrightarrow{k}}\uf604}^{2}\u3009\ue89e\underset{\times {V}_{k}}{\to}\ue89e\underset{\overrightarrow{k}}{\otimes}\ue89e\uf603V\ue8a0\left(\overrightarrow{k}\right)\ue89e{\uf603{\stackrel{~}{\chi}}_{\overrightarrow{k}}\uf604}^{2}\u3009,& \left(57\right)\end{array}$
The cost of O_{V,binary }may be expressed in term of the four standard reversible arithmetical operations, addition, subtraction, division, and multiplication, which each cost (poly(l)) primitive gates. The first steps ADD adds (N) pairs of two bits n_{{right arrow over (x)},σ=1 }+n_{{right arrow over (x)},σ=1 }and costs (N) arithmetic operations. The second step FFT is a ddimensional FastFourierTransform on (N) binary numbers and requires (N log (N)) arithmetic operations. The third step computes the absolutevaluesquared of (N) binary numbers, and uses (N) arithmetic operations. The last step multiplies each {tilde over (x)}_{k}^{2 }with the corresponding V_{k}, and costs (N) arithmetic operations. This last step may actually be avoided by rescaling the time parameter e^{−iA}^{j}^{t}→e^{−iA}^{j}^{V}^{k}^{t }in Eq. (55). Thus the total cost of O_{V,binary }is (N log (N) poly(l))=(N log(1/ϵ)). Using one query to O_{V,binary}, O_{V,binary}^{†}, and (N log(1/ϵ)) primitive quantum gates, one may thus apply e^{−iVt }with a phase error (ϵ) for a fixed value of t.
A. Application to Quantum Chemistry in the PlaneWave Basis
The Hamiltonian that generates timeevolution for a state φ(t) of interacting electrons in d=3 dimension consists of three operators: the electron kinetic energy T, the electronnuclei potential energy U, and the electronelectron potential energy V. It has been demonstrated that this electronic structure Hamiltonian is a special case of the general Hubbard Hamiltonian of Eq. (46). In the planewave basis,
$\begin{array}{cc}\phantom{\rule{4.4em}{4.4ex}}\ue89ei\ue89e\phantom{\rule{0.3em}{0.3ex}}\ue89e\partial \uf603\psi \ue8a0\left(t\right)\u3009=H\ue89e\uf603\psi \ue8a0\left(t\right)\u3009\ue89e\text{}\ue89e{H}_{P}=\frac{1}{2}\ue89e\sum _{\stackrel{>}{p},\sigma}\ue89e\phantom{\rule{0.3em}{0.3ex}}\ue89e{\uf603{\stackrel{>}{k}}_{\stackrel{>}{p}}\uf604}^{2}\ue89e{c}_{\stackrel{>}{p},\sigma}^{\u2020}\ue89e{c}_{\stackrel{>}{p},\sigma}+\frac{4\ue89e\phantom{\rule{0.3em}{0.3ex}}\ue89e\pi}{\Omega}\ue89e\sum _{\stackrel{>}{p}\ne \stackrel{>}{q}\ue89e\text{}\ue89ej,\sigma}\ue89e\phantom{\rule{0.3em}{0.3ex}}\ue89e({\zeta}_{j}\ue89e\frac{{e}^{i\ue89e\phantom{\rule{0.3em}{0.3ex}}\ue89e{\stackrel{>}{k}}_{\stackrel{>}{q}\stackrel{>}{p}}\ue89e{\stackrel{>}{R}}_{j}}}{{\uf603{\stackrel{>}{k}}_{\stackrel{>}{p}\stackrel{>}{q}\ue89e\phantom{\rule{0.3em}{0.3ex}}}\uf604}^{2}}\ue89e{c}_{\stackrel{>}{p},\sigma}^{\u2020}\ue89e{c}_{\stackrel{>}{p},\sigma}+\frac{2\ue89e\phantom{\rule{0.3em}{0.3ex}}\ue89e\pi}{\Omega}\ue89e\sum _{\left(\stackrel{>}{p},\sigma \right)\ne \left(\stackrel{>}{q},{\sigma}^{\prime}\right)\ue89e\text{}\ue89e\stackrel{>}{v}\ne 0}\ue89e\phantom{\rule{0.3em}{0.3ex}}\ue89e\frac{{c}_{\stackrel{>}{p},\sigma}^{\u2020}\ue89e{c}_{\stackrel{>}{q},{\sigma}^{\prime}}^{\u2020}\ue89e{c}_{\stackrel{>}{q}+\stackrel{>}{v},{\sigma}^{\prime}}\ue89e{c}_{\stackrel{>}{p}\stackrel{>}{v},\sigma}}{{\uf603{\stackrel{>}{k}}_{\stackrel{>}{v}}\uf604}^{2}},& \left(58\right)\end{array}$
where k_{{right arrow over (p)}}=2π{right arrow over (p)}/Ω^{1/3}, {right arrow over (p)} ∈ [−N^{1/3},N^{1/3}]^{3}, r_{{right arrow over (p)}}={right arrow over (p)}(Ω/N)^{1/3}, Ω represents the volume of the simulation, and n_{j }is the nuclear charge of the j^{th }nucleus. Whereas T is diagonal here, one may find an alternate basis where U and V are diagonal. This is the dual basis, defined through the unitary transform FFFT of Eq. (44). In this basis, let us define the state φ_{D}(t)=FFFT^{†}φ(t), which evolves under the Hamiltonian H, which is of exactly that of Eq. (46), with coefficients
$\begin{array}{cc}\stackrel{~}{T}\ue8a0\left(\stackrel{>}{p}\right)=\frac{1}{2}\ue89e{\uf603{\stackrel{>}{k}}_{\stackrel{>}{p}}\uf604}^{2}\ue89e\text{}\ue89eU\ue8a0\left(\stackrel{>}{p}\right)=\frac{4\ue89e\phantom{\rule{0.3em}{0.3ex}}\ue89e\pi}{\Omega}\ue89e\sum _{\stackrel{>}{v}\ne 0,\text{}\ue89ej\ue89e\text{}}\ue89e\phantom{\rule{0.3em}{0.3ex}}\ue89e\frac{{\zeta}_{j}\ue89e\phantom{\rule{0.8em}{0.8ex}}\ue89e\mathrm{cos}\ue8a0\left[{\stackrel{>}{k}}_{\stackrel{>}{v}}\xb7\left({\stackrel{>}{R}}_{\stackrel{>}{j}}{\stackrel{>}{r}}_{\stackrel{>}{p}}\right)\right]}{{\uf603{\stackrel{>}{k}}_{\stackrel{>}{v}}\uf604}^{2}},\text{}\ue89eV\ue8a0\left(\stackrel{>}{s}\right)=\frac{2\ue89e\phantom{\rule{0.3em}{0.3ex}}\ue89e\pi}{\Omega}\ue89e\sum _{\stackrel{>}{v}\ne 0}\ue89e\frac{\mathrm{cos}\ue8a0\left[{\stackrel{>}{k}}_{\stackrel{>}{v}}\xb7{\stackrel{>}{r}}_{\stackrel{>}{s}}\right]}{{\uf603{\stackrel{>}{k}}_{\stackrel{>}{v}}\uf604}^{2}}.& \left(59\right)\end{array}$
Thus the cost of timeevolution by the electronic structure Hamiltonian e^{−iHt }using the interaction picture is given by Eq. (53). The only parameter that depends on the problem is the normalization factor
$\begin{array}{cc}{\alpha}_{T}=\sum _{\stackrel{>}{p},\sigma}\ue89e\phantom{\rule{0.3em}{0.3ex}}\ue89e\frac{{\uf603{\stackrel{>}{k}}_{\stackrel{>}{p}}\uf604}^{2}}{2}=\ue52e\ue8a0\left({\int}_{0}^{{N}^{1/3}}\ue89e\frac{{p}^{2}}{{\Omega}^{2/3}}\ue89e\left(4\ue89e\phantom{\rule{0.3em}{0.3ex}}\ue89e\pi \ue89e\phantom{\rule{0.3em}{0.3ex}}\ue89e{p}^{2}\right)\ue89e\mathrm{dp}\right)=\ue52e\ue8a0\left(\frac{{N}^{5/3}}{{\Omega}^{2/3}}\right).\ue89e\phantom{\rule{0.2em}{0.2ex}}& \left(60\right)\end{array}$
Thus the total gate complexity of timeevolution under the assumption of constant density (i.e. N/Ω∈O(1)) is
$\begin{array}{cc}{C}_{\mathrm{TDS}}=\ue52e\ue8a0\left(\frac{{N}^{8/3}}{{\Omega}^{2/3}}\ue89et\ue89e\phantom{\rule{0.3em}{0.3ex}}\ue89e\mathrm{poly}\ue89e\phantom{\rule{0.3em}{0.3ex}}\ue89e\mathrm{log}\ue8a0\left(N,{\alpha}_{T},t,\u03f5\right)\right)=\stackrel{~}{\ue52e}\ue8a0\left({N}^{2}\ue89et\ue89e\phantom{\rule{0.3em}{0.3ex}}\ue89e\mathrm{log}\ue8a0\left(1/\u03f5\right)\right).& \left(61\right)\end{array}$
In contrast, the cost of simulation in the planewave dual basis by prior art Ryan Babbush, Nathan Wiebe, Jarrod McClean, James McClain, Hartmut Neven, and Garnet Kin Chan, “Low depth quantum simulation of electronic structure,” arXiv preprint arXiv:1706.00023, 2017, applies the ‘Qubitization’ technique (see Guang Hao Low and Isaac L Chuang, “Hamiltonian simulation by qubitization,” arXiv preprint arXiv:1610.06546, 2016) and has gate complexity that scales like Õ(N^{11/3}t), and also has a polynomial dependence on the nuclear charge n_{j}. Embodiments of the disclosed method outperform this, and notably depends only polylogarithmically on the nuclear charges within the problem.
VI. Application to Sparse Hamiltonian Simulation
In this section, a complexitytheoretic perspective of the improvements that are enabled by simulation with the truncated Dyson series and simulation in the interaction picture are submitted. One can do so by evaluating the query complexity for the simulation of sparse Hamiltonians H. Such Hamiltonians of dimension N are called dsparse if there are at most d=(polylog(N)) nonzero entries in every row, and the position and values of these entries may be efficiently output, in say a binary representation, by some classical circuit of size (polylog(N)). This abstract model is useful in quantum complexity theory as a natural generalization these classical circuits leads to unitary quantum oracles that can be queried to access the same information, but now in superposition. With this model, one can achieve inSection VIA timedependent simulation with a squareroot improvement with respect to the sparsity parameter, and gate complexity scaling with the average instead of worstcase rateofchange ∥H∥. By moving to the interaction picture in Section VI B, one can find more efficient timeindependent simulation algorithms for diagonallydominant Hamiltonians.
This model assumes that the Hamiltonian is input to the simulation routine through two oracles: O_{f }and O_{H}. O_{H }is straight forward; it provides the values of the matrix elements of the Hamiltonian given a time index t and indices x, y to for the row and column of H as follows
O_{H}t,x,y,0)=t,x,y,H_{xy}(t) (62)
O_{f }provides the locations of the nonzero matrix elements in any given row or column of H. Specifically, let ƒ(x,j) give the column index of the j^{th }nonzero matrix element in row x if it exists and an appropriately chosen zero element if it does not. In particular, let r_{t,j }be the list of column indices of these nonzero matrix elements in row J. One can then define, with a timeindex t,
O_{ƒ}t,x,j=t,x,ƒ_{t}(x,j). (63)
A. Simulation of Sparse TimeDependent Hamiltonians
Applying the truncated Dyson series simulation algorithm to sparse Hamiltonians requires us to synthesize HAMT in Eq. (7) from these oracles O_{H}, O_{F}. This is possible by a straightforward construction.
Lemma 4 (Synthesis of HAMT from sparse Hamiltonian oracles). Let an N×N timedependent dsparse Hamiltonian H(s) with maxnorm ∥H∥_{max}:=max_{s}∥H(s)∥_{max }be defined on the interval s ∈[0, t] to n_{p }bits of precision. Then
$\begin{array}{cc}(\u30080\ue89e{\ue85c}_{a}\ue89e\otimes {I}_{s})\ue89e\mathrm{HAM}\ue89e\text{}\ue89eT\ue8a0\left({\uf6030\u3009}_{a}\otimes {I}_{s}\right)=\sum _{t}\ue89e\phantom{\rule{0.3em}{0.3ex}}\ue89e\uf603t\u3009\ue89e{\u3008t\uf604}_{d}\otimes \frac{H\ue8a0\left(t\right)}{d\ue89e{\uf605H\uf606}_{m\ue89e\phantom{\rule{0.3em}{0.3ex}}\ue89e\mathrm{ax}}},& \left(64\right)\end{array}$
can be implemented with O(1) queries to O_{f }and O_{H}, and (poly(n_{p})+log (N)) primitive gates.
Proof. Let U_{col}, U_{row }be the following unitary transformations
$\begin{array}{cc}{U}_{\mathrm{col}}\ue89e{\uf603t\u3009}_{d}\ue89e{\uf603k\u3009}_{s}\ue89e{\uf6030\u3009}_{a}:={\uf603t\u3009}_{d}\ue89e\uf603{\chi}_{k}\ue8a0\left(t\right)\u3009=\frac{1}{\sqrt{d}}\ue89e\sum _{p\in {r}_{k}}\ue89e{\uf603t\u3009}_{d}\ue89e{\uf603k\u3009}_{s}\ue89e{\uf603p\u3009}_{{a}_{1}}\ue89e\left(\sqrt{\frac{{H}_{k,p}^{*}\ue8a0\left(t\right)}{{\uf605H\uf606}_{\mathrm{ma}\ue89e\phantom{\rule{0.3em}{0.3ex}}\ue89ex}}}\ue89e{\uf6030\u3009}_{{a}_{2}}+\sqrt{1\frac{\uf603{H}_{k,p}\ue8a0\left(t\right)\uf604}{{\uf605H\uf606}_{\mathrm{ma}\ue89e\phantom{\rule{0.3em}{0.3ex}}\ue89ex}}}\ue89e{\uf6031\u3009}_{{a}_{2}}\right),& \left(65\right)\\ {\u30080\uf604}_{a}\ue89e{\u3008j\uf604}_{s}\ue89e{\u3008t\uf604}_{d}\ue89e{U}_{\mathrm{row}}^{\u2020}:=\u3008{\stackrel{\_}{\chi}}_{j}\ue8a0\left(t\right)\uf604\ue89e{\u3008t\uf604}_{d}=\frac{1}{\sqrt{d}}\ue89e\sum _{q\in {r}_{j}}\ue89e\left(\sqrt{\frac{{\delta}_{j,q}\ue89e{H}_{j,q}\ue8a0\left(t\right)+\left(1{\delta}_{j,q}\right)\ue89e{H}_{j,q}\ue8a0\left(t\right)}{{\uf605H\uf606}_{\mathrm{ma}\ue89e\phantom{\rule{0.3em}{0.3ex}}\ue89ex}}}\ue89e{\u30080\uf604}_{{a}_{2}}+\sqrt{1\frac{\uf603{H}_{q,j}\ue8a0\left(t\right)\uf604}{{\uf605H\uf606}_{\mathrm{ma}\ue89e\phantom{\rule{0.3em}{0.3ex}}\ue89ex}}}\ue89e{\u30082\uf604}_{{a}_{2}}\right)\ue89e{\u3008j\uf604}_{{a}_{1}}\ue89e{\u3008q\uf604}_{s}\ue89e{\u3008t\uf604}_{d},& \left(66\right)\\ \phantom{\rule{4.4em}{4.4ex}}\ue89e\u3008{\stackrel{\_}{\chi}}_{j}\ue8a0\left(t\right){\chi}_{k}\ue8a0\left(t\right)\u3009=\frac{\sqrt{{H}_{j,k}\ue8a0\left(t\right)\ue89e{H}_{k,j}^{*}\ue8a0\left(t\right)}}{d\ue89e{\uf605H\uf606}_{\mathrm{ma}\ue89e\phantom{\rule{0.3em}{0.3ex}}\ue89ex}}=\frac{{H}_{j,k}}{d\ue89e{\uf605H\uf606}_{\mathrm{ma}\ue89e\phantom{\rule{0.3em}{0.3ex}}\ue89ex}}.& \left(67\right)\end{array}$
Let φ=Σ_{t,k}α_{t,k}t_{d}k)_{s}. One then has that
$\begin{array}{cc}\begin{array}{c}\left[\uf6030\u3009\ue89e{\u30080\uf604}_{a}\otimes \right]\ue89e{U}_{\mathrm{row}}^{\u2020}\xb7{U}_{\mathrm{col}}\ue89e{\uf603t\u3009}_{d}\ue89e\uf603\psi \u3009\ue89e{\uf6030\u3009}_{a}=\ue89e\uf6030\u3009\ue89e{\u30080\uf604}_{a}\otimes \sum _{{t}^{\prime},j}\ue89e\phantom{\rule{0.3em}{0.3ex}}\ue89e\left(\uf603{t}^{\prime}\u3009\ue89e{\u3008{t}^{\prime}\uf604}_{d}\otimes \uf603j\u3009\ue89e{\u3008j\uf604}_{s}\right)\\ \ue89e\sum _{t,k}\ue89e\phantom{\rule{0.3em}{0.3ex}}\ue89e{a}_{t,k}\ue89e{U}_{\mathrm{row}}^{\u2020}\xb7{U}_{\mathrm{col}}\ue89e{\uf603t\u3009}_{d}\ue89e{\uf603k\u3009}_{s}\ue89e{\uf6030\u3009}_{a}.\\ =\ue89e{\uf6030\u3009}_{a}\ue89e\sum _{{t}^{\prime},j}\ue89e{\uf603{t}^{\prime}\u3009}_{d}\ue89e{\uf603j\u3009}_{s}\ue89e\sum _{t,k}\ue89e{a}_{t,k}\ue8a0\left(\u3008{\chi}_{j}\ue8a0\left({t}^{\prime}\right)\ue89e\uf603{\u3008{t}^{\prime}\uf604}_{d})\ue89e({\uf603t\u3009}_{d}\uf604\ue89e{\chi}_{k}\ue8a0\left(t\right)\u3009\right)\\ =\ue89e\sum _{j,k}\ue89e{a}_{t,k}\ue89e{\uf6030\u3009}_{a}\ue89e{\uf603t\u3009}_{d}\ue89e{\uf603j\u3009}_{s}\ue89e\frac{{H}_{j,k}\ue8a0\left(t\right)}{d\ue89e{\uf605H\uf606}_{\mathrm{ma}\ue89e\phantom{\rule{0.3em}{0.3ex}}\ue89ex}}=\frac{{\uf6030\u3009}_{s}\ue89eH\ue8a0\left(t\right)\ue89e\uf603\psi \u3009}{d\ue89e{\uf605H\uf606}_{\mathrm{ma}\ue89e\phantom{\rule{0.3em}{0.3ex}}\ue89ex}}.\end{array}& \left(68\right)\end{array}$
As this result holds for any input state φ, the choice HAMT=U_{row}^{†}·U_{col }satisfies Eq. (64).
The query cost then follows from the fact that U_{col }and U_{row}^{†}can be implemented using O(1) calls to O_{ƒ}. In particular, U_{col }can be prepared in the following steps:
$\begin{array}{cc}\phantom{\rule{4.4em}{4.4ex}}\ue89e\uf603t\u3009\ue89e\uf603k\u3009\ue89e\uf6030\u3009\mapsto \uf603t\u3009\ue89e\uf603k\u3009\ue89e\frac{1}{\sqrt{d}}\ue89e\sum _{\ue54b=1}^{d}\ue89e\phantom{\rule{0.3em}{0.3ex}}\ue89e\uf603\ue54b\u3009\ue89e\uf6030\u3009\ue89e\text{}\ue89e\phantom{\rule{4.4em}{4.4ex}}\ue89e\ue89e\uf603t\u3009\ue89e\uf603k\u3009\ue89e\frac{1}{\sqrt{d}}\ue89e\sum _{p\in {r}_{k}}\ue89e\phantom{\rule{0.3em}{0.3ex}}\ue89e\uf603p\u3009\ue89e\uf6030\u3009\ue89e\text{}\ue89e\phantom{\rule{4.4em}{4.4ex}}\ue89e\ue89e\uf603t\u3009\ue89e\uf603k\u3009\ue89e\phantom{\rule{0.3em}{0.3ex}}\ue89e\frac{1}{\sqrt{d}}\ue89e\sum _{p\in {r}_{k}}\ue89e\uf603p\u3009\ue89e\uf603{H}_{k,p}\ue8a0\left(t\right)\u3009\ue89e\uf6030\u3009\ue89e\text{}\mapsto \uf603t\u3009\ue89e\uf603k\u3009\ue89e\phantom{\rule{0.3em}{0.3ex}}\ue89e\frac{1}{\sqrt{d}}\ue89e\sum _{p\in {r}_{k}}\ue89e\uf603p\u3009\ue89e\uf603{H}_{k,p}\ue8a0\left(t\right)\u3009\ue89e\phantom{\rule{0.3em}{0.3ex}}\ue89e\left(\sqrt{\frac{{H}_{k,p}^{*}\ue8a0\left(t\right)}{{\uf605H\uf606}_{\mathrm{ma}\ue89e\phantom{\rule{0.3em}{0.3ex}}\ue89ex}}}\ue89e\uf6030\u3009+\sqrt{1\frac{\uf603{H}_{k,p}\ue8a0\left(t\right)\uf604}{{\uf605H\uf606}_{\mathrm{ma}\ue89e\phantom{\rule{0.3em}{0.3ex}}\ue89ex}}}\ue89e\uf6031\u3009\right)\ue89e\text{}\ue89e\ue89e\uf603t\u3009\ue89e\uf603k\u3009\ue89e\phantom{\rule{0.6em}{0.6ex}}\ue89e\frac{1}{\sqrt{d}}\ue89e\sum _{p\in {r}_{k}}\ue89e\uf603p\u3009\ue89e\phantom{\rule{0.3em}{0.3ex}}\ue89e\left(\sqrt{\frac{{H}_{k,p}^{*}\ue8a0\left(t\right)}{{\uf605H\uf606}_{\mathrm{ma}\ue89e\phantom{\rule{0.3em}{0.3ex}}\ue89ex}}}\ue89e\uf6030\u3009+\sqrt{1\frac{\uf603{H}_{k,p}\ue8a0\left(t\right)\uf604}{{\uf605H\uf606}_{\mathrm{ma}\ue89e\phantom{\rule{0.3em}{0.3ex}}\ue89ex}}}\ue89e\uf6031\u3009\right)\ue89e\uf6030\u3009={U}_{\mathrm{col}}\ue89e\uf603t\u3009\ue89e\uf603k\u3009\ue89e\uf6030\u3009.& \left(69\right)\end{array}$
Therefore accessing U_{col }unitaries requires O(1) queries to the fundamental oracles as claimed, along with a arithmetic circuit, of size polynomial in the number of bits used to represent H_{k,p}(t), for computing trigonometric functions of the magnitudes of the complexvalued matrix elements as well as their arguments. The argument that U_{row}^{†}requires O(1) queries follows in exactly the same manner, but with an additional final step that swaps the s and α_{1 }registers.□
Once HAMT, is obtained, the complexity of simulation follows directly from previous results.
Theorem 4 (Simulation of sparse timedependent Hamiltonians). Let an N×N timedependent Hamiltonian H(s) for s ∈ [0, t] with maxnorm, ∥H∥_{max}:=max_{s}∥H(s)∥_{max }and average rateofchange
$\u3008\uf605\stackrel{.}{H}\uf606\u3009=\frac{1}{t}\ue89e{\int}_{0}^{t}\ue89e\uf605\frac{\mathrm{dH}\ue8a0\left(s\right)}{\mathrm{ds}}\uf606\ue89e\phantom{\rule{0.2em}{0.2ex}}$
ds be encoded in the oracles O_{H }and O_{ƒ }from Eqs. (62) and (63) to n_{p }bits of precision. Let α=d∥H∥_{max }and =tα. Then the timeordered evolution, operator [exp(−iƒ_{0}^{t}H(s)ds)] may be approximated with error ϵ using
$\begin{array}{cc}\phantom{\rule{4.4em}{4.4ex}}\ue89e\mathrm{Queries}\ue89e\phantom{\rule{0.8em}{0.8ex}}\ue89e\mathrm{to}\ue89e\phantom{\rule{0.8em}{0.8ex}}\ue89e{O}_{H}\ue89e\phantom{\rule{0.8em}{0.8ex}}\ue89e\mathrm{and}\ue89e\phantom{\rule{0.8em}{0.8ex}}\ue89e{O}_{f}\ue89e\text{:}\ue89e\phantom{\rule{0.8em}{0.8ex}}\ue89e\ue52e\ue89e\phantom{\rule{0.8em}{0.8ex}}\ue89e\left(\tau \ue89e\frac{\mathrm{log}\ue8a0\left(\tau /\u03f5\right)}{\mathrm{log}\ue89e\phantom{\rule{0.3em}{0.3ex}}\ue89e\mathrm{log}\ue8a0\left(\tau /\u03f5\right)}\right).& 1.\\ \phantom{\rule{4.4em}{4.4ex}}\ue89e\mathrm{Qubits}\ue89e\text{:}\ue89e\phantom{\rule{0.8em}{0.8ex}}\ue89e\ue52e\ue8a0\left({n}_{p}+\mathrm{log}\ue89e\phantom{\rule{0.8em}{0.8ex}}\ue89e\left(\frac{{N}_{\tau}}{\u03f5\ue89e\phantom{\rule{0.3em}{0.3ex}}\ue89e{\alpha}^{2}}\ue89e\left(\frac{\u3008\uf605\stackrel{.}{H}\uf606\u3009}{\alpha}+\frac{{\uf605H\uf606}^{2}}{{\alpha}^{2}}\right)\right)\right).& 2.\\ \mathrm{Primitive}\ue89e\phantom{\rule{0.8em}{0.8ex}}\ue89e\mathrm{gates}\ue89e\text{:}\ue89e\phantom{\rule{0.8em}{0.8ex}}\ue89e\ue52e\ue8a0\left(\mathrm{poly}\ue8a0\left({n}_{p}\right)+\mathrm{log}\ue8a0\left(\frac{{N}_{\tau}}{\u03f5\ue89e\phantom{\rule{0.3em}{0.3ex}}\ue89e{\alpha}^{2}}\ue89e\left(\frac{\u3008\uf605\stackrel{.}{H}\uf606\u3009}{\alpha}+\frac{{\uf605H\uf606}^{2}}{{\alpha}^{2}}\right)\right)\ue89e\frac{\mathrm{log}\ue8a0\left(\tau /\u03f5\right)}{\mathrm{log}\ue89e\phantom{\rule{0.3em}{0.3ex}}\ue89e\mathrm{log}\ue8a0\left(\tau /\u03f5\right)}\right).& 3.\end{array}$
Proof. From Theorem 1, the number of qubits requried is (n_{s}+n_{a}+n_{d}+log log (1/ϵ)). Values for these parameters are obtained from the construction of HAMT in Lemma 4, which also requires an additional n_{p }qubits for the bits of precision to which H is encoded. In this construction, n_{α}=(n_{s})=(log (N)). n_{d }is obtained from the number of timediscretization points required by Theorem 1. Simulation for time t is implemented by simulating segments of duration ⊖(α^{−1}). As there are ⊖(tα) segments, one can rescale ϵ→⊖(ϵ/(tα)).□
This is a quadratic improvement in sparsity d. Furthermore, instead of scaling with the worstcase (log (max_{s}∥{dot over (H)}(s)∥)), one can obtain scaling with average rateofchange ∥{dot over (H)}∥. If one further assumes that the computed matrix elements H_{j,k }are not exact, the number of bit of precision scales like n_{p}=(log (∥H∥t/ϵ)). Note several generic improvements to Theorem 4 are possible, but will not be pursued further as they are straightforward. For instance, if ∥H(t)∥_{max }as a function of time is known, one may use step sizes of varying size by encoding each segment t ∈[t_{j},t_{j+1}] with the largest max_{t∈[t}_{j}_{,t}_{j+1}_{]}∥H(t)∥_{max}, rather than the worstcase max_{t}∥H(t)∥_{max}.
B. Simulation of Sparse TimeIndependent Hamiltonians in the Interaction Picture
This section concerns timeindependent dsparse Hamiltonians H=A +B where A is diagonal and B_{kk}=0 for all k. In particular, consider the case of diagonally dominant Hamiltonians, where ∥A∥≥d∥B∥_{max}. Given norms for each of these terms ∥A∥ and ∥B∥_{max}, it is straightforward to simulate timeevolution e^{−iHt }in the Schrödinger picture. For instance, using the truncated Taylor series approach in Eq. (38), one obtains a query complexity of (t(d∥B∥_{max}+∥A∥) polylog(t,d,∥A∥,∥B∥_{max}, ϵ)). By instead simulating H_{I}(t)=e^{iAt}Be^{−iAt }in the interaction picture, the dependence on ∥A∥ can be removed, which is particularly useful in cases of strong diagonal dominance ∥A∥≥d∥B∥_{max}, of which the Hubbard model with longranged interactions in Section V is an example. Similar to the disclosed results for timedependent sparse Hamiltonian simulation in Section VIA, this is easily proven by mapping the input oracles O_{H }and O_{ƒ}for matrix values and positions to the oracles of Theorem 3 for the more general result.
Theorem 5 (Simulation of sparse diagonally dominant Hamiltonians). Let an N×N timeindependent Hamiltonian H for s ∈[0,t] with maxnorm ∥A∥ for the diagonal component and maxnorm ∥B∥_{max }for the offdiagonal component be encoded in the oracles O_{H }and O_{ƒ }from Eqs. (62) and (63) to n_{p }bits of precision. Let α_{b}=d∥B∥_{max}. Then the timeevolution operator e^{−iHt }may be approximated with error ϵ using
$\begin{array}{cc}\phantom{\rule{4.4em}{4.4ex}}\ue89e\mathrm{Queries}\ue89e\phantom{\rule{0.8em}{0.8ex}}\ue89e\mathrm{to}\ue89e\phantom{\rule{0.8em}{0.8ex}}\ue89e{O}_{H}\ue89e\phantom{\rule{0.8em}{0.8ex}}\ue89e\mathrm{and}\ue89e\phantom{\rule{0.8em}{0.8ex}}\ue89e{O}_{f}\ue89e\text{:}\ue89e\phantom{\rule{0.8em}{0.8ex}}\ue89e\ue52e\ue89e\phantom{\rule{0.8em}{0.8ex}}\ue89e\left({\alpha}_{B}\ue89et\ue89e\frac{\mathrm{log}\ue8a0\left({\alpha}_{B}\ue89e\tau /\u03f5\right)}{\mathrm{log}\ue89e\phantom{\rule{0.3em}{0.3ex}}\ue89e\mathrm{log}\ue8a0\left({\alpha}_{B}\ue89e\tau /\u03f5\right)}\right).& 1.\\ \phantom{\rule{4.4em}{4.4ex}}\ue89e\mathrm{Qubits}\ue89e\text{:}\ue89e\phantom{\rule{0.8em}{0.8ex}}\ue89e\ue52e\ue8a0\left({n}_{p}+\mathrm{log}\ue8a0\left(N\right)+\mathrm{log}\ue89e\phantom{\rule{0.8em}{0.8ex}}\ue89e\left(\frac{t}{\u03f5}\ue89e\frac{\uf605A\uf606}{{\alpha}_{B}}\right)\right).& 2.\\ \mathrm{Primitive}\ue89e\phantom{\rule{0.8em}{0.8ex}}\ue89e\mathrm{gates}\ue89e\text{:}\ue89e\phantom{\rule{0.8em}{0.8ex}}\ue89e\ue52e\ue8a0\left(\left(\mathrm{log}\ue8a0\left(N\right)+\mathrm{poly}\ue8a0\left({n}_{p}\right)\ue89e\mathrm{log}\ue8a0\left(\frac{t}{\u03f5}\ue89e\frac{\uf605A\uf606}{{\alpha}_{B}}\right)\right)\ue89e{\alpha}_{B}\ue89et\ue89e\frac{\mathrm{log}\ue8a0\left({\alpha}_{B}\ue89et/\u03f5\right)}{\mathrm{log}\ue89e\phantom{\rule{0.3em}{0.3ex}}\ue89e\mathrm{log}\ue8a0\left({\alpha}_{B}\ue89et/\u03f5\right)}\right)\ue89e\phantom{\rule{0.8em}{0.8ex}}\ue89e\mathrm{primitive}\ue89e\phantom{\rule{0.8em}{0.8ex}}\ue89e\mathrm{gates}.& 3.\end{array}$
Proof. This follows immediately by combining the query complexity of Lemma 3 to e^{−iAt }and HAMT that encodes the Hamiltonian H_{I}(t)=e^{iAt}Be^{−iAt }in the rotating frame, with the query complexity of the approach in Theorem 3 for synthesizing these oracles using the input oracles O_{H }and O_{ƒ}. One possible decomposition of HAMT is
$\begin{array}{cc}\mathrm{HAM}\ue89e\text{}\ue89eT=\left(\sum _{m=0}^{M1}\ue89e\phantom{\rule{0.3em}{0.3ex}}\ue89e\uf603m\u3009\ue89e{\u3008m\uf604}_{d}\otimes {I}_{a}\otimes {e}^{\mathrm{iArm}/M}\right)\ue89e\left({I}_{d}\otimes {O}_{B}\right)\ue89e\left(\sum _{m=0}^{M1}\ue89e\uf603m\u3009\ue89e{\u3008m\uf604}_{d}\otimes {I}_{a}\otimes {e}^{\mathrm{iArm}/M}\right),\text{}\ue89e\phantom{\rule{4.4em}{4.4ex}}\ue89e\left({\u30080\uf604}_{a}\otimes {I}_{s}\right)\ue89e{O}_{B}\ue8a0\left({\u30080\uf604}_{a}\otimes {I}_{s}\right)=\frac{B}{{\alpha}_{B}},& \left(70\right)\end{array}$
where α_{B}=d∥B∥_{max}, =(α_{B}^{−1}) and
$M=\ue52e\ue8a0\left(\frac{t}{\u03f5}\ue89e\frac{\uf605A\uf606}{{\alpha}_{B}}\right).$
Note that with this construction,
$\left({\u30080\uf604}_{a}\otimes {I}_{s}\right)\ue89e\mathrm{HAM}\ue89e\text{}\ue89eT\ue8a0\left({\u30080\uf604}_{a}\otimes {I}_{s}\right)={\Sigma}_{m=0}^{M1}\ue89e\phantom{\rule{0.3em}{0.3ex}}\ue89e\uf603m\u3009\ue89e{\u3008m\uf604}_{d}\otimes \frac{{H}_{I}\ue8a0\left(\tau \ue89e\phantom{\rule{0.3em}{0.3ex}}\ue89em/M\right)}{d\ue89e{\uf605B\uf606}_{m\ue89e\phantom{\rule{0.3em}{0.3ex}}\ue89e\mathrm{ax}}}.$
First, synthesize (Σ_{m=0}^{M−1}∥mm_{d}⊗e^{iAτm/M}) using (1) queries to the input oracles O_{H}, and (log(N)+n_{p}log(M)) primitive gates. Since A is diagonal, e^{−iAt }can be simulated for any t>0 using only two queries. This is implemented by the following steps:
$\begin{array}{cc}\uf603k\u3009\ue89e\uf6030\u3009\mapsto \uf603k\u3009\ue89e\uf603k\u3009\ue89e\uf6030\u3009\ue89e\underset{{O}_{H}}{\mapsto}\ue89e\uf603k\u3009\ue89e\uf603k\u3009\ue89e\uf603{H}_{k,k}\u3009\ue89e\uf6030\u3009\mapsto \uf603k\u3009\ue89e\uf603k\u3009\ue89e{e}^{{\mathrm{iH}}_{k,k}\ue89e\mathrm{Zt}}\ue89e\uf6030\u3009={e}^{{\mathrm{iH}}_{k,k}\ue89et}\ue89e\uf603k\u3009\ue89e\uf603k\u3009\ue89e\uf6030\u3009\ue89e\underset{{O}_{H}^{\u2020}}{\mapsto}\ue89e{e}^{{\mathrm{iH}}_{k,k}\ue89et}\ue89e\uf603k\u3009\ue89e\uf603k\u3009\ue89e\uf6030\u3009\to {e}^{{\mathrm{iH}}_{k,k}\ue89et}\ue89e\uf603k\u3009\ue89e\uf6030\u3009.& \left(71\right)\end{array}$
Step one uses n_{s}=(log (N)) CNOT gates to copy the computational basis state k. Step three applies (n_{p}) phase rotation with angle controlled by the bits of H_{k,k}, and the value of t, which is given beforehand. Subsequently, (Σ_{m=0}^{M−1}mm_{d}⊗I_{α}⊗e^{iAtm/M}) may be implemented by a sequence of rotations with angles increasing in a geometric series, and each controlled by a different qubit in the d register, e.g. controlled−e^{iAτ2}^{−1}^{/M},e^{iAτ2}^{−2}^{/M},e^{iAτ2}^{−n+p}^{/M}. Naively, this requires (log (M)) queries. However, it is only necessary to compute H_{k,k} once as the entire sequence of controlledphases may be applied after step three. Similarly, it is only necessary to copy the computational basis state (1) times.
Second, synthesize O_{B }using (1) queries to the input oracles O_{H }and O_{ƒ}. How this is done should be clear from Lemma 4, by omitting the timeindex, and preparing the state 00_{α}_{2}+11_{α}_{2 }in Eq. (65) when the input indices k=p. This has gate complexity (poly(n_{p})+log (N)). Thus e^{−iAτ }and HAMT combined have query complexity (1) to O_{H }and O_{ƒ}, and gate complexity (poly(n_{p}) log (M)+log (N)). By substituting into Lemma 3, one can obtain the stated results.□
This provides a formal proof that the query complexity of simulating a Hamiltonian, within the interaction picture, is independent of the magnitude of the diagonal elements of the Hamiltonian.
VII. General Embodiments
In this section, example methods for performing aspects of the disclosed embodiments are disclosed. The particular embodiments described should not be construed as limiting, as the disclosed method acts can be performed alone, in different orders, or at least partially simultaneously with one another. Further, any of the disclosed methods or method acts can be performed with any other methods or method acts disclosed herein.
FIG. 6 is an example method for a timedependent simulation algorithm as disclosed herein.
In FIG. 6, the inputs are shown at 610. The inputs include quantum state subroutine (S); for computing (H(t)), evolution time (t), number of segments (r), order (K). At 612, a determination is made as to whether any segments remain. If so, the process proceeds to 614, where a robust oblivious amplitude amplication procedure is applied. If the procedure of 614 is done, then the process proceeds to the next segment; if the procedure 614 is not done, then the procedure proceeds to 616. At 616, a weighted quantum superposition is prepared over the first K terms in the Dyson series of time evolution operator. At 618, a quantum superposition is prepared over each of the times that the Hamiltonian is evaluated at. At 620, subroutine (S) is applied to construct a state that stores the columns of a rescaling of the kth term in the dyson series, for each k=0 . . . K. At 622, subroutine (S) is applied to construct a state that stores the rows of a rescaling of the kth term in the dyson series, for each k=0 . . . , K. At 624, linear combinations of unitary methods are used to conditionally implement the evolution to the input state of 614.
FIG. 7 is an example method for an example interaction picture simulation method as disclosed herein.
At 710, a quantum state, Hamiltonian function, and evolution time are input. At 712, a classical computer is used to construct an interaction picture Hamiltonian by conjugating each term in the offdiagonal Hamiltonian by the evolution according to the diagonal terms. At 714, the interaction picture version of the Hamiltonian is simulated on a quantum computer using a timedependent Hamiltonian simulation method (e.g., such as that of FIG. 1). At 716, the interaction picture is transformed back on the quantum computer by counterrotating the quantum state. At 718, the resulting quantum state is output.
FIG. 8 is an example method simulating chemistry and Hubbard models as disclosed herein. At 810, a quantum state, Hamiltonian function, and an evolution time are input. At 182, a classical computer is used to transform to a frame where a potential operator is diagonal. At 814, a simulation method (e.g., the method of FIG. 2) is applied to a quantum state by the classical computer. And, at 816, an evolved quantum state is output (e.g., by the classical computer).
FIG. 9 is an example method of a compression method as disclosed herein. At 910, numerous values are input, including one or more of (a) a classical variable j=1; (b) a success quantum register A; (c) target quantum register S; (d) counter quantum register B; (e) counter quantum register C; (f) a set of J quantum circuits Q_{1},Q_{2},Q_{j}, where Q_{j }applies matrix M_{j }on register S conditional register A being in a state that indicates ‘success. At 912, it is determined whether j=j+1. If not, then the method proceeds to 914. At 914, and conditional on B and C being the initial quantum state, Q_{j }is applied. Then, at 916, j is incremented by 1, and conditional on the A being in a ‘success quantum state, B is decremented and C is incremeneted with quantum arithmetic. 916 leads back into 912, until it is determined that j=j+1. If so, at 918, the method outputs a quantum circuit that applies on register S the product of matrices M_{1 }. . . M_{k }conditional on register A and C being in the initial quantum state, and register b being in the quantum state k.
FIG. 14 illustrates an example method for performing a quantum simulation in accordance with an embodiment of the disclosed technology.
At 1410, a quantum computer is configured to simulate a quantum system, wherein a Hamiltonian in the simulation is represented in the interaction picture. At 1412, a simulation of the quantum system is performed using the quantum computer.
In certain implementations, the simulation of the quantum system is a subroutine that is repeated two or more times. In some implementations, the simulation is performed using linear combinations of matrices. For instance, in some cases, the simulation uses linear combinations of unitaries performed on a diagonally dominant matrix (e.g., using linear combinations of unitaries performed on the diagonally dominant components of the diagonally dominant matrix). In certain implementations, the quantum system is modelled by a Hubbard model. In particular implementations, the quantum system describes a physical chemical system or molecule. In some implementations, the Hamiltonian is sparse and the simulation uses a state of an auxillary qubit to encode matrix elements of the Hamiltonian instead of using graph decomposition techniques. In certain implementations, the simulation is performed by compressing ancillas for quantum simulation of a timedependent Hamiltonian.
FIG. 15 illustrates a further example method for performing a quantum simulation in accordance with an embodiment of the disclosed technology.
At 1510, a quantum algorithm is implemented on a quantum computer for simulating a general sparse timedependent quantum system. In this embodiment, the quantum algorithm does not use graph decomposition techniques.
In certain implementations, a Hamiltonian used in the simulation is represented in the interaction picture. In some implementations, the simulation uses linear combinations of unitaries performed on a diagonally dominant matrix. In particular implementations, a Hamilltonian used in the simulation is represented in the interaction picture and is chosen to be a diagonal matrix. In certain implementations the quantum system is modelled by a Hubbard model, a physical chemical system, or a molecule. In further implementations, the simulation includes compressing ancillas used to index a time for the simulation.
VIII. Example Computing Environments
FIG. 10 illustrates a generalized example of a suitable classical computing environment 1000 in which aspects of the described embodiments can be implemented. The computing environment 1000 is not intended to suggest any limitation as to the scope of use or functionality of the disclosed technology, as the techniques and tools described herein can be implemented in diverse generalpurpose or specialpurpose environments that have computing hardware.
With reference to FIG. 10, the computing environment 1000 includes at least one processing device 1010 and memory 1020. In FIG. 2, this most basic configuration 1030 is included within a dashed line. The processing device 1010 (e.g., a CPU or microprocessor) executes computerexecutable instructions. In a multiprocessing system, multiple processing devices execute computerexecutable instructions to increase processing power. The memory 1020 may be volatile memory (e.g., registers, cache, RAM, DRAM, SRAM), nonvolatile memory (e.g., ROM, EEPROM, flash memory), or some combination of the two. The memory 1020 stores software 1080 implementing tools for peforming any of the disclosed simulation techniques in the quantum computer as described herein. The memory 1020 can also store software 1080 for synthesizing, generating, or compiling quantum circuits for performing the described simulation techniques using quantum computing devices as described herein.
The computing environment can have additional features. For example, the computing environment 1000 includes storage 1040, one or more input devices 1050, one or more output devices 1060, and one or more communication connections 1070. An interconnection mechanism (not shown), such as a bus, controller, or network, interconnects the components of the computing environment 1000. Typically, operating system software (not shown) provides an operating environment for other software executing in the computing environment 1000, and coordinates activities of the components of the computing environment 1000.
The storage 1040 can be removable or nonremovable, and includes one or more magnetic disks (e.g., hard drives), solid state drives (e.g. flash drives), magnetic tapes or cassettes, CDROMs, DVDs, or any other tangible nonvolatile storage medium which can be used to store information and which can be accessed within the computing environment 1000. The storage 1040 can also store instructions for the software 1080 implementing any of the disclosed simulation techniques in a quantum computing device. The storage 1040 can also store instructions for the software 1080 for generating and/or synthesizing any of the described techniques, systems, or quantum circuits.
The input device(s) 1050 can be a touch input device such as a keyboard, touchscreen, mouse, pen, trackball, a voice input device, a scanning device, or another device that provides input to the computing environment 1000. The output device(s) 1060 can be a display device (e.g., a computer monitor, laptop display, smartphone display, tablet display, netbook display, or touchscreen), printer, speaker, or another device that provides output from the computing environment 1000.
The communication connection(s) 1070 enable communication over a communication medium to another computing entity. The communication medium conveys information such as computerexecutable instructions or other data in a modulated data signal. A modulated data signal is a signal that has one or more of its characteristics set or changed in such a manner as to encode information in the signal. By way of example, and not limitation, communication media include wired or wireless techniques implemented with an electrical, optical, RF, infrared, acoustic, or other carrier.
As noted, the various methods, techniques for controlling a quantum computing device, circuit design techniques, or compilation/synthesis techniques can be described in the general context of computerreadable instructions stored on one or more computerreadable media. Computerreadable media are any available media (e.g., memory or storage device) that can be accessed within or by a computing environment. Computerreadable media include tangible computerreadable memory or storage devices, such as memory 1020 and/or storage 1040, and do not include propagating carrier waves or signals per se (tangible computerreadable memory or storage devices do not include propagating carrier waves or signals per se).
Various embodiments of the methods disclosed herein can also be described in the general context of computerexecutable instructions (such as those included in program modules) being executed in a computing environment by a processor. Generally, program modules include routines, programs, libraries, objects, classes, components, data structures, and so on, that perform particular tasks or implement particular abstract data types. The functionality of the program modules may be combined or split between program modules as desired in various embodiments. Computerexecutable instructions for program modules may be executed within a local or distributed computing environment.
An example of a possible network topology 1100 (e.g., a clientserver network) for implementing a system according to the disclosed technology is depicted in FIG. 11. Networked computing device 1120 can be, for example, a computer running a browser or other software connected to a network 1112. The computing device 1120 can have a computer architecture as shown in FIG. 10 and discussed above. The computing device 1120 is not limited to a traditional personal computer but can comprise other computing hardware configured to connect to and communicate with a network 1112 (e.g., smart phones, laptop computers, tablet computers, or other mobile computing devices, servers, network devices, dedicated devices, and the like). Further, the computing device 1120 can comprise an FPGA or other programmable logic device. In the illustrated embodiment, the computing device 1120 is configured to communicate with a computing device 1130 (e.g., a remote server, such as a server in a cloud computing environment) via a network 1112. In the illustrated embodiment, the computing device 1120 is configured to transmit input data to the computing device 1130, and the computing device 1130 is configured to implement a technique for controlling a quantum computing device according to any of the disclosed embodiments and/or a circuit generation/compilation/synthesis technique for generating quantum circuits for performing any of the techniques disclosed herein. The computing device 1130 can output results to the computing device 1120. Any of the data received from the computing device 1130 can be stored or displayed on the computing device 1120 (e.g., displayed as data on a graphical user interface or web page at the computing devices 1120). In the illustrated embodiment, the illustrated network 1112 can be implemented as a Local Area Network (LAN) using wired networking (e.g., the Ethernet IEEE standard 802.3 or other appropriate standard) or wireless networking (e.g. one of the IEEE standards 802.11a, 802.11b, 802.11g, or 802.11n or other appropriate standard). Alternatively, at least part of the network 1112 can be the Internet or a similar public network and operate using an appropriate protocol (e.g., the HTTP protocol).
Another example of a possible network topology 1200 (e.g., a distributed computing environment) for implementing a system according to the disclosed technology is depicted in FIG. 12. Networked computing device 1220 can be, for example, a computer running a browser or other software connected to a network 1212. The computing device 1220 can have a computer architecture as shown in FIG. 12 and discussed above. In the illustrated embodiment, the computing device 1220 is configured to communicate with multiple computing devices 1230, 1231, 1232 (e.g., remote servers or other distributed computing devices, such as one or more servers in a cloud computing environment) via the network 1212. In the illustrated embodiment, each of the computing devices 1230, 1231, 1232 in the computing environment 1200 is used to perform at least a portion of a technique for controlling a quantum computing device according to any of the disclosed embodiments and/or a circuit generation/compilation/synthesis technique for generating quantum circuits for performing any of the techniques disclosed herein. In other words, the computing devices 1230, 1231, 1232 form a distributed computing environment in which aspects of the techniques as disclosed herein and/or quantum circuit generation/compilation/synthesis processes are shared across multiple computing devices. The computing device 1220 is configured to transmit input data to the computing devices 1230, 1231, 1232, which are configured to distributively implement such as process, including performance of any of the disclosed methods or creation of any of the disclosed circuits, and to provide results to the computing device 1220. Any of the data received from the computing devices 1230, 1231, 1232 can be stored or displayed on the computing device 1220 (e.g., displayed as data on a graphical user interface or web page at the computing devices 1220). The illustrated network 1212 can be any of the networks discussed above with respect to FIG. 11.
With reference to FIG. 13, an exemplary system for implementing the disclosed technology includes computing environment 1300. In computing environment 1300, a compiled quantum computer circuit description (including quantum circuits configured to perform any of the disclosed simulation techniques as disclosed herein) can be used to program (or configure) one or more quantum processing units such that the quantum processing unit(s) implement the circuit described by the quantum computer circuit description.
The environment 1300 includes one or more quantum processing units 1302 and one or more readout device(s) 1308. The quantum processing unit(s) execute quantum circuits that are precompiled and described by the quantum computer circuit description. The quantum processing unit(s) can be one or more of, but not limited to: (a) a superconducting quantum computer; (b) an ion trap quantum computer; (c) a faulttolerant architecture for quantum computing; and/or (d) a topological quantum architecture (e.g., a topological quantum computing device using Majorana zero modes). The precompiled quantum circuits, including any of the disclosed circuits, can be sent into (or otherwise applied to) the quantum processing unit(s) via control lines 1306 at the control of quantum processor controller 520. The quantum processor controller (QP controller) 1320 can operate in conjunction with a classical processor 1310 (e.g., having an architecture as described above with respect to FIG. 10) to implement the desired quantum computing process. In the illustrated example, the QP controller 1320 further implements the desired quantum computing process via one or more QP subcontrollers 1304 that are specially adapted to control a corresponding one of the quantum processor(s) 1302. For instance, in one example, the quantum controller 1320 facilitates implementation of the compiled quantum circuit by sending instructions to one or more memories (e.g., lowertemperature memories), which then pass the instructions to lowtemperature control unit(s) (e.g., QP subcontroller(s) 1304) that transmit, for instance, pulse sequences representing the gates to the quantum processing unit(s) 1302 for implementation. In other examples, the QP controller(s) 1320 and QP subcontroller(s) 1304 operate to provide appropriate magnetic fields, encoded operations, or other such control signals to the quantum processor(s) to implement the operations of the compiled quantum computer circuit description. The quantum controller(s) can further interact with readout devices 1308 to help control and implement the desired quantum computing process (e.g., by reading or measuring out data results from the quantum processing units once available, etc.)
With reference to FIG. 13, compilation is the process of translating a highlevel description of a quantum algorithm into a quantum computer circuit description comprising a sequence of quantum operations or gates, which can include the circuits as disclosed herein (e.g., the circuits configured to perform one or more simulations as disclosed herein). The compilation can be performed by a compiler 1322 using a classical processor 1310 (e.g., as shown in FIG. 10) of the environment 1300 which loads the highlevel description from memory or storage devices 1312 and stores the resulting quantum computer circuit description in the memory or storage devices 1312.
In other embodiments, compilation and/or verification can be performed remotely by a remote computer 1360 (e.g., a computer having a computing environment as described above with respect to FIG. 10) which stores the resulting quantum computer circuit description in one or more memory or storage devices 1362 and transmits the quantum computer circuit description to the computing environment 1300 for implementation in the quantum processing unit(s) 1302. Still further, the remote computer 1300 can store the highlevel description in the memory or storage devices 1362 and transmit the highlevel description to the computing environment 1300 for compilation and use with the quantum processor(s)). In any of these scenarios, results from the computation performed by the quantum processor(s) can be communicated to the remote computer after and/or during the computation process. Still further, the remote computer can communicate with the QP controller(s) 1320 such that the quantum computing process (including any compilation, verification, and QP control procedures) can be remotely controlled by the remote computer 1360. In general, the remote computer 1360 communicates with the QP controller(s) 1320, compiler/synthesizer 1322, and/or verification tool 1323 via communication connections 1350.
In particular embodiments, the environment 1300 can be a cloud computing environment, which provides the quantum processing resources of the environment 1300 to one or more remote computers (such as remote computer 1360) over a suitable network (which can include the internet).
IX. Summation
The dynamics of interest in many condensed matter systems occur in the lowenergy subspace of the Hamiltonian. This work demonstrates that by simulating quantum dynamics in the interaction picture, the cost of simulation on a quantum computer for a large class of Hamiltonians can be made to scale with the lowenergy component. This is especially useful when the dynamics at high energy scales are, in a certain sense, simple, but nevertheless are necessary as they couple to low energy scales. This represents a significant advance compared to stateofart timeindependent quantum simulation algorithms that generally scales with the spectral norm of the full Hamiltonian. Finding further improvements in simulation algorithms specialized for lowenergy subspaces remains a major open problem. Embodiments of the disclosed approach are complementary to alternatives based on spectral gap amplification, and the possibility of combining these methods is an interesting future direction.
More generally, the complexity of timeindependent quantum simulation for generic Hamiltonians, given minimal information, appears to have been resolved. Thus future advancements, such as this work, will likely focus on exploiting the detailed structure of Hamiltonians of interest. For example, the result for timedependent sparse Hamiltonian simulation has scaling (d∥H∥_{max}t), but improvement to (√{square root over (d∥H∥_{max}∥H∥_{1}t)}) is straightforward, if the induced onenorm ∥H∥_{1 }of the Hamiltonian is also known beforehand. The promise of results in similar directions is exemplified by recent work that exploit the geometric locality of interactions, or the sizes of different terms in a Hamiltonian. The challenge will be finding characterizations of Hamiltonians that are sufficiently specific so as to enable a speedup, yet sufficiently general so as to include problems of practical and scientific value.
X. Appendix
A. The Truncated Taylor Series Algorithm
The truncated Taylor series simulation algorithm was a major advance in quantum simulation for its conceptual simplicity and computational efficiency. The original algorithm is motivated by truncating the Taylor expansion of the timeevolution operator at degree K.
$\begin{array}{cc}{e}^{\mathrm{iHt}}=1\mathrm{iHt}+\frac{{\left(\mathrm{iHt}\right)}^{2}}{2!}+\frac{{\left(\mathrm{iHt}\right)}^{3}}{3!}\ue89e\phantom{\rule{0.8em}{0.8ex}}\ue89e\dots =\underset{{\stackrel{\_}{R}\ue89e\phantom{\rule{0.3em}{0.3ex}}}_{K}}{\underset{\uf613}{\sum _{k=0}^{K}\ue89e\phantom{\rule{0.3em}{0.3ex}}\ue89e\frac{{\left(\mathrm{iHt}\right)}^{k}}{k!}}}+\underset{{R}_{K}}{\underset{\uf613}{\sum _{k=K+1}^{\infty}\ue89e\frac{{\left(\mathrm{iHt}\right)}^{k}}{k!}}}.& \left(72\right)\end{array}$
Assuming that t>0 and that the truncation order K≥2∥H∥t, the norms of R_{K }and the remainder term R_{K }are bounded by
$\begin{array}{cc}\phantom{\rule{4.4em}{4.4ex}}\ue89e\uf605{\stackrel{\_}{R}}_{K}\uf606=\uf605{e}^{\mathrm{iHt}}{R}_{K}\uf606\le 1+\uf605{R}_{K}\uf606,\text{}\ue89e\uf605{R}_{K}\uf606\le \sum _{k=K+1}^{\infty}\ue89e\phantom{\rule{0.3em}{0.3ex}}\ue89e\frac{{\left(\uf605H\uf606\ue89et\right)}^{k}}{k!}\le \frac{{\left(\uf605H\uf606\ue89et\right)}^{K+1}}{\left(K+1\right)!}\ue89e\sum _{k=K+2}^{\infty}\ue89e{\left(1/2\right)}^{kK1}=\frac{2\ue89e{\left(\uf605H\uf606\ue89et\right)}^{K+1}}{\left(K+1\right)!}.& \left(73\right)\end{array}$
Thus any unitary quantum circuit TTS that acts jointly on registers α,b,s and applies the nonunitary operator (00_{ab}⊗I_{s})TTS(00_{ab}⊗I_{s})≈R_{K }approximates the timeevolution operator with error δ and failure probability p given by
$\begin{array}{cc}\phantom{\rule{4.4em}{4.4ex}}\ue89e\delta =\uf605{e}^{\mathrm{iHt}}{\stackrel{\_}{R}}_{K}\uf606=\uf605{R}_{K}\uf606\le \frac{2\ue89e{\left(\uf605H\uf606\ue89et\right)}^{K+1}}{\left(K+1\right)!},\text{}\ue89ep\le 1\underset{{\uf603\psi \u3009}_{s}}{\mathrm{min}}\ue89e{\uf603\frac{{\stackrel{\_}{R}}_{K}\ue89e{\uf603\psi \u3009}_{s}}{1+\uf605{R}_{K}\uf606}\uf604}^{2}=1\underset{{\uf603\psi \u3009}_{s}}{\mathrm{min}}\ue89e{\uf603\frac{\left({e}^{\mathrm{iHt}}{R}_{K}\right)\ue89e{\uf603\psi \u3009}_{s}}{1+\uf605{R}_{K}\uf606}\uf604}^{2}\le 1{\uf603\frac{1\uf605{R}_{K}\uf606}{1+\uf605{R}_{K}\uf606}\uf604}^{2}=4\ue89e\uf605{R}_{K}\uf606=4\ue89e\delta .& \left(74\right)\end{array}$
Solving Eq. (74) for ∥H∥t=(1) gives the required truncation order
$K=\ue52e\ue8a0\left(\frac{\mathrm{log}\ue8a0\left(1/\delta \right)}{\mathrm{log}\ue89e\phantom{\rule{0.3em}{0.3ex}}\ue89e\mathrm{log}\ue8a0\left(1/\delta \right)}\right).$
The simulation algorithm TTS in FIG. 4 is obtained by constructing two oracles. HAM_{K}, which applies positive integer powers of (−iH)^{k }up to k=K, and COEF_{K}, which prepares a quantum state that selects these terms with the right coefficients. HAM_{K }will require additional ancilla registers, which can be indexed with {right arrow over (α)}and {right arrow over (b)}. Note that the gate and space complexity in the truncated Taylor series algorithm is dominated by that of HAM_{K}.
$\begin{array}{cc}\left({\u30080\uf604}_{\stackrel{>}{a}}\otimes {I}_{s\ue89e\stackrel{>}{b}}\right)\ue89e{\mathrm{HAM}}_{K}\ue8a0\left({\u30080\uf604}_{\stackrel{>}{a}}\otimes {I}_{s\ue89e\stackrel{>}{b}}\right):=\sum _{k=0}^{K}\ue89e\phantom{\rule{0.3em}{0.3ex}}\ue89e\uf603k\u3009\ue89e{\u3008k\uf604}_{\stackrel{>}{b}}\otimes {\left(\mathrm{iH}\right)}^{k},\text{}\ue89e{\mathrm{COEF}}_{K}\ue89e{\uf6030\u3009}_{\stackrel{>}{b}}:=\frac{1}{\sqrt{\beta}}\ue89e\sum _{k=0}^{K}\ue89e\sqrt{\frac{{t}^{k}}{k!}}\ue89e{\uf603k\u3009}_{\stackrel{>}{b}},\text{}\ue89e\beta =\sum _{k=0}^{K}\ue89e\frac{{t}^{k}}{k!}\le {e}^{t}.& \left(75\right)\end{array}$
FIG. 4 depicts the quantum circuit representation of (top, left) an example implementation of HAM_{K }from Eq. (75) using K queries to controlledHAM, (top, right) a single step of the truncated Taylor series algorithm before oblivious amplitude amplification, and (bottom) a single step of timeevolution by the truncated Taylor series algorithm from Eq. (78). Note that β=2 as a singleround of oblivious amplitude amplification is used. Thin horizontal lines depict singlequbit registers. Filled circles depict a unitary controlled on the 0 state.
The original algorithm implements HAM_{K }using K queries to controlledHAM
CHAM:=11_{b}⊗I_{as}+00_{b}⊗HAM (76)
with K copies of registers α and b. The state k_{{right arrow over (b)}}=0^{⊗k}1^{⊗K−k }that selects desired powers of H is encoded in unary, and so COEF_{K }may be implemented using (K) primitive gates. Up to a proportionality factor β, the unitaries of Eq. (75) allow us to implement the desired linear combination R_{K }for simulating timeevolution.
$\begin{array}{cc}{\mathrm{TTS}}_{\beta}:=\left({\mathrm{COEF}}_{K}^{\u2020}\otimes {I}_{\overrightarrow{a}\ue89es}\right)\ue89e{\mathrm{HAM}}_{K}\ue8a0\left({\mathrm{COEF}}_{K}\otimes {I}_{\overrightarrow{a}\ue89es}\right)\ue89e\text{}\ue89e(\u30080\ue89e{\ue85c}_{\overrightarrow{a}\ue89eb}\ue89e\otimes {I}_{s})\ue89e{{\mathrm{TTS}}_{\beta}(\ue85c0\u3009}_{\overrightarrow{a}\ue89eb}\otimes {I}_{s}=\frac{{\stackrel{\_}{R}}_{K}}{\beta}\approx \frac{{e}^{\mathrm{iHt}}}{\beta}.& \left(77\right)\end{array}$
As R_{k }is close to unitary, the success probability ≈1/β^{2 }may be boosted using oblivious amplitude amplification. When β=2, a single round of oblivious amplitude amplification suffices to boost the success probability to 1−(δ). Thus, one can choose ln 2≤t=(1) such that β=2. If one desires t<ln 2, β may be decreased by appending a singlequbit ancilla and noting that 0e^{iθX}0=cos θ≤1. Thus simulation is accomplished with the circuit
TTS=TTS_{β=2}(REF_{{right arrow over (α)}b}⊗I_{s})TTS_{β=2}^{†}(REF_{{right arrow over (α)}b}⊗I_{s}(TTS_{β=2},REF_{{right arrow over (α)}b}=I_{{right arrow over (α)}b}−200_{{right arrow over (α)}b}. (78)
This approximates timeevolution by e^{−iHt }with error ∥(0_{{right arrow over (α)}b}⊗I_{s})TTS(0_{{right arrow over (α)}b}⊗I_{s})−e^{−iHt}∥=(δ). In order to simulate evolution e^{−iHT }by longer times T>t, one can apply TTS^{T/t}− here t=⊖(1) is chosen such that T/t is an integer. The overall error
ϵ=∥(0_{{right arrow over (α)}b}⊗I_{s})TTS^{T/t}(0_{{right arrow over (α)}b}⊗I_{s})−e^{−iHT}∥=(Tδ), (79)
and success probability 1−(ϵ) may thus be controlled by choosing the error of each segment to
$\mathrm{be}\ue89e\phantom{\rule{0.8em}{0.8ex}}\ue89e\delta =\ue52e\ue8a0\left(\frac{\u03f5}{T}\right).$
This requires a truncation order of
$K=\ue52e\ue8a0\left(\frac{\mathrm{log}\ue8a0\left(\alpha \ue89e\phantom{\rule{0.3em}{0.3ex}}\ue89eT/\u03f5\right)}{\mathrm{loglog}\ue8a0\left(\alpha \ue89e\phantom{\rule{0.3em}{0.3ex}}\ue89eT/\u03f5\right)}\right).$
One may drop the implicit assumption that ∥H∥≤1, by rescaling H→H/α, for some normalization constant α≥∥H∥. Thus simulation of e^{−iHt }requires
$\ue52e\ue8a0\left(\alpha \ue89e\phantom{\rule{0.3em}{0.3ex}}\ue89eT\ue89e\frac{\mathrm{log}\ue8a0\left(\alpha \ue89e\phantom{\rule{0.3em}{0.3ex}}\ue89eT/\u03f5\right)}{\mathrm{loglog}\ue8a0\left(\alpha \ue89e\phantom{\rule{0.3em}{0.3ex}}\ue89eT/\u03f5\right)}\right)$
queries to CHAM. Note that the gate cost of all queries to COEF_{K }at
$\ue52e\ue8a0\left(\alpha \ue89e\phantom{\rule{0.3em}{0.3ex}}\ue89eT\ue89e\frac{\mathrm{log}\ue8a0\left(\alpha \ue89e\phantom{\rule{0.3em}{0.3ex}}\ue89eT/\u03f5\right)}{\mathrm{loglog}\ue8a0\left(\alpha \ue89e\phantom{\rule{0.3em}{0.3ex}}\ue89eT/\u03f5\right)}\right)$
and that of REF_{{right arrow over (α)}b }at
$\ue52e\ue8a0\left({n}_{a}\ue89e\alpha \ue89e\phantom{\rule{0.3em}{0.3ex}}\ue89eT\ue89e\frac{\mathrm{log}\ue8a0\left(\alpha \ue89e\phantom{\rule{0.3em}{0.3ex}}\ue89eT/\u03f5\right)}{\mathrm{loglog}\ue8a0\left(\alpha \ue89e\phantom{\rule{0.3em}{0.3ex}}\ue89eT/\u03f5\right)}\right),$
is typically dominated by the gate cost of all applications of CHAM.
The ancilla overhead of the truncated Taylor series algorithm, at
${n}_{s}+\ue52e\ue8a0\left({n}_{a}\ue89e\frac{\mathrm{log}\ue8a0\left(1/\u03f5\right)}{\mathrm{loglog}\ue8a0\left(1/\u03f5\right)}\right)$
qubits, may be significantly improved by choosing the sequence of unitaries in the compression gadget Lemma 2 of Section III C to be U_{j}=−iHAM. This straightforwardly furnishes the following result.
Corollary 1 (Hamiltonian simulation by a compressed truncated Taylor series). Let a timeindependent Hamiltonian H be encoded in standardform with normalization α and n_{s}+n_{α }qubits, as per Eq. (5). Then the truncated Taylor series algorithm approximates the timeevolution operator e^{−iHt }for any αt≤ln2 to error ϵ using
$1.\ue89e\phantom{\rule{0.8em}{0.8ex}}\ue89e\mathrm{Queries}\ue89e\phantom{\rule{0.8em}{0.8ex}}\ue89e\mathrm{to}\ue89e\phantom{\rule{0.8em}{0.8ex}}\ue89e\mathrm{HAM}\ue89e\text{:}\ue89e\phantom{\rule{0.8em}{0.8ex}}\ue89e\ue52e\ue8a0\left(\frac{\mathrm{log}\ue8a0\left(1/\u03f5\right)}{\mathrm{loglog}\ue8a0\left(1/\u03f5\right)}\right).\text{}\ue89e2.\ue89e\phantom{\rule{0.8em}{0.8ex}}\ue89e\mathrm{Qubits}\ue89e\text{:}\ue89e\phantom{\rule{0.8em}{0.8ex}}\ue89e{n}_{s}+\ue52e\ue8a0\left({n}_{a}+\mathrm{loglog}\ue8a0\left(1/\u03f5\right)\right).\text{}\ue89e3.\ue89e\phantom{\rule{0.8em}{0.8ex}}\ue89e\mathrm{Primitive}\ue89e\phantom{\rule{0.8em}{0.8ex}}\ue89e\mathrm{gates}\ue89e\text{:}\ue89e\phantom{\rule{0.8em}{0.8ex}}\ue89e\ue52e\ue8a0\left(\left({n}_{a}+\mathrm{loglog}\ue8a0\left(1/\u03f5\right)\right)\ue89e\frac{\mathrm{log}\ue8a0\left(1/\u03f5\right)}{\mathrm{loglog}\ue8a0\left(1/\u03f5\right)}\right).$
For longertime simulations e^{−iHT }of duration T>t, Corollary l is applied αT/ln (2) times, each with error
$\ue52e\ue8a0\left(\frac{\u03f5}{\alpha \ue89eT}\right).$
This leads a query complexity
$\ue52e\ue8a0\left(\alpha \ue89eT\ue89e\frac{\mathrm{log}\ue8a0\left(\alpha \ue89eT/\u03f5\right)}{\mathrm{loglog}\ue8a0\left(\alpha \ue89eT/\u03f5\right)}\right).$
Though the compressed algorithm is still worse than the quantum signal processing approach, which uses n_{s}+(n_{α}) qubits, the technique is applicable to simulating timedependent Hamiltonians, as demonstrated in Section III.
XI. Error from Truncating and Discretizing the Dyson Series
In this section, the proof of Lemma 1 is completed for the error from truncating the Dyson series at order K, and the error from approximating its terms, which are timeordered integrals, with Riemann sums. These results provide a rigorous upper bound on the error of timedependent Hamiltonian simulation. Let D_{k }be the k^{th }term in the Dyson expansion, and let B_{k }be the Riemann sum of D_{k }with each dimension discretized into M=t/Δ segments.
$\begin{array}{cc}\ue533\ue8a0\left[{e}^{i\ue89e{\int}_{0}^{t}\ue89eH\ue8a0\left(s\right)\ue89e\mathrm{ds}}\right]=\sum _{k=0}^{\infty}\ue89e\phantom{\rule{0.3em}{0.3ex}}\ue89e{\left(i\right)}^{k}\ue89e{D}_{k}=\underset{M\to \infty}{\mathrm{lim}}\ue89e\sum _{k=0}^{\infty}\ue89e\phantom{\rule{0.3em}{0.3ex}}\ue89e\frac{{\left(\mathrm{it}\right)}^{k}}{{M}^{k}}\ue89e{B}_{k},\text{}\ue89e{D}_{k}:=\frac{1}{k!}\ue89e{\int}_{0}^{t}\ue89e\cdots \ue89e{\int}_{0}^{t}\ue89e\ue533\ue8a0\left[H\ue8a0\left({t}_{1}\right)\ue89e\phantom{\rule{0.8em}{0.8ex}}\ue89e\cdots \ue89e\phantom{\rule{0.8em}{0.8ex}}\ue89eH\ue8a0\left({t}_{k}\right)\right]\ue89e{d}^{k}\ue89et,\text{}\ue89e{B}_{k}:=\sum _{0\le {m}_{k}<\cdots <{m}_{1}<M}\ue89eH\ue8a0\left({m}_{k}\ue89e\Delta \right)\ue89e\phantom{\rule{0.8em}{0.8ex}}\ue89e\cdots \ue89e\phantom{\rule{0.8em}{0.8ex}}\ue89eH\ue8a0\left({m}_{1}\ue89e\Delta \right).& \left(80\right)\end{array}$
It is now proven that the bounds on the ϵ_{1 }term, which is the error due to truncating the Dyson series at order K.
Lemma 5. Let H(s): →^{N×N }be differentiable on the domain [0,t]. For any ϵ_{1 }∈ [0, 2^{−e}], an approximation to the time ordered operator exponential of −iH(s) can be constructed such that
$\uf605\ue533\ue8a0\left[{e}^{i\ue89e{\int}_{0}^{t}\ue89eH\ue8a0\left(s\right)\ue89e\mathrm{ds}}\right]\sum _{k=0}^{K}\ue89e\phantom{\rule{0.3em}{0.3ex}}\ue89e{\left(i\right)}^{k}\ue89e{D}_{k}\uf606\le {\u03f5}_{1},\text{}\ue89e1.\ue89e\phantom{\rule{0.8em}{0.8ex}}\ue89eK\ge \lceil 1+\frac{2\ue89e\mathrm{ln}\ue8a0\left(1/{\u03f5}_{1}\right)}{\mathrm{lnln}\ue8a0\left(1/{\u03f5}_{1}\right)+1}\rceil $
$2.\ue89e\phantom{\rule{0.8em}{0.8ex}}\ue89e{\mathrm{max}}_{s}\ue89e\uf605H\ue8a0\left(s\right)\uf606\ue89et\le \mathrm{ln}\ue89e\phantom{\rule{0.3em}{0.3ex}}\ue89e2$
if one takes
Proof of Lemma 5. We start by bounding D_{k}∥.
$\begin{array}{cc}\uf605{D}_{k}\uf606=\frac{1}{k!}\ue89e\uf605{\int}_{0}^{t}\ue89e\cdots \ue89e{\int}_{0}^{t}\ue89e\ue533\ue8a0\left[H\ue8a0\left({t}_{1}\right)\ue89e\phantom{\rule{0.8em}{0.8ex}}\ue89e\cdots \ue89e\phantom{\rule{0.8em}{0.8ex}}\ue89eH\ue8a0\left({t}_{k}\right)\right]\ue89e{d}^{k}\ue89et\ue89e\phantom{\rule{0.em}{0.ex}}\uf606\le \hspace{1em}\ue89e\hspace{1em}\hspace{1em}\frac{1}{k!}\ue89e\uf605{\int}_{0}^{t}\ue89e\cdots \ue89e{\int}_{0}^{t}\ue89e\prod _{j=1}^{k}\ue89e\phantom{\rule{0.3em}{0.3ex}}\ue89e\uf605H\ue8a0\left({t}_{j}\right)\uf606\ue89e{d}^{k}\ue89et\uf606\le \frac{{\left({t\ue89e\mathrm{max}}_{s}\ue89e\uf605H\ue8a0\left(s\right)\uf606\right)}^{k}}{k!}.& \left(81\right)\end{array}$
At this point, the proof is identical to the timeindependent case as max_{s }∥H(s)∥ is independent of time. Thus using Stirling'"'"'s approximation and assuming K≥2 max_{s }∥H(s)∥t,
$\begin{array}{cc}{\u03f5}_{1}=\uf605\ue533\ue8a0\left[{e}^{i\ue89e{\int}_{0}^{t}\ue89eH\ue8a0\left(s\right)\ue89e\mathrm{ds}}\right]\sum _{k=0}^{K}\ue89e\phantom{\rule{0.3em}{0.3ex}}\ue89e{\left(i\right)}^{k}\ue89e{D}_{k}\uf606\le \sum _{k=K+1}^{\infty}\ue89e\phantom{\rule{0.3em}{0.3ex}}\ue89e\uf605{D}_{k}\uf606\le \sum _{k=K+1}^{\infty}\ue89e\frac{{\left({t\ue89e\mathrm{max}}_{s}\ue89e\uf605H\ue8a0\left(s\right)\uf606\right)}^{k}}{k!}\le \frac{{\left({t\ue89e\mathrm{max}}_{s}\ue89e\uf605H\ue8a0\left(s\right)\uf606\right)}^{K+1}}{\left(K+1\right)!}\ue89e\sum _{k=K+2}^{\infty}\ue89e{\left(1/2\right)}^{kK1}=\frac{{\left({t\ue89e\mathrm{max}}_{s}\ue89e\uf605H\ue8a0\left(s\right)\uf606\right)}^{K+1}}{\left(K+1\right)!}\le {\left(\frac{\left({\mathrm{te}\ue89e\mathrm{max}}_{s}\ue89e\uf605H\ue8a0\left(s\right)\uf606\right)}{K+1}\right)}^{K+1}& \left(82\right)\end{array}$
Now one can find that this in turn is less than ϵ_{1 }if max_{s }∥H(s)∥te<min{ln(1/ϵ_{1}),eln2}≤1 given that, ϵ_{1}≤2^{−e }and
$\begin{array}{cc}K\ge \mathrm{max}\ue89e\left\{1+\frac{\mathrm{ln}\ue8a0\left(1/{\u03f5}_{1}\right)}{W\ue8a0\left(\frac{\mathrm{ln}\ue8a0\left(1/{\u03f5}_{1}\right)}{{\mathrm{max}}_{s}\ue89e\uf605H\ue8a0\left(s\right)\uf606\ue89e\mathrm{te}}\right)},2\ue89e\underset{s}{\mathrm{max}}\ue89e\uf605H\ue8a0\left(s\right)\uf606\ue89e\uf603t\uf604\right\},& \left(83\right)\end{array}$
where W is the LambertW function. Using the fact that for x≥1, W(x)≥(ln(x)+1)/2 and ln(eln2)<1 one can obtain the simpler bound
$\begin{array}{cc}K=\lceil 1+\frac{2\ue89e\mathrm{ln}\ue8a0\left(1/{\u03f5}_{1}\right)}{\mathrm{lnln}\ue8a0\left(1/{\u03f5}_{1}\right)+1}\rceil =\ue52e\ue8a0\left(\frac{\mathrm{ln}\ue8a0\left(1/{\u03f5}_{1}\right)}{\mathrm{lnln}\ue8a0\left(1/{\u03f5}_{1}\right)}\right).& \left(84\right)\end{array}$
It is now proven that the bounds on the ϵ_{2 }term for the error from approximating the Dyson series with its Riemann sum.
Lemma 6. Let H(s): →^{N×N }be differentiable on the domain [0,t]. Let one also define the quantities
$\u3008\uf605\stackrel{.}{H}\uf606\u3009:=\frac{1}{t}\ue89e{\int}_{0}^{t}\ue89e\uf605\frac{\mathrm{dH}\ue8a0\left(s\right)}{\mathrm{ds}}\uf606$
ds. For integer K≥0 and ϵ_{2}>0,
$\uf605\sum _{k=0}^{K}\ue89e{\left(i\right)}^{k}\ue89e{D}_{k}\sum _{k=0}^{K}\ue89e{\left(i\ue89e\frac{t}{M}\right)}^{k}\ue89e{B}_{k}\uf606\le {\u03f5}_{2},$
by choosing any M such that
$1.\ue89e\phantom{\rule{0.8em}{0.8ex}}\ue89eM\ge \frac{{t}^{2}}{{\u03f5}_{2}}\ue89e4\ue89e{e}^{{\mathrm{max}}_{s}\ue89e\uf605H\ue8a0\left(s\right)\uf606\ue89et}\ue8a0\left(\u3008\uf605\stackrel{.}{H}\uf606\u3009+{\mathrm{max}}_{s}\ue89e{\uf605H\ue8a0\left(s\right)\uf606}^{2}\right),\text{}\ue89e2.\ue89e\phantom{\rule{0.8em}{0.8ex}}\ue89eM\ge {K}^{2}.$
Proof of Lemma 6. One first expands the timeordered evolution operator using the Dyson series. One can then examine the error incurred in evaluating a given order of the Dyson series for a small hypercubic region of sidelength Δ=t/M. One can then upper bound the maximum number of such hypercubes within the allowed volume and use the triangle inequality to argue that the error is the product of the number of such hypercubes and the maximum error per hypercube.
Since H(s) is a differentiable function it holds from Taylor'"'"'s theorem that for any δ«1 and computational basis states x,y,
$\begin{array}{cc}\u3008x\uf604\ue89eH\ue8a0\left(s+\delta \right)\ue89e\uf603y\u3009=\u3008x\uf604\ue89eH\ue8a0\left(s\right)\ue89e\uf603y\u3009+\delta \ue89e\u3008x\uf604\ue89e\stackrel{.}{H}\ue8a0\left(s\right)\ue89e\uf603y\u3009+o\ue8a0\left(\underset{s}{\mathrm{max}}\ue89e{\uf605\stackrel{.}{H}\ue8a0\left(s\right)\uf606}_{\mathrm{max}}\ue89e\delta \right).& \left(85\right)\end{array}$
Since computational basis states form a complete orthonormal basis it follows through norm inequalities that
H(s+δ)=H(s)+δH(s)+o(∥H(s)∥_{max }N^{2}δ). (86)
One then has from Taylor'"'"'s theorem and the triangle inequality that
$\begin{array}{cc}\uf605H\ue8a0\left(s+\delta \right)H\ue8a0\left(s\right)\uf606=\uf605\sum _{j=1}^{r}\ue89e\left[H\ue8a0\left(s+j\ue89e\delta /r\right)H\ue8a0\left(s+\left[j1\right]\ue89e\delta /r\right)\right]\uf606\le \uf605\sum _{j=1}^{r}\ue89e\stackrel{.}{H}\ue8a0\left(s+\left[j1\right]\ue89e\delta /r\right)\ue89e\delta /r\uf606+r\ue8a0\left[o\ue8a0\left(\underset{s}{\mathrm{max}}\ue89e{\uf605\stackrel{.}{H}\ue8a0\left(s\right)\uf606}_{\mathrm{max}}\ue89e{N}^{2}\ue89e\delta /r\right)\right]\le {\int}_{0}^{\delta}\ue89e\uf605\stackrel{.}{H}\ue8a0\left(s\right)\uf606\ue89e\mathrm{ds}+r\ue8a0\left[o\ue8a0\left(\underset{s}{\mathrm{max}}\ue89e{\uf605\stackrel{.}{H}\ue8a0\left(s\right)\uf606}_{\mathrm{max}}\ue89e{N}^{2}\ue89e\delta /r\right)\right].& \left(87\right)\end{array}$
Since this equation holds for all r, it also holds in the limit as r approaches infinity. Therefore
∥H(s+δ)−H(s)∥≤∫_{0}^{δ}∥{dot over (H)}(s)∥ds. (88)
Next, let one consider the error in approximating the integral over a hypercube to lowest order and let one define the hypercube to be C with x_{1}, . . . , x_{q }being the corner of the hypercube with smallest norm. First note that in general if A_{j }and B_{j }are a sequence of bounded operators and ∥·∥ is a submultiplicative norm then it is straight forward to show using an inductive argument that for all positive integer q.
$\begin{array}{cc}\uf605\prod _{j=1}^{q}\ue89e\phantom{\rule{0.3em}{0.3ex}}\ue89e{A}_{j}\prod _{j=1}^{q}\ue89e{B}_{j}\uf606\le \sum _{k=1}^{q}\ue89e\left(\prod _{j=1}^{k1}\ue89e\uf605{A}_{j}\uf606\right)\ue89e\uf605{A}_{k}{B}_{k}\uf606\ue89e\left(\prod _{j=k+1}^{q}\ue89e\uf605{B}_{j}\uf606\right).& \left(89\right)\end{array}$
By applying this in combination with Eq. (88) to region C, the error induced is
$\begin{array}{cc}\uf605{\int}_{C}\ue89eH\ue8a0\left({x}_{1}+{y}_{1}\right)\ue89e\phantom{\rule{0.8em}{0.8ex}}\ue89e\dots \ue89e\phantom{\rule{0.8em}{0.8ex}}\ue89eH\ue8a0\left({x}_{q}+{y}_{q}\right)\ue89e{d}^{q}\ue89ey{\Delta}^{q}\ue89e\prod _{j=1}^{q}\ue89eH\ue8a0\left({x}_{j}\right)\uf606\le {\int}_{C}\ue89e\uf605\prod _{j=1}^{q}\ue89eH\ue8a0\left({x}_{j}+{y}_{j}\right)\prod _{j=1}^{q}\ue89eH\ue8a0\left({x}_{j}\right)\uf606\ue89e{\mathrm{dy}}^{q},\le \sum _{k=1}^{q}\ue89e\left(\prod _{j=1}^{k1}\ue89e{\int}_{0}^{\Delta}\ue89e\uf605H\ue8a0\left({x}_{j}+s\right)\uf606\ue89e\mathrm{ds}\right)\ue89e\phantom{\rule{0.3em}{0.3ex}}\ue89e\left({\int}_{0}^{\Delta}\ue89e{\int}_{0}^{s}\ue89e\uf605\stackrel{.}{H}\ue8a0\left({x}_{k}+y\right)\uf606\ue89e\mathrm{dyds}\right)\ue89e\phantom{\rule{0.3em}{0.3ex}}\ue89e\left(\prod _{j=k+1}^{q}\ue89e{\int}_{0}^{\Delta}\ue89e\uf605H\ue8a0\left({x}_{j}\right)\uf606\ue89e\mathrm{ds}\right)\le \sum _{k=1}^{q}\ue89e\left(\prod _{j=1}^{k1}\ue89e{\int}_{0}^{\Delta}\ue89e\uf605H\ue8a0\left({x}_{j}+s\right)\uf606\ue89e\mathrm{ds}\right)\ue89e\left({\int}_{0}^{\Delta}\ue89e{\int}_{0}^{s}\ue89e\uf605\stackrel{.}{H}\ue8a0\left({x}_{k}+y\right)\uf606\ue89e\mathrm{dyds}\right)\ue89e\phantom{\rule{0.3em}{0.3ex}}\ue89e\left(\prod _{j=k+1}^{q}\ue89e{\int}_{0}^{\Delta}\ue89e\alpha \ue89e\phantom{\rule{0.3em}{0.3ex}}\ue89e\mathrm{ds}\right)\le {\left(\alpha \right)}^{q1}\ue89e\Delta \ue89e\sum _{k=1}^{q}\ue89e\left(\prod _{j\ne k}^{q}\ue89e{\int}_{0}^{\Delta}\ue89e\mathrm{ds}\right)\ue89e\left({\int}_{0}^{\Delta}\ue89e\uf605\stackrel{.}{H}\ue8a0\left({x}_{k}+s\right)\uf606\ue89e\mathrm{ds}\right)& \left(90\right)\\ \phantom{\rule{4.4em}{4.4ex}}\ue89e\mathrm{where}\ue89e\phantom{\rule{0.8em}{0.8ex}}\ue89e\alpha :={\mathrm{max}}_{s}\ue89e\uf605H\ue8a0\left(s\right)\uf606.& \phantom{\rule{0.3em}{0.3ex}}\end{array}$
There are two regions in the problem. The first region, which one calls the bulk, is the region that satisfies all the constraints of the problem namely bulk :={(t_{1}, . . . , t_{q}):└t_{1}/Δ┘>. . . >└t_{q}/Δ┘}. Thus for any index x_{1}, . . . , x_{q }to a hypercube in the bulk, the ordering of terms H(x_{1}+t_{1}) . . . H(x_{q}+t_{q}) in the integrand of Eq. (90) is fixed. The second region is called the boundary which is the region in which the hypercubes used in the Riemann sum would stretch outside the allowed region for the integral. Since one can approximate the integral to be zero on all hyperc