Low-cost quantum circuits for classically intractable instances of the Hamiltonian dynamics simulation problem

We develop circuit implementations for digital-level quantum Hamiltonian dynamics simulation algorithms suitable for implementation on a reconfigurable quantum computer, such as trapped ions. Our focus is on the codesign of a problem, its solution, and quantum hardware capable of executing the solution at the minimal cost expressed in terms of the quantum computing resources used, while demonstrating the solution of an instance of a scientifically interesting problem that is intractable classically. The choice for Hamiltonian dynamics simulation is due to the combination of its usefulness in the study of equilibrium in closed quantum mechanical systems, a low cost in the implementation by quantum algorithms, and the difficulty of classical simulation. By targeting a specific type of quantum computer and tailoring the problem instance and solution to suit physical constraints imposed by the hardware, we are able to reduce the resource counts by a factor of 10 in a physical-level implementation and a factor of 30–60 in a fault-tolerant implementation over state-of-the-art.


I. INTRODUCTION
Quantum supremacy is a computational experiment designed to demonstrate a computational capability of a quantum machine that cannot be matched by a classical computer.It is highly relevant to this paper, since we too focus on a quantum computation of the size that cannot be performed via classical means.Quantum supremacy is an important milestone in the development of quantum computers.Multiple IT giants are targeting quantum supremacy [1], and it is perhaps reasonable to anticipate that a successful demonstration may be obtained within at most a few years.
Once quantum supremacy is demonstrated, a next step is to go beyond what the supremacy would have achieved.The supremacy experiment proposed in [2], in particular, reduces to the execution of a random quantum circuit on a quantum computer that is too large for a classical computer to cope with simulating.In this work, rather than focusing on an artificial problem designed purely for demonstrating quantum supremacy [2], we focus on the selection of a known computational problem and a specific input instance, and developing a quantum circuit computing the answer for the selected problem/instance pair such that this quantum circuit relies on the least quantum resources and can be suitable for the execution on near-term quantum computers, while, to the best of our knowledge, the problem/instance pair requires a classically intractable computation.Such a computation constitutes a qualitative step forward, where a quantum computer can now be thought of as being a tool in the solution of a problem rather than the focus of the study.A more advanced demonstration past the one we are reporting in this work could target the solution of a problem with a commercial value.
In the following, we introduce the problem and the specific instance of this problem that we are proposing as satisfying the conditions outlined in the previous paragraph.We next describe our solution, including numerous techniques used to improve quantum computing resources used.Finally, we discuss the quality of our solution, expressed in terms of quantum circuit gate count and depth, and viewed through the lens of comparisons to prior and similar-spirited work.We stress that we develop complete and fully specified quantum circuits as a part of this study.

II. PROBLEM
We consider the problem of simulating Hamiltonian dynamics for time t, accurate to within the error ε, where the target Hamiltonian H is defined as follows: where σ i x , σ i y , and σ i z denote Pauli x, y, and z matrices acting on the qubit i, d i ∈ [−1, 1] are chosen uniformly at random, G is a graph describing the two-qubit interactions, E(G) is the set of its edges, and n is the number of qubits.Such Hamiltonian is known as the Heisenberg Hamiltonian over a graph G with a random disorder in the z direction.It has been studied in [3][4][5] in the context of many-body localization.
We chose n to be in the range 50 to 100.This is because the largest quantum circuit (vector state as opposed to full unitary) simulations demonstrated to date were achieved with dozens of qubits over a circuit depth of about 40 [6,7].This means that with anywhere more than 50 qubits and depth in excess of, say, 200, the problem of quantum circuit simulation may likely become intractable for a classical computer with the simulation techniques such as [6,7].Note that our circuits, despite numerous optimizations applied, remain substantively deeper than those considered in [6,7] motivating our choice to consider reducing the number of qubits from as many as 144 for a shallow circuit with depth 27 [6] down to 50 at the cost of significantly extending the anticipated length of the computation.We also note that the underlying qubit-to-qubit connectivity pattern in our circuits does not appear to allow breaking the qubit interaction graph into two components by cutting a small number of edges, which lies at the core of efficient simulations such as [6,7].The smallest number of qubits we chose to consider, 50, exceeds the number 42 used in the best high-depth state-vector type simulation [8].Finally, we highlight that the largest numerical simulation of the kind of Hamiltonian we consider is restricted to just 22 qubits [3], whereas our work focuses on the Hamiltonians over at least 50 qubits.
We chose the evolution time t = 2d, where d is the diameter of graph G, similarly to [9].The motivation behind such choice is as follows: it takes time at most d [10] for quantum information to propagate from any node in the underlying graph G to any other node, as such, one may expect to enter a highly entangled simulation regime by the simulation time of d.This means that our simulation spends at least half the time evolving in the regime that we believe is difficult to simulate classically.Note that due to the use of the product formula approach [11,12] in our simulations, circuits for other selection of time t can be developed as effortlessly as changing the number of times a certain block operation is applied.
A previous study [9] showed that the product formula algorithm(s) [11,12] for Hamiltonian simulation with a heuristic bound yields best practical results.We note that choosing a small diameter graph G, such as what we do next, is natural, given the discussion in the previous paragraph.However, small diameter graphs require large number of edges, and the number of edges in the graph G is directly proportional to the circuit complexity of the single stage of product formula.Thus, the balance between graph diameter and the number of edges has to be chosen carefully so as to minimize quantum computational resources while maximizing the expected classical difficulty of the simulation.
We chose graph G to be a minimal distance k-regular graph over n nodes.We select k such that the effort required to develop the individual two-qubit gates for all nk 2 interactions prescribed by the graph in a technology such as the trapped ions is not large [13], while keeping the graph distance small (resulting in a short enough time t for the evolution before Hamiltonian dynamics simulation enters a regime that is believed to be classically difficult) and the number of qubits as large as possible (while keeping it to between 50 and 100).In practice, we selected the following values of k: 3, 4, 5, 6, and 7.The specific regular graphs G considered in our work can be described by the respective degree-diameter-nodes triple (k−d−n), as follows: (3−5−70) Alegre-Fiol-Yebra graph [14], (4−4−98) graph by Exoo [15], (5−3−72) graph by Exoo [15], and (7−2−50) Hoffman-Singleton graph [16].The largest known 6-regular distance-2 graph has 32 nodes (less than 50 targeted in our work), and the largest known 6-regular distance-3 graph has 110 nodes (more than 100 targeted in this work) [17].Therefore, we did not consider 6-regular graphs.Note that our goal was to select a graph with a large number of nodes, small degree, and small distance.A plenty of such graphs can be developed so long as one chooses degree and distance parameters above the minimums known for a given n, selects n smaller than the known maximum for the fixed degree and diameter, and expands the attention to graphs other than the regular kind.Our selection of graphs is very restrictive so as to narrow down the set of specific graphs explicitly considered in our work to a few.We believe the performance of the Hamiltonian simulation over similar graphs to those considered can be similar.
We chose the approximation error ε = 0.001, measured as the spectral norm distance between the target ideal evolution and the evolution obtained by running our quantum simulation circuits, as in [9].This choice of the error value is largely arbitrary, but some choice is necessary to explicitly construct the circuits.
cnot count optimal implementation of the Heisenberg interaction up to a global phase of e −iπ/4 .This is a special case of the circuit shown in Fig. 6 of [20].Note that this circuit requires 3 cnot gates.
FIG. 2. Implementation of the Heisenberg interaction, optimal in the number of real-valued degrees of freedom, up to a global phase of e −ia .Implementation of the controlled-z a can be substituted from Figure 3.

III. SOLUTION
Suzuki-Trotter formula based algorithms with empirical bound showed themselves as a potent candidate among quantum algorithms that simulate Hamiltonian dynamics [9,11,12].We therefore chose to focus exclusively on this type of algorithms.Specifically, for a Hamiltonian of the form H = j α j H j , we approximate the evolution operator according to where λ := −it/r and with p k := 1/(4 − 4 1/(2k−1) ) for k > 1 [11].We selected 4th (k=2) and 6th (k=3) order formula versions since in our case they achieved the best results.For a given order formula, the algorithm applies the operation S 4 /S 6 r times, where S 4 /S 6 contains 10/50 repetitions of the product of exponentials of the individual terms of the Hamiltonian H, and r is the number of repetitions of S 4 /S 6 for a given selection of the evolution time t, desired accuracy ε, underlying graph G, and the selection of random disorders d i (1).We note that for a fixed ε, graph G, and a given selection of random disorders, r is proportional to a growing function of t.Explicit bounds on the value r are O(t 1+1/4 ) for the 4th order formula and O(t 1+1/6 ) for the 6th order formula [12].
Since the selection of t in our implementation is t = 2d, it is important to minimize the diameter of the underlying graph G, which explains our focus on the low-diameter graphs.The cost of S 4 /S 6 can be described as 10X/50X, where X is the gate cost of the implementation of a single stage of the product of individual terms in the target Hamiltonian.In our work, we obtained physical-level implementations of the X stage with the cost of 3nk 2 CNOT gates and the twoqubit depth 3k to 3k+3, and fault-tolerant cost of k(4n+log , where parameter a is defined as the depth of the approximation of an r z (see below for details).This means that all figures of merit are linearly dependent on the degree k of the graph G, and therefore to achieve the best performance, the degree k must be minimized.
To optimize the depth of our circuits, we employed a customized version of the implementation of Vizing's theorem [18], which guarantees circuit depth k or k+1 for the implementation of a single stage of the product of exponentials of the target Hamiltonian.Our modification includes additional heuristic that randomly reorders the list of edges of the graph G in an attempt to find a depth-k layout when a depth k+1 layout is found, though in practice we found this to be of limited use.Our modification also readily accepts a manual input in case a depth-k layout is known, , where although we chose not to employ this method for the circuits we consider in this paper, since in general the problem of finding a depth-k layout, if at all exists, is NP-complete [19] and thus it is unlikely that the appropriate manual input would be known for a generic graph.We implemented the Heisenberg interaction, exp[−ia( σ i x σ j x + σ i y σ j y + σ i z σ j z )], using two different circuits, depending on whether we focus on saving quantum resources in a pre-fault tolerant (physical-level) or a fault-tolerant implementation.In particular, for the physical-level implementation, we used the circuit shown in Figure 1 in order to minimize the most expensive two-qubit gates, whereas for the fault-tolerant implementation, we used the circuit shown in Figure 2 in order to minimize the most expensive t gates.By directly synthesizing the Heisenberg interaction, compared to the standard Pauli-matrix basis approach, for instance employed in [9], we save 50% of the cost in the pre-fault tolerant implementation (evidenced through the reduction of the cnot gate count from 6 down to 3) and almost 66% of the cost in the fault-tolerant implementation (evidenced through the reduction from 3 r z gates down to 1 r z and 4 t gates).Because our construction directly implements the Heisenberg interaction, as opposed to an approximate implementation with XX, YY, and ZZ interactions that arise from the standard approach, our implementation also performs better at the algorithmic level, i.e., we do not need as large value of r as in the standard approach to keep the overall error level down that incurs from the approximate Pauli-basis implementation.
We laid out the Heisenberg interaction terms in the circuit implementation of the product formula algorithm in alternate orders-forward and reverse-to ensure we obtain maximal gain from the application of the circuit optimizer [21].We furthermore modified the original optimizer [21] to ensure it can natively handle controlled-z a = cz a gates and apply the merging rule cz a (x, y)cz b (x, y) → cz a+b (x, y).Also implemented was the merging of two Heisenberg interactions as per circuit implementation in Figure 1 as an optimization rule, and the capability of being able to handle classically controlled quantum gates.The quality of optimization by the automated optimizer ranged from 7% to 14% in the cnot gate count reduction, and 16% to 24% in the r z count reduction, with the simulations over lower degree graphs yielding a better quality of optimization.Indeed, circuit implementations over graphs with lower degree have a larger proportion of gates at the edges of the circuit implementation of the Hamiltonian terms, and these are the types of gates that often admit optimizations by techniques such as [21].
For fault-tolerant implementations, we decided to break the error budget evenly between the algorithmic errors that arise from the product formula algorithm and the approximation errors that arise from the approximation of r z gates in the Clifford+t basis.We distributed the approximation error budget evenly across all r z gates in the given circuit.We employed three optimization strategies for fault-tolerant implementations: the mixing unitaries approach detailed in [22], the application of equal-angle r z rotations through computing input weight [24], and a combination of gridsynth [26,27] and repeat-until-success (RUS) [25] strategies for the approximation r z gates by Clifford+t circuits.These three strategies are detailed in the next three paragraphs.
The mixing unitaries approach [22] relies on generating a set of four approximating circuits for a given r z (θ) gate, and then applies a randomly drawn approximation from the set of four for every occurrence of the r z (θ) gate in the circuit, subject to a certain probability distribution.We found out through simulation experiments that in practice the mixing unitaries approach achieves less than theoretically possible quadratic improvement in the error.To study practical performance of the mixing strategy, we chose to investigate an amenable sample case of (5−1−6) graph in depth.We varied the per-gate error budget by tweaking the power p of (ε approx./(N rz )) p , where ε approx.is the approximation error budget and N rz is the number of r z gates in the given circuit, since the number of t gates for an approximation sequence scales linearly in the logarithm of the inverse of per-gate error level.An extensive numerical investigation showed that choosing p between 0.86 and 0.89 works well for the direct synthesis method and the values of p between 0.65 and 0.68 represent the advantage by the mixing approach.We therefore chose to use p = (0.65 + 0.68)/2 = 0.665 as the per-gate approximation power for the mixing unitaries method, for all cases we considered.We believe this is a safe assumption since in the larger circuits the mixing strategy is expected to give better results as it takes more time to properly average out the errors.Employing the mixing unitaries strategy allows an estimated t count savings of 25% to 33%, depending on whether the power p = 0.875 or p = 1 is considered as the starting point.
In our implementation of the two-qubit Hamiltonian interaction, we lay out the respective circuitry in parallel with the help of Vizing's theorem.This results in the parallel application of as many as m≤ n 2 r z (θ) gates with equal rotation angles [23].Such transformation can be accomplished directly at the cost of m•Cost(r z (θ)) t gates, but a better approach is to induce this set of gates via the calculation of the input weight sum and the application of r z (2 log(m) θ), r z (2 log(m) −1 θ), ..., r z (2θ), and r z (θ), to the binary digits (qubits) of the integer number ( log(m) +1qubit ket) representing the input weight of the m qubits needing the application of r z (θ) gates, at the cost of 4m + log(m)•Cost(r z (θ)) + Const (where Const is a small constant) t gates [24].Indeed, for the implementations we considered Cost(r z (θ)) ≈ 50 t gates, and thus the saving in the t count is substantial.To induce the input weight calculation at the cost of at most 4m t gates, we use at most m full and half adders.The application of each reduces the number of bits yet to be added by at least 1, and we perform lower digit summations first while doing so in parallel.Since both half and full adder need one relative phase Toffoli gate to be computed (costing 4 t gates each) and only Clifford gates and one measurement to be uncomputed [24], the t count in the calculation of the input weight and proper reset of ancillae does not exceed 4m.The overall optimization of the t count from applying parallel r z gates through the calculation of the input weight ranges from 51% to 60%, being better for simulations of higher degree graphs.Indeed, for higher degree graphs, the product formula algorithm expels a higher fraction of resources on implementing Hamiltonian interactions compared to the resources spent on implementing the random disorder.
We use gridsynth [26,27] to synthesize optimal single-qubit Clifford+t circuits implementing individual r z gates.A better strategy relies on the RUS circuits [25] that use ancilla, measurement, and feedforward to reduce the t count in the implementation of a single r z gate by a factor of 2.5 on average.However, a fully automated implementation of the RUS strategy is unavailable, and the implementation we have is labor-intensive.Thus, we do not compute RUS circuits explicitly, but rather report the respective expected gate counts.The RUS implementation can be easily included into our software as an external and independent package.The result of the overall optimization of the Hamiltonian dynamics simulation circuits by applying the RUS approach ranges between 44% and 50%, with optimization quality favoring smaller order graphs.Indeed, those have smaller parts composed with half and full adders that are not optimized by the RUS.
We note that optimizations described in the last two paragraphs reduce the t count at the cost of introducing a number of cnot gates.We keep track of the cnot gates to make sure their cost does not overwhelm that of the t gates.
To determine the empirical bound on the number r of iterations of the product formula for our problem that aims to address the graphs (3−5−70), (4−4−98), (5−3−72), and (7−2−50), we considered random regular graphs with k = 3, 4, 5, and 7, generated using random matching approach [28], where the number of vertices n ranges from k+1 to 12.For the cases when nk is odd and it is impossible to generate the corresponding random k-regular graph (the number of edge ends can only be even), we take a degree-k random graph with n+1 vertices, and remove a randomly selected vertex as well as all edges leading to it.We then insert edges connecting k/2 non-overlapping pairs of the resulting k degree-(k−1) vertices, chosen at random, provided that the introduction of the new edges does not lead to a multi-edged graph.
Figure 4 and Figure 5 show, for pre-fault tolerant and fault-tolerant implementations, respectively, the r-scaling for the 4th and 6th order formulas for degree k = 3, 4, 5, and 7 random regular graphs.By performing least square linear fitting on the log-log scale we determined the following scaling of r in the pre-fault tolerant case, 4th order formula r 3 =116.3n 0.169 , r 4 =66.4n 0.331 , r 5 =30.1n 0.602 , in the pre-fault tolerant case, 6th order formula r 3 =12.5n 0.564 , r 4 =7.05n0.759 , r 5 =4.24n 0.883 , r 7 =7.11n 0.476 , while for the fault-tolerant case, 4th order formula r 3 =126n 0.215 , r 4 =80.5n 0.328 , r 5 =34.6n 0.620 , and for the fault-tolerant case, 6th order formula r 3 =15.9n 0.493 , r 4 =7.82n 0.751 , r 5 =5.26n 0.826 , r 7 =7.61n 0.490 , where r k is the empirical bound for the degree-k graph with evolution time t = 2d = 10, 8, 6, 4 for graphs with k = 3, 4, 5, 7, respectively.We did not include the scaling for the 4th order formula and k=7, since there were too few points to study and there did not seem to be enough stability in the data.Detailed results showing concrete gate counts for individual cases are available in Table II and Table III for the 4th and 6th order formulas, respectively.Also shown are the gate counts expected for RUS approach [25], where, for each r z approximation, we expect a factor 2.5 reduction in the t count at the cost of the addition of at most (t-count+1) cnot gates and one ancilla.

IV. RESULTS
Table I shows gate counts in the "post-supremacy" Hamiltonian dynamics simulation circuits we are proposing.We show resource counts targeted for both physical-level and fault-tolerant implementations.We note that the counts were obtained while optimizing the circuit depth.We found that they may be improved by about 20% if we optimize the gate counts themselves instead.
For the purpose of comparison to prior work reporting detailed gate counts, we focus on our best result, the Hamiltonian simulation over the (3−5−70) graph.In quantum chemistry, the simulation of FeMoco (the primary cofactor of nitrogenase, which is an enzyme used in the nitrogen fixation) required the circuit with 10 14 t gates over 111 qubits [29].In comparison, our fault-tolerant circuits with 6.8×10 6 and 1.3×10 7 t gates (4th and 6th order product formulas, correspondingly) are significantly shorter, while relying on a comparable number of qubits, 131, and corresponding to solving a task of a similar classical complexity.To factor a 1,024-digit integer number, [30] constructed a circuit with 5.7×10 9 t gates spanning 3,132 qubits.Our circuits are orders of magnitude shorter and operate over a much smaller number of qubits.Recent work [31] required 10 8 t gates to solve a problem in the study of solid-state electronic structure, whereas our t count is only 6.8×10 6 .Finally, the task of a very similar complexity (70 qubits, ε = 0.001, Hamiltonian on a cycle [9]) is solved using 4×10 8 t gates or 7×10 6 cnot gates, showing the advantage of our approach by a factor of 60 in the t count and a factor of 10 in the cnot count.In fact, the circuits we developed to simulate the Hamiltonian evolution are so short, that we hope they may be possible to execute on pre-fault tolerant quantum computers: indeed, all cnot counts are within less than several million.
The two-qubit depth of our quantum circuits is very small, making them particularly suitable for implementations over QIPs with limited T 1 /T 2 coherence times.Our shortest circuit has the two-qubit gate depth of only 25,333 (Table I).
We believe circuit implementations reported in this paper could be particularly relevant to the near-term quantum information processors based on the trapped ions technology.This is because the low-diameter graphs, suitable for the kinds of experiments we considered, have connectivities that require long-range interactions between qubits.Leveraging all-to-all connectivity of the trapped ions QIP, we expect no overhead cost in shuttling quantum information around, a serious point to consider, as has been pointed out in [32].Naturally-provided long coherence time of ion qubits should also be a boon for a pre-fault tolerant implementation.

V. CONCLUSION
In this paper, we synthesized short quantum circuits that aim to solve a scientifically interesting problem.Specifically, we considered the Heisenberg Hamiltonian simulation with a random disorder on a suite of graphs with small diameter.Compared to the previous state of the art, our work shows significant gate count savings, a short circuit depth, while natively relying on the qubit-to-qubit connectivity suitable for the implementation over a reconfigurable quantum computer, such as the trapped ions one.Specifically, we reported a circuit for simulating Hamiltonian dynamics over (3−5−70) graph for the time t = 10 and accurate to within the error ε = 0.001 using at most 648,885 cnot gates in depth 25,333 in a pre-fault tolerant implementation and 6,751,395 t gates in a fault-tolerant implementation.We believe the problem instance we considered is intractable for classical computers and yet the quantum resource estimates are very low, showing the promise for solving interesting problems by a quantum computer in a not-too-distant future.

FIG. 4 .FIG. 5 .
FIG.4.Parameter r scaling data for pre-fault tolerant implementation as a function of the system size n for 4th (left) and 6th (right) order formulas.The black squares, orange circles, and blue triangles denote k = 3, 4, and 5 graphs, respectively.The error bars denote one standard deviation.The solid lines are the best fit power law curves (4) and (5) for 4th and 6th orders, respectively.

TABLE I .
Gate counts for our proposed experiments.

TABLE II .
Gate counts for different graphs using 4th order PF.Graph (k−d−n) Pre-fault tolerant cost Fault-tolerant cost (Gridsynth) Fault-tolerant cost (RUS)

TABLE III .
Gate counts for different graphs using 6th order PF.Graph (k−d−n) Pre-fault tolerant cost Fault-tolerant cost (Gridsynth) Fault-tolerant cost (RUS)