Quantum simulation of tunneling in small systems.

A number of quantum algorithms have been performed on small quantum computers; these include Shor's prime factorization algorithm, error correction, Grover's search algorithm and a number of analog and digital quantum simulations. Because of the number of gates and qubits necessary, however, digital quantum particle simulations remain untested. A contributing factor to the system size required is the number of ancillary qubits needed to implement matrix exponentials of the potential operator. Here, we show that a set of tunneling problems may be investigated with no ancillary qubits and a cost of one single-qubit operator per time step for the potential evolution, eliminating at least half of the quantum gates required for the algorithm and more than that in the general case. Such simulations are within reach of current quantum computer architectures.

Quantum simulations on quantum computers are one of a set of algorithms that give exponential improvement in computational resources relative to the best classical algorithm1234.

Algorithms for studying many types of quantum field theory have been considered in both analog form in which a quantum Hamiltonian (typically manybody or multiple spin) is mapped either directly or via a suitable pulse-sequence to a computational Hamiltonian56789101112131415161718192021 and digital form in which a quantum system's Hamiltonian is split into free and interacting operators, then, using Trotter's formula, is simulated on a quantum computer2223.

Small quantum simulations have already been realized on NMR2425262728293031, atomic3233, ion trap3435363738 and photonic3940 quantum computers in both analog and digital forms.

A concerted effort has also been made to investigate digital quantum particle simulations3 for the simulation of chemical dynamics4142. Algorithms for state preparation43, the simulation of temporal dynamics41 and the measurement of observables44 have all been developed. However, this type of simulation has remained untested due to the large number of gates and/or ancillary qubits needed to compute the kinetic and potential operators4145. In this paper, we focus on reducing the number of gates required for the simulation of temporal dynamics to the bare minimum, while still simulating interesting physics.

For a review of the current state and outlook for quantum simulation, see4647484950 and51.

The standard digital quantum simulation algorithm for a particle on a one-dimensional grid345 encodes the position efficiently in n = log2 N qubits, where N is the size of the lattice of discretized particle locations, x = kΔx, k = 0, ... ,N − 1. The method uses a split operator approach to integrate a Schrödinger equation with a time-independent Hamiltonian that is first-order accurate in the time step, Δt3: Higher order methods that give more accurate time integration have been developed525354, but methods of order higher than two require more gates per time step. For simplicity, we will only consider first-order methods here. States in the qubit Hilbert space where |ψ〉 ∈ {|0〉, |1〉}, represent particle location in the binary representation, |x〉 = Σ∈0,..., 2|ψ〉. The matrix exponential for the kinetic operator is calculated using a quantum Fourier transform (QFT) where F† is a discrete Fourier transform operator and T is diagonal with entries proportional to −q2/2m, and q denotes the Fourier mode wavenumber. The resulting periodic, shift invariant unitary transform, e−, gives an approximation to the time evolution on the lattice due to the kinetic energy operator that is accurate to N'th order in space. This leads to the digital quantum particle simulation algorithm: The QFT takes of order n2 gates to calculate55 and general algorithms implementing the diagonal T and V operators require ancillary qubits56, although it has been shown that the quadratic kinetic energy operator may be computed with n2 two-qubit gates with no ancillary qubits45.

The point of this paper is to describe a quantum particle simulation that can be used as a proof-of-principle demonstration. To do this, it must be implemented in current quantum computer architectures where the number of qubits is limited and the number of gates that may be performed before decoherence destroys the computation is also limited. We consider the special case of square-well potentials. We show that a set of square-well potentials may be implemented with a sole, single-qubit operator and no ancillary qubits (see Methods). This virtually eliminates the calculation of the potential operator. Thus, simulations of important physical phenomena such as tunneling and the evolution of quasi-stable states may be performed with considerably fewer gates than simulations requiring arbitrary potentials. Additionally, we present a reduction of the number of gates required for the diagonal part of the kinetic operator in the case of 4 or fewer qubits.

Results

A Two Qubit Tunneling Simulation

Let us consider the smallest possible tunneling simulation. An n = 2 qubit simulation of N = 4 lattice points may be performed with the circuit

One Time Step of a Two-qubit Digital Quantum Single-particle Tunneling Simulation.

In this circuit and below, Ω and Φ are controlled phase gates applied to qubits i and j and H and Z are a Hadamard gate and a Z-rotation on qubit i, respectively (see Methods). In total, each time step in this tunneling simulation requires 10 operations, 7 single qubit operations and 3 two-qubit operations. The circuit for a single time step is shown above. This circuit implements a double-well potential with the gate acting on the lowest order qubit.

In Fig. 1, we plot lattice occupation probabilities from two simulations with a double-well potential with Δt = 1/10: one a free-particle simulation with v = 0 and the other a tunneling simulation with v = 10. This value of Δt traded off accuracy for gate number (see Methods section for a discussion of simulation errors), but the qualitative dynamics did not change even for more accurate (and costly) simulations. The initial state was |ψ〉 = |01〉, corresponding to a particle in one of the wells. The free-particle probability distribution spreads across all lattice points as it evolves, whereas the particle tunnels from the well at lattice point 1 (|01〉) to the well at lattice point 3 (|11〉) in the tunneling simulation. These results show that differences in the evolution of the probability distribution are evident within 4 time steps. Thus, such a proof-of-principle simulation may be implemented on a quantum computer with 4 × 10 = 40 gates (28 single-qubit and 12 two-qubit).

Figure 1

Particle probability distributions as a function of time for the first four steps of a two qubit simulation for v = 0 (free particle) and v = 10 (particle in double well).

The double-well potential and gray scale used to plot the probabilities are depicted to the right.

Multi-Qubit Quantum Tunneling Simulations

Larger simulations require more gates. For instance, a three-qubit simulation requires 6 gates per QFT (3 singlequbit and 3 two-qubit), 6 gates for the diagonal kinetic energy operator (3 singlequbit and 3 two-qubit), and one single-qubit gate for the potential as shown in the circuit diagram below.

where is shown, representatively, for a double-well potential, but other square-well potentials could be generated by acting on different qubits. The QFT and diagonal kinetic energy operators are

In Fig. 2, we show results from a three-qubit simulation with a double-well potential, where each well is resolved with two lattice points. The time step Δt = 1/5 and v = 5. Here, a tradeoff was also necessary to find a simulation that captured the oscillatory and tunneling time scales and had few gates. The initial state was |ψ〉 = |110〉. Because the initial state only occupies half of one well, oscillatory dynamics are visible within the well. After a few time steps, the oscillatory state tunnels between wells. Oscillatory dynamics are evident within 4 time steps, but 5 time steps are required before the oscillatory state tunnels appreciably to the second well and 7 time steps are needed to see oscillation of the tunneled state. Thus, between 4 × 19 = 76 or 7 × 19 = 133 gates would be required in order to see interesting tunneling effects in such a simulation. Four-qubit simulations can be envisaged using circuits based on similar methodology, although the number of gates would most likely be prohibitive for current quantum computer architectures.

Figure 2

Particle probability distribution as a function of time over ten time steps in a three-qubit double-well simulation.

The potential schematic is shown to the right, as in Fig. 1. Gray scales are as in Fig. 1. The double-well potential has two lattice points per well and, for the given initial state, an oscillation is induced in one of the wells then tunnels to the other well.

Discussion

Only n = log2 N qubits are required for an N lattice-point particle simulation, therefore this algorithm is efficient in the number of qubit resources required. In general, n2 gates are needed for the kinetic operator (see Methods), but for fewer than n = 4 qubits, further efficiency is possible using the basis approach that we outlined. Furthermore, as noted above, our algorithm virtually eliminates the calculation of the potential operator. With very few qubits, interesting tunneling dynamics may be simulated with a gate count that is within reach of current quantum architectures.

Much recent work has been done to understand the resources necessary to perform fully error-corrected quantum particle simulations using various error correction schemes57. We note that, as in the non-fault-tolerant case, the simulations presented here would also require the fewest resources in fully error corrected quantum architectures, since the number of logical qubits and gates is as small as possible in the error-corrected case as well.

Methods

Square-well Potentials

A set of square-well potentials may be implemented with a sole single-qubit operator and no ancillary qubits. To see this, consider the single-qubit Z-rotation on the highest order qubit where v is a parameter, a superscript indicates the qubit to which the operator is applied and σ is the Pauli z-matrix The operator e−, when acting on a lattice state, implements a square-well potential by rotating qubit states with |0〉 (|1〉 resp.) highest order qubit with positive (negative resp.) phase velocity v. The single-qubit operator acting on the next highest order qubit implements a double-square-well potential, and so on, with the last potential in this series implementing a Dirac comb-like potential.

By simulating this class of square-well potentials we reduce the complexity of the potential calculation from a potentially large number of gates and ancillary qubits to just one single-qubit operation. In contrast, the phase kickback algorithm used in41 requires a total of 2n qubits (n for the quantum state and n ancilla qubits) as well as n2 extra gates for a QFT on n ancilla qubits.

Although the main point of this paper is the simplification of the simulation by reducing the potential computation to one gate, further gate reductions can be made to the kinetic operator when considering few qubit systems. This can be done by forming a diagonal Hamiltonian using a basis of diagonal operators. This method is only efficient for n ≤ 4 and scales as 2 − 1. Therefore, for n > 4, ancillary qubits or the Benenti-Strini (BS)45 method (that scales as n2) should be used. However, for small systems of qubits, forming a basis of diagonal operators takes fewer gates. Note that for n = 2, 3 and 4, we would need 4, 9 and 16 (resp.) operators for BS, but at most 3, 7 and 15 (resp.) are required using a basis along the diagonal. Furthermore, fewer multi-qubit operations are necessary with this method, although a three-qubit operator may be necessary for three-qubit simulations and three- and four-qubit operators for four-qubit simulations. Here, we use the method to construct the diagonal operator D for the kinetic energy operator, however, arbitrary potentials could also be constructed this way. Because our goal is to find interesting simulations with few gates, we do not pursue more complex potentials here.

In our two-qubit simulations (see Results), we use a double well potential, . The QFT may be computed with the operators55 where H is a Hadamard operator on the i'th qubit and the controlled-phase gate Ω01 = diag(1, 1, 1, ω) with ω = exp(2πi/4). Note that F† results in a bitswapped Fourier transform (in our notation, F is the inverse Fourier transform matrix).

The kinetic operator is then K = FDF†, with (Note that this operator is also bit-swapped and we have taken m = 1/2).

The diagonal operation D may be achieved up to an overall phase with the operator where the single-qubit operators and Φ01 is a controlled-phase operator on qubits 0 and 1, where . The coefficients (c0 = −1, c1 = −4 and c2 = 4) in the unitary operators Z0, Z1 and Φ01 (resp.) were obtained by noting that the vectors (1, 1,−1,−1), (1,−1,1,−1) and (1, 1, 1,−1) that form their diagonal elements are a basis for zero mean vectors (i.e. neglecting the constant phase proportional to (1, 1, 1, 1)) for .

To calculate the QFT in the three-qubit case, the unitary operators are the single-qubit Hadamard gates and the two-qubit controlled-phase operators Ω = Ω(ω), where ω = exp(2πi/8) (see circuit diagram in Results). To calculate the quadratic diagonal operator in the three-qubit case (see Results), the unitary operators are the single-qubit operators and the two-qubit controlled-phase operators where and Note that, in principle, a seventh three-qubit operator would also be necessary proportional to the vector (1, 1, 1, 1, 1, 1, 1,–1) along the diagonal, but for the (bit-swapped) diagonal of the 8 lattice point simulation, (0, 16, 4, 4, 1, 9, 9, 1), its coefficient is identically zero. This circuit requires a total of 19 gates (10 single-qubit and 9 two-qubit).

Once we have computed the unitary operators for one time step of the kinetic, U(Δt), and potential, U(Δt), evolution, it is straightforward to calculate the exact Hamiltonian evolution: Uexact(Δt) = exp(log(U(Δt))+log(U(Δt))). We then determine RMS errors relative to the exact evolution for the two and three-qubit simulations. These are given in Table 1. The simulations that we present in Figs. 1 and 2 have a total, integrated error of 64% and 20%, respectively. For both simulations, we checked that state evolution was qualitatively similar to the exact simulations for the chosen values of Δt. Clearly, for high-precision simulations, smaller errors would be desirable. To achieve this, one would use smaller values for Δt at the expense of an increase in gate count.

Table 1

RMS Error per time step for a range of values of Δt for the 2 and 3 qubit simulations described in Results

Δt	0.2	0.1	0.05	0.025	0.0125
Error 2 qubit	0.40	0.16	0.06	0.03	0.01
Error 3 qubit	0.08	0.02	0.007	0.003	0.002