Producing Energy Eigenstates of the ${\rm H}_{2}$ Molecule by Classically Emulated Quantum Simulation

Quantum information theory has been investigated for more than half a century[1]. Among others the Shor’s factoring algorithm[2] has caused great motivation to realize quantum computers. In the last decade, realization of quantum computers of the universal type have developed so much[3] that several studies have been done[4] using quantum computers of the universal type. Especially, problems in quantum physics and quantum chemistry would be natural objects of study by quantum computers. In quantum physics and quantum chemistry it is an important task to identify energy eigenstates of quantum systems such as condensed matters, atoms and molecules. Variational Quantum Eigensolver(VQE)[5] is a famous algorithm to prepare a ground state of a quantum system that uses both classical computers and quantum computers. Variational Quantum Deflation(VQD)[6] is an application of VQE that prepare exited states of a quantum system.

The ${\rm H}_{2}$ molecule is the simplest molecule and has been studied in terms of quantum computers [7, 8, 9] The ${\rm H}_{2}$ molecule can be described by a two-qubits Hamiltonian and energy eigenstate of this system have been studied by quantum simulation and by NMR[10].

Recently, the author of the present paper has proposed a quantum algorithm that produces energy eigenstates of a Hamiltonian by classically emulated quantum simulation[11]. Our procedure uses a proper time evolution by a Hamiltonian controlled by concatenated ancilla qubits. Starting from an adequate initial state we can attain the corresponding energy eigenstate. Energy eigenstates have been produced for a few-qubits Hamiltonian that derives from the (1+1)-dimensional lattice Schwinger model[12]. The purpose of this paper is to apply the above quantum algorithm to the ${\rm H}_{2}$ molecule that is more symmetric than the model that derives from the (1+1)-dimensional Schwinger model. We produce energy eigenstates of the ${\rm H}_{2}$ molecule by classically emulated digital quantum simulation without using VQE nor VQD. We do not use algebraic calculation nor classical-quantum hybrid algorithm. Our algorithm uses ancilla qubits and resembles the quantum phase estimation algorithm [16]. Our algorithm also resembles the energy difference estimation algorithm[17], where a phase difference that is caused by an energy difference is used.

Our starting point is the following Hamiltonian of the ${\rm H}_{2}$ molecule written in atomic units[7, 8, 9]

$$H=a_{0}I_{0}I_{1}+a_{1}Z_{0}I_{1}+a_{2}I_{0}Z_{1}+a_{3}Z_{0}Z_{1}+a_{4}X_{0}X_{1},$$

(1)

where $\{a_{i}\}$ are constants depending on the distance of the two nuclei, and $I$, $X$ and $Z$ represent the identity operator and the Pauli matrices

$$\sigma_{x}=\left(\begin{array}[]{cc}0&1\\ 1&0\end{array}\right),\quad\sigma_{z}=\left(\begin{array}[]{cc}1&0\\ 0&-1\end{array}\right),$$

(2)

that operate on each electron specified by a 0 or 1 subscript in the molecular orbital. The Hamiltonian (1) has been obtained from the ordinary quantum mechanical Hamiltonian by the second quantization and representing fermionic variables by Pauli matrices[13, 14].

In the matrix form we have

$$H=\left(\begin{array}[]{cccc}a_{0}+a_{1}+a_{2}+a_{3}&0&0&a_{4}\\ 0&a_{0}+a_{1}-a_{2}-a_{3}&a_{4}&0\\ 0&a_{4}&a_{0}-a_{1}+a_{2}-a_{3}&0\\ a_{4}&0&0&a_{0}-a_{1}-a_{2}+a_{3}\end{array}\right).$$

(3)

We use the values $a_{0}=-1.04391,a_{1}=0.42045,a_{2}=-0.42405,a_{3}=-0.01150$ and $a_{4}=0.179005$ for the internuclear distance 0.70Å[10]. The energies of the ground state and exited states are algebraically calculated to be $E_{0}=-1.89215,E_{1}=-1.23442,E_{2}=-0.876405$ and $E_{3}=-0.172668[{\rm E_{h}}]$, where $1[{\rm E_{h}}]=27.211386[{\rm eV}]$. In the following, we assume that we do not know these values.

We explain twirling operations[11, 15] that produce a specific energy eigenstate $|E_{i}\rangle$ from a superposition of energy eigenstate $|\psi\rangle=\Sigma_{k}c_{k}|E_{k}\rangle$. We introduce ancilla qubits in the initial state $|0\rangle$, We Hadamard transform an ancilla qubit as $|0\rangle\rightarrow{1\over\sqrt{2}}(|0\rangle+|1\rangle)$. We apply time evolution of the Schrödinger equation for the Hamiltonian controlled by the ancilla qubit to the physical state $|\psi\rangle=\Sigma_{k}c_{k}|E_{k}\rangle$. By setting the time period adequately, only a specific energy eigenstate $|E_{i}\rangle$ can acquire the phase $2\pi$ by the time evolution. After the time evolution we again Hadamard transform the ancilla qubit. Only the state of ancilla qubit that associate with the specific energy eigenstate $|E_{i}\rangle$ can return to the original state

$$|E_{i}\rangle{|0\rangle+|1\rangle\over\sqrt{2}}\xrightarrow{H}|E_{i}\rangle|0\rangle.$$

(4)

Another energy eigenstate $|E_{k}\rangle$ acquire a non-trivial phase factor $e^{i\theta_{k}}(\neq 1)$ by the time evolution and the ancilla qubit that associate with the energy eigenstate $|E_{k}\rangle$ does not return to the initial state $|0\rangle$ by the Hadamard transformation

$$|E_{k}\rangle{|0\rangle+e^{i\theta_{k}}|1\rangle\over\sqrt{2}}\longarrownot\longrightarrow|E_{k}\rangle|0\rangle.$$

(5)

We concentrate on the physical state that is associated with the ancilla qubit that remains to be $|0\rangle$. We call the physical states that are associated with ancilla qubit $|0\rangle$ active states. After the twirling operation, the relative amplitude of the specific state $|E_{i}\rangle$ increases in the active state. Using the new physical state the next twirling operation is prescribed. We adopt the new time period $\tau_{1}={\pi\over 2E_{i1}}$ calculated by $E_{i1}=\langle\psi_{1}|H|\psi_{1}\rangle$, where $|\psi_{1}\rangle$ is the new active physical state. Introducing a new ancilla qubit, we can carry out the next twirling operation and the relative amplitude of the specific state $|E_{i}\rangle$ increases furthermore. Repeating this procedure, only the specific state $|E_{i}\rangle$ survives in the active state. Thus we can produce the specific state $|E_{i}\rangle$ from the initial physical state $|\psi\rangle=\Sigma_{k}c_{k}|E_{k}\rangle$.

We start from an initial physical state $|\psi_{i0}\rangle$, where the subscript $i$ means that $|\psi_{i0}\rangle$ can lead to an energy eigenstate $|E_{i}\rangle$ and the subscript $0$ means that no twirling operations have been done. We evaluate $E_{i0}=\langle\psi_{i0}|H|\psi_{i0}\rangle$. We choose the time period $\tau_{i0}={\pi\over 2E_{i0}}$ to use the identity $ie^{-{\pi{i}\over 2}}=1$. If we had chosen $\tau_{i0}={2\pi\over E_{i0}}$ we had used the identity $e^{-2\pi{i}}=1$. We carry out the first twirling operation to obtain the first approximate state $|\psi_{i1}\rangle$ of $|E_{i}\rangle$ in the active state. We evaluate $E_{i1}=\langle\psi_{i1}|H|\psi_{i1}\rangle$ and using $E_{i1}$ and a fresh ancilla qubit $|0\rangle$, we carry out the second twirling operation to obtain the second approximate state $|\psi_{i2}\rangle$ of $|E_{i}\rangle$ in the active state. Repeating this procedure the physical state $|\psi_{ij}\rangle$ will approach the state $|E_{i}\rangle$ as the number of the twirling operations $j$ increases. In Fig.1 we show the $j$-th block of the twirling operations that we use in simulation.

************************************** Fig.1 ********************************************

We show in Fig.2 a simulation result that we start from the initial physical state $|1\rangle_{0}|0\rangle_{1}$, where the subscripts 0 and 1 mean the $0$th qubit and the $1$st qubit, respectively. We see that the expectation values of $H$ cannot be distinguished from the theoretical value $E_{0}=-1.89215$ after the first twirling operation. After the first twirling operation, the expectation values of $H$ seem hardly change. This will mean that after the first twirling operation the physical state reaches a steady state and the ground state $|E_{0}\rangle$ has been produced.

********************************** Fig.2 ********************************

We show in Fig.3 a simulation result that we start from the initial phisical state $|-\rangle_{0}|+\rangle_{1}$, where $|\pm\rangle={|0\rangle\pm|1\rangle\over\sqrt{2}}$. At the initial point($j=0$) the expectation value ${}_{1}\langle+|_{0}\langle-|H|-\rangle_{0}|+\rangle_{1}$ takes a near value of $E_{1}$ by chance. The state $|-\rangle_{0}|+\rangle_{1}$, however, is not an eigenstate of $H$ and by the first twirling operation($j=1$) the expectation value of $H$ is separated from the value $E_{1}$. The initial physical state $|-\rangle_{0}|+\rangle_{1}$ is expanded as $|-\rangle_{0}|+\rangle_{1}=c_{0}|E_{0}\rangle+c_{1}|E_{1}\rangle+c_{2}|E_{2}\rangle+c_{3}|E_{3}\rangle$. By the first twirling operation, the relative amplitudes of $c_{0},c_{2}$ and $c_{3}$ in the active physical state diminish and the expectation value of $H$ falls. We see that the expectation values of H cannot be distinguished from the theoretical value $E_{1}=-1.23442$ after the number of the twirling operations $j=4$, which will mean that the first excited state $|E_{1}\rangle$ has been produced by the fourth twirling operation.

********************************** Fig.3 ********************************

We show in Fig.4 a simulation result that we start from the initial physical state $|+\rangle_{0}|-\rangle_{1}$. We see that after the number of the twirling operations $j=4$ the seond exited state $|E_{2}\rangle$ has been produced in the active state.

********************************** Fig.4 ********************************

We show in Fig.5 a simulation result that we start from the initial physical state $|0\rangle_{0}|1\rangle_{1}$. We see that after the number of the twirling operations $j=1$ the third exited state $|E_{3}\rangle$ has been produced in the active state.

********************************** Fig.5 ********************************

We show in Tbale I two typical simulation results that we start from the initial state $|0\rangle_{0}|0\rangle_{1}$. There are two cases that the expectation values of $H$ can approach to $E_{1}$ or $E_{2}$. The initial physical state $|0\rangle_{0}|0\rangle_{1}$ will be a point of a watershed.

Table I  Expectation values of $H$ for the initial physical state $|0\rangle_{0}|0\rangle_{1}$. The active state means the number of the active states that have survived among the $10^{7}$ initial states.

	$j=0$	$j=1$	$j=2$	$j=3$	$j=4$	$j=5$	$j=6$
$|0\rangle_{0}|0\rangle_{1}$	-1.05538	-1.05633	-0.94812	-0.90006	-0.89804	-0.87099	-0.87118
active state	$10^{7}$	6135059	108699	425	287	258	218
$|0\rangle_{0}|0\rangle_{1}$	-1.05547	-1.051434	-1.18508	-1.17334	-1.23233	-1.23379	-1.23400
active state	$10^{7}$	3989101	1637809	1233808	518612	401809	390440

We have studied energy eigenstates of the Hamiltonian of the ${\rm{H}_{2}}$ molecule. We have used concatenated ancilla qubits to carryout the twirling operations repeatedly. Starting from adequate initial physical states, we have produced energy eigenstates of the Hamiltonian of the ${\rm{H}_{2}}$ molecule without algebraic calculation. We have chosen five typical initial physical states. Any physical state $|\psi\rangle=\sum_{k}c_{k}|E_{k}\rangle$ will attain one of the four energy eigenstates by the twirling operations. The ground state and the highest energy eigenstate will be produced by a smaller number of the twirling operations. It requires a greater number of the twirling operations to produce the first and the second exited states. The Hamiltonian of the ${\rm H}_{2}$ molecule is symmetric, and in general there is no one-to-one correspondence between the initial physical state and the energy eigenstate of the end. The initial physical state $|0\rangle_{0}|0\rangle_{1}$ is one of points of a watershed leading to $|E_{1}\rangle$ or $|E_{2}\rangle$. Our procedure is applicable to any Hamiltonian that is represented by Pauli matrices.

The author thanks National Institute of Technology, Gunma College, where this study has been done.