Quantum Computation and Visualization of Hamiltonians using Discrete Quantum Mechanics and IBM QISKit

TL;DR

This paper explores methods to translate various quantum systems into a form compatible with quantum computers, specifically using matrix formulations and visualization techniques, to enhance quantum simulation capabilities.

Contribution

It introduces two approaches—position basis and Gaussian basis—for mapping systems onto quantum computers and visualizes wave functions for better understanding.

Findings

Successful translation of quantum systems into matrix form for quantum algorithms

Visualization of wave functions aligns with theoretical continuous operator results

Demonstrates the applicability of discrete quantum mechanics in quantum computing

Abstract

Quantum computers have the potential to transform the ways in which we tackle some important problems. The efforts by companies like Google, IBM and Microsoft to construct quantum computers have been making headlines for years. Equally important is the challenge of translating problems into a state that can be fed to these machines. Because quantum computers are in essence controllable quantum systems, the problems that most naturally map to them are those of quantum mechanics. Quantum chemistry has seen particular success in the form of the variational quantum eigensolver (VQE) algorithm, which is used to determine the ground state energy of molecular systems. The goal of our work has been to use the matrix formulation of quantum mechanics to translate other systems so that they can be run through this same algorithm. We describe two ways of accomplishing this using a position basis approach and a Gaussian basis approach. We also visualize the wave functions from the eigensolver and make comparisons to theoretical results obtained with continuous operators.

Equations

$$\ell (a)=\frac{2a-1-n}{2},a\in \{1,2,\dots ,n\}$$

$$⟨j|X_{pos}|k⟩=\sqrt{\frac{2\pi }{n}}\ell (j){\delta }_{j,k}$$

$$⟨j|F|k⟩=\frac{1}{\sqrt{n}}{e}^{\frac{2\pi i}{n}\ell (j,n)\ell (k,n)}$$

$$P_{pos}=F^{†}{X}_{pos}F$$

$$⟨j|A|k⟩=\sqrt{j}{\delta }_{j,k-1}$$

$$X_{en}=\frac{1}{\sqrt{2}}(A^{†}+A),P_{en}=\frac{i}{\sqrt{2}}(A^{†}-A)$$

$${\hat{H}}_{ho}=\frac{{\hat{p}}^2}{2m}+\frac{1}{2}m{\omega }^2{\hat{x}}^2$$

$$H_{ho}=\frac{P^2}{2}+\frac{X^2}{2}$$

$$H_{aho}=\frac{P^2}{2}+\frac{X^2}{2}-\alpha X^3$$

$$H=\frac{P^2}{2}+V(X)$$

$$A=\frac{1}{\sqrt{2}}(X+iP),A^{†}=\frac{1}{\sqrt{2}}(X-iP)$$

$$H_{ho}=A^{†}A+I/2$$

$$E(n)=\hbar\omega\bigg{(}n+\frac{1}{2}\bigg{)}=n+\frac{1}{2}$$

$$H_{ho}=\frac{1}{2}\big{(}A^{\dagger}A+AA^{\dagger}\big{)}$$

$$\left[a,a^{†}\right]=1$$

$$\left[A,A^{†}\right]=I$$

$$\left[A,A^{†}\right]=\left(\begin{array}{cccccc}1& 0& 0& \dots & 0& 0\\ 0& 1& 0& \dots & 0& 0\\ 0& 0& 1& \dots & 0& 0\\ ⋮& ⋮& ⋮& \ddots & ⋮& ⋮\\ 0& 0& 0& \dots & 1& 0\\ 0& 0& 0& \dots & 0& -(n-1)\end{array}\right)$$

$$H_{ho}^{\prime }=\frac{1}{2}(X^2+P^2-\tilde{I}+I)=A^{†}A+I/2$$

$$V_{cub}=X^2-\alpha X^3$$

$$V_{qua}=X^2+\beta X^4$$

$$W_{cub}(n)=\hbar\omega_{0}\bigg{(}n+\frac{1}{2}\bigg{)}-\frac{5\lambda^{2}\hbar^{2}}{12m\omega_{0}^{4}}\bigg{(}n^{2}+n+\frac{11}{30}\bigg{)}$$

$$W_{quar}(n)=\hbar\omega_{0}\bigg{(}n+\frac{1}{2}\bigg{)}+\frac{3\lambda\hbar^{2}}{8m\omega_{0}^{2}}\bigg{(}n^{2}+n+\frac{1}{2}\bigg{)}\\ -\frac{\lambda^{2}\hbar^{2}}{64m^{2}\omega_{0}^{5}}\bigg{(}17n^{3}+\frac{51}{2}n^{2}+\frac{59}{2}n+\frac{21}{2}\bigg{)}$$

$$C|0⟩=0,C|1⟩=|0⟩;$$

$$C^{†}|0⟩=|1⟩,C^{†}|1⟩=0$$

$$C=\left(\begin{array}{cc}0& 1\\ 0& 0\end{array}\right),C^{†}=\left(\begin{array}{cc}0& 0\\ 1& 0\end{array}\right)$$

$$|0⟩=\left(\begin{array}{c}1\\ 0\end{array}\right),|1⟩=\left(\begin{array}{c}0\\ 1\end{array}\right)$$

$${\tilde{O}}_B=O_B\otimes I_F$$

$${\tilde{O}}_F=I_B\otimes O_F$$

$$H=H_0+g(U_0+q\left[a^{†},a\right]);$$

$$\displaystyle H_{0}=\frac{1}{2}\big{(}p^{2}+\omega_{0}q^{2}+\omega_{0}\left[a^{\dagger},a\right]\big{)};$$

$$U_0={\omega }_0{q}^3+\frac{1}{2}g{q}^4$$

$$H=\omega_{0}\Big{(}A^{\dagger}A+\frac{1}{2}I\Big{)}\\ +\frac{1}{2}\Big{(}2\omega_{0}gX^{3}+g^{2}X^{4}+\big{(}\omega_{0}I+2gX\big{)}\left[C^{\dagger},C\right]\Big{)}$$

Full text

I Introduction

Once they are fully implemented, quantum computers will likely occupy a role similar to graphics cards in the overall computing landscape. Where graphics cards are very well suited to performing parallel computations, quantum computers can be used to perform calculations which span a vast amount of variables. The computational power of a quantum computer rises exponentially with its available qubits, outstripping supercomputers at a mere 50.

This project aims to leverage quantum computation to calculate the ground state energy of some Hamiltonians. In order to translate this program so that it can be fed to the quantum computer, we need to use the language of linear algebra. We will represent the operators of quantum mechanics with matrices and the states as vectors. To create and manipulate these objects I utilized Mathematica. For the quantum computation I utilized the QISKit Python library made by IBM. It allows us to load in the matrices we created and run quantum programs on them using either simulators or IBM’s own quantum computers over the Cloud.

II Operators in the Position Basis

In the position basis, we work on a lattice with $n$ sites indexed $\frac{1-n}{2},\frac{3-n}{2},\frac{5-n}{2},\dots,\frac{n-3}{2},\frac{n-1}{2}$. To make notation simpler, we introduce a mapping from the natural numbers to our lattice:

[TABLE]

The eigenvectors of our Hamiltonians will represent the probability density for a particle to be found on a particular site of the lattice.

We start by constructing the position operator as an $n\times n$ matrix:

[TABLE]

So the position operator is a traceless diagonal matrix made from our lattice indices with a scaling factor in front. To construct the momentum operator, we will first need a matrix for the discrete Fourier transform:

[TABLE]

Where $\omega=e^{\frac{2\pi i}{n}}$ is a primitive n-th root of unity. The exact form of the Fourier transform matrix depends on the lattice that we choose. The momentum operator is then

[TABLE]

We proceeded in a fashion similar to the approach described by Jagannathan and Vasudevan in [1]. In their paper they directly defined both their position and momentum operators and then detailed the relations between them and the finite Fourier transform. Instead we created the same position operator and Fourier transform and then used them to create a momentum operator that would satisfy all the same relations. We chose this path because it works for any size of matrix operator; the method described in the paper works only for odd values, and we need even values to implement these operators on a quantum computer, as we will explain later.

III Operators in the Harmonic Oscillator Energy Basis

In the energy basis, the eigenfunctions of our Hamiltonians will represent the probability amplitudes for our system to have a certain energy. We start this time by constructing a discrete version of the annihilation operator as the $n\times n$ matrix

[TABLE]

Then our position and momentum operators are

[TABLE]

where $i=\sqrt{\-1}$.

Additional work on discrete quantum mechanics can be found in papers by Meurice [2], Singh et. al. [3], Macridin et. al. [4], and Berenstein [5].

IV Building Hamiltonians

Once we have made our operators in the desired basis, constructing our Hamiltonians proceeds as usual. The Hamiltonian of the quantum harmonic oscillator is written as

[TABLE]

Setting the constants $\hbar$, $\omega$, and $m$ to 1 and replacing the usual continuous operators with our discrete ones, we get

[TABLE]

Similarly, the Hamiltonian for the anharmonic oscillator with a cubic term is given by

[TABLE]

In general, Hamiltonians can be expressed as the sum of a kinetic term and a potential:

[TABLE]

In certain cases it is advantageous to construct our Hamiltonians from the raising and lowering operators instead. In this case we can use the inverses of equations 6 to get these operators from the momentum and position operators:

[TABLE]

Then our harmonic oscillator Hamiltonian, for example, can be written as:

[TABLE]

Where $I$ is the $n\times n$ identity matrix.

V Comparison - Harmonic Oscillator

We can compare the position and energy bases by implementing a few simple Hamiltonians in each and looking the spectra of the matrices. We start with the simple harmonic oscillator formulated using the position and momentum operators as in equation 8.

The energy levels of the harmonic oscillator are given by

[TABLE]

We observe that the position basis spectrum matches the exact expression very well for the first $3/4$ of the energies, but then diverges. The energy basis spectrum matches the exact expression, but has a kink in the middle. This happens because the Hamiltonian matrix in the energy basis has a degenerate eigenvalue. We could proceed by simply throwing out this “junk” eigenvalue, but when we get to more complicated Hamiltonians it could be more difficult to pick out.

While playing around with the different bases, I noticed that if we instead express the Hamiltonian using the ladder operators as in equation 12, its spectrum matches the exact expression for any value of n. To investigate this, I tried substituting equations 6 into equation 8 and expanding to get

[TABLE]

The continuum ladder operators satisfy the commutation relation

[TABLE]

If we take this to mean that our finite ladder operators should satisfy

[TABLE]

then equation 14 would match equation 12 exactly. However, if we compute the commutator for the matrix ladder operators we find that it is

[TABLE]

This makes sense because matrix commutators are traceless. So what we have created here is a “traceless identity” matrix, which we’ll call $\tilde{I}$. As the size of this matrix increases, it resembles the identity operator more and more. To get the spectrum to match our exact expression, we can use the expression

[TABLE]

Interestingly, if we try to use the form of equation 12 with the position basis ladder operators, the spectrum will also have a (close to) duplicate eigenvalue in the middle.

So it seems like we can intuit a tentative rule of thumb: If we want to work with ladder operators, we should use the energy basis; if we want to work with position and momentum operators, we should instead use the position basis.

Additional work on the discrete quantum harmonic oscillator can be found in papers by Atakishiyev et. al. [6], Lorente [7], Barker et. al. [8], and Aunola [9].

We also tested some perturbative Hamiltonians, namely the cubic and quartic anharmonic oscillators. We can compare the spectra of these systems to perturbative expressions formulated by Heisenberg. [10] The potentials are:

[TABLE]

[TABLE]

To implement these in the energy basis as just discussed, we simply append the perturbative part of each potential, written in terms of the position operator, to the basic quantum harmonic oscillator Hamiltonian written in terms of the ladder operators.

The perturbative energy expressions are:

[TABLE]

[TABLE]

The following figures show a comparison between the energy basis, position basis, and perturbative expressions:

In the cubic case, the two different spectra and the perturbative expression agree to within 1% for over half the eigenvalues. In the quartic case, the agreement is a slightly worse 5% for those same values. The agreement improves as we decrease the coupling to the perturbative component.

VI Supersymmetric Oscillator

Now we present another consistency check with a more complicated potential, namely the supersymmetric anharmonic oscillator. For the normal harmonic oscillator we use bosonic operators. Since bosons can exist in an infinite number of excited states, their matrix operators can be extended indefinitely. On the other hand, fermions can only be in two states, which we’ll call $\ket{0}$ and $\ket{1}$. The lowering and raising operators, $C$ and $C^{\dagger}$ respectively, must obey the following relations:

[TABLE]

Whereas with bosons the only “illegal” operation is lowering the ground state, with fermions you can’t raise the single excited state either. These conditions are satisfied by the $2\times 2$ matrices:

[TABLE]

The state vectors for fermions are then simply:

[TABLE]

In the simplest supersymmetric system, we have a boson paired with a fermion. In the context of matrix operators, this pairing is achieved by extending our original operators using Kronecker products. [11] For a bosonic operator:

[TABLE]

And for a fermionic operator:

[TABLE]

Here $I_{b}$ and $I_{f}$ are the identity matrices that are the same size as our bosonic and fermionic spaces, respectively. From now on we will use $A$, $C$, etc. to refer to the new extended operators. We will work with extended energy basis ladder operators, and the other operators can be built from these as described before.

Yu Musin [12] uses the following Hamiltonian:

[TABLE]

Here $q$ and $p$ are the position and momentum operators, $\omega_{0}$ and $g$ are constants, and $a$ and $a^{\dagger}$ are actually the fermionic ladder operators. $H_{0}$ is the supersymmetric harmonic oscillator and the rest is the anharmonic perturbation. Expanding, rewriting in terms of our matrix operators and replacing the non-perturbative part with ladder operators:

[TABLE]

Musin also calculated a perturbative expression for the energy levels of this Hamiltonian, which depends on $n_{B}$ and $n_{F}$, respectively the energy levels of our boson and fermion111There was a typo in the original paper: the $g$ in the second term should have been $g^{2}$. Equation 33 reflects this change.:

[TABLE]

When we compare the spectrum of equation 32 with the corresponding values of equation 33, we see again that they match quite closely for the first half of the eigenvalues:

VII Quantum Computing

VII-A Brief Overview

The basic operational units of quantum computers are qubits. Qubits in turn are two-state quantum systems, be they atomic spins or photon polarizations. The quantum computer consists of the qubits, the connections between them, and the apparatus used to manipulate their states. In most cases the qubits are very sensitive to mechanical and thermal perturbations, which can cause them to decohere and lose their quantum properties. To prevent this from happening the whole machine is encased in a multistage dilution refrigerator to keep it cold and housed in a shock absorbent chamber to keep it very still. Even with all of these precautions, qubits still typically decohere in a matter of microseconds.

Calculations are performed on a quantum computer by acting on single or multiple qubits via quantum logic gates and thereby changing their states. Since these are quantum systems, they can exhibit superposition and entanglement, which are key to their computational power. When a qubit is in superposition it has a probability of being measured in either of its two states. Entanglement refers to a state involving two qubits which can’t be described independently. If we measure one of them, we automatically know the state of the other.

We can conceptualize a qubit as a complex two-component vector:

[TABLE]

Here $\alpha$ and $\beta$ represent the probability amplitudes for the qubit being in its two states, so the sum of their squares must be one. With this in hand, we can think of single-qubit quantum gates as $2\times 2$ unitary matrices. They must be unitary to conserve probability. This formalism also allows us to think of the state of a qubit as existing on the surface of a unit sphere called the Bloch sphere:

Using the two angles $\theta$ and $\phi$, we can express $\alpha$ and $\beta$ as

[TABLE]

So quantum logic gates simply move the point representing the state of the qubit along the Bloch sphere. Some single-qubit gates include rotations about the x-, y- and z-axes, the familiar Pauli matrices, and the Hadamard gate which can put a qubit into superposition. Arguably the most important 2-qubit gate is the controlled-not or “CNOT” gate, which is commonly used to entangle two neighboring qubits.

Quantum gates are used to build up quantum circuits, which take our set of qubits to a desired state. Typically a circuit will end with all the qubits being measured to get the result. Since we are dealing with probabilistic objects, we can’t get an “answer” from a single measurement; we need to run the circuit many times and histogram all the results to get a probability distribution. The peaks in the distribution tell us what the outcome of the quantum program.

VII-B The Variational Quantum Eigensolver (VQE) Algorithm

The VQE algorithm [13] tries to approximate the ground state energy of a given matrix Hamiltonian through a hybrid classical/quantum optimization process. We start by initializing our qubits with a trial state populated with guessed parameters. The quantum computer then evaluates the expectation value of this state with the given Hamiltonian. Using the returned value, a classical optimization algorithm varies the parameters. These two steps cycle until a user-set goal is set, either for number of steps or stability of the calculated energy.

In actuality, the algorithm is not run on the whole Hamiltonian at once, since in most cases it can’t be represented as a single layer of quantum gates. To be fed into the VQE algorithm the Hamiltonian matrix must be decomposed into a sum of tensor products of Pauli matrices. This can be done because the Pauli matrices and the identity matrix form a complete basis for complex $2x2$ matrices, and their complete set of tensor products does the same for larger matrices.

VII-C Scripts

My workflow included three major steps:

Create Hamiltonian matrix (Mathematica) 2. 2.

Decompose Hamiltonian matrix into sum of Pauli tensor products (Mathematica) 3. 3.

Run VQE algorithm (IBM QISKit in Python)

If we limit ourselves to a small quantity of qubits - no more than 10 or so - we can use Mathematica to build the operators and make Hamiltonians with them. For each different type of Hamiltonian I wrote a script that would build the matrix and export it to a file as an array. I then wrote a second script that would take the matrix file as an input and output a list of Pauli tensor products and their corresponding coefficients.

The second script works by generating a list of all possible Pauli matrix tensor product for a given number of qubits, then taking the trace of the product of each with the Hamiltonian. This trace, multiplied by a normalization factor, is the coefficient we want. Each Pauli tensor with a non-zero coefficient is then written to an output file. It is this file which is read by the VQE script.

The VQE script’s main inputs are listed below:

Hamiltonian matrix as list of Pauli tensors in file 2. 2.

variational form (trail circuit) defined in separate Python script 3. 3.

quantum depth (number of times to repeat circuit in a single cycle) 4. 4.

maximum number of optimization cycles 5. 5.

device to use, either local simulator or IBM’s devices via the Cloud

Once it runs, the VQE script outputs a convergence plot, a circuit diagram depicting the trial function, and a log file. The outputs of runs on the same type of Hamiltonian are automatically saved to their own directory, and the logs of different runs are appended to the same file. All the runs I’ve done so far have been on a local simulator instead of IBM’s real quantum devices. For example, a run of a 6-cubit anharmonic oscillator Hamiltonian with a simple trial function, quantum depth of 3 and just 100 optimization steps took about 45 minutes to finish and did not converge completely to the exact energy.

We saw from limited runs that using the local simulator will require high-performance computing (HPC) for large runs. The harmonic and anharmonic oscillator Hamiltonian runs converged correctly when using a very simple circuit. More complicated Hamiltonians, involving supersymmetry and/or finite temperature, will require additional study. Example outputs for the harmonic oscillator, showing a convergence to within 2% of the expected ground state energy, can be seen in figure 7 below.

For additional information on discrete quantum mechanics using quantum computers, see the recent paper by Klco et. al. [14].

VIII Conclusion

Though the results are still rather fresh, we’ve seen some success running simple Hamiltonians on quantum computers using simple trial functions. Moving to bigger and/or more complicated Hamiltonians presents various challenges. Due to the nature of the tensor product, a composite system with just a few simple parts could have a monstrously large matrix Hamiltonian. Making larger matrices will likely become too time-consuming for Mathematica, so we’ll have to migrate our scripts to a lower level programming language. This would still leave us with the issue of efficiently storing and accessing these matrices, which will quickly outgrow the space on a typical laptop. The greatest challenge in fully exploiting the VQE algorithm will be designing good trial functions, which often need to be tailored to specific Hamiltonians. Even if we find one that works well for a given Hamiltonian, its effectiveness will not necessarily persist as we scale up the matrix size. We have also been working on implementing a Hamiltonian that describe finite temperature effects on a supersymmetric system, based on work by Das et al [15], as well as a matrix effective model studied by Kashiwa et. al. [16]. The former proved difficult to implement because of the complexity of the theory and the fact that the size of the matrix Hamiltonian grows more rapidly than usual as we increase the size of the boson components.

IX Acknowledgements

This project was supported in part by the Brookhaven National Laboratory (BNL), Computational Science Initiative under the BNL Supplemental Undergraduate Research Program (SURP) . The author would like to thank Michael McGuigan for many stimulating discussions, as well as past and present CSI interns for their help setting up the Mathematica and Python code. We acknowledge use of the IBM Q for this work. The views expressed are those of the authors and do not reflect the official policy or position of IBM or the IBM Q team. Michael McGuigan is supported from DOE HEP Office of Science DE-SC0019139: Foundations of Quantum Computing for Gauge Theories and Quantum Gravity.