文献情報

• タイトル：Quantum simulation of the Hubbard model on a graphene hexagon: Strengths of IQPE and noise constraints
• 著者：Mohammad Mirzakhani, Kyungsun Moon
• 書誌情報(DOI)：https://doi.org/10.48550/arXiv.2506.05031

概要

ハバード模型は電子同士の強い相互作用を記述する物理学的に重要なモデルですが、系のサイズが大きくなると古典コンピューターで解くことが指数関数的に困難になります。この基底エネルギーを推定する手法としてQPE(Quantum Phase Estimation / 量子位相推定)アルゴリズムが有望視されていますが、多くの量子ビットを要するため現状のNISQデバイス上で実行することは困難です。

本論文では、QPEを改良した量子アルゴリズムであるIQPE(Iterative Quantum Phase Estimation / 反復量子位相推定)アルゴリズムを利用し、グラフェンのハバード模型における基底エネルギー探索への有効性を、Qiskitを用いたシミュレーションにより確認しました。また、ノイズを加えたシミュレーションや量子コンピューター実機での検証も行い、IQPEのノイズ耐性や実機における実用性についても検証しました。

なお、本論文ではIQPEに関する検証のほかにも、断熱時間発展を利用した基底状態の物理量の計算を行っています。こちらの記事ではIQPEの検証部分に焦点を当てるため、物理量の計算に関する説明は割愛します。

背景

ここでは、ハバード模型と既存の量子アルゴリズムが抱える課題について解説します。

ハバード模型は、強相関電子系の複雑な振る舞いを表現するための代表的なモデルです。このモデルは磁性体、超伝導体、モット絶縁体転移[1]といった重要な現象を捉えており、物性物理学の探求において不可欠なものとなっています。しかしながら、ハバード模型は系のサイズが大きくなるにつれ解くことが指数関数的に困難になるという特徴をもっており、22サイトの格子上の17電子という小さな系が、現時点で厳密解が求められている最大のサイズ[2]です。

QPEアルゴリズムはこうした系の基底エネルギーを計算するための有望な量子アルゴリズムです。しかし、実行には多くの量子ビットを要するため、量子ビットが少なくエラーの多い現状のNISQデバイスでは実行が困難とされています。NISQデバイスでも実行可能な量子アルゴリズムとしてVQE(Variational Quantum Eigensolver/変分量子固有値ソルバー)[3,4]やQAOA((Quantum Approximate Optimization Algorithm/量子近似最適化アルゴリズム)[5]といった古典量子ハイブリッドアルゴリズムが存在しますが、これらは古典コンピューターにおけるエラー緩和の計算負荷が大きく、性能に限界があります。

手法

本論文では、QPEをNISQデバイスでも実行可能な形に改良したIQPEアルゴリズムのハバード模型に対する有用性を調べます。ハバード模型の対象として、図1のような形の6原子からなるグラフェン片を採用します。

グラフェンのハバード模型

ハミルトニアン

6原子グラフェンのハミルトニアン $\mathcal H$ を次の式で定義します：

$$
\mathcal H=\sum_{i,\sigma}\epsilon_{i\sigma}n_{i\sigma}-\gamma_0\sum_{\langle i,j\rangle,\sigma}(a^\dagger_{i\sigma}a_{j\sigma}+a_{j\sigma}^\dagger a_{i\sigma})+U_0\sum_in_{i\uparrow}n_{i\downarrow}\\n_{i\sigma}=a_{i\sigma}^\dagger a_{i\sigma}
$$

ここで、 $a_{i\sigma}^\dagger,a_{i\sigma}$ は生成演算子、消滅演算子と呼ばれる演算子で、それぞれ位置 $i(=0,\dots,5)$ 、スピン $\sigma(=\uparrow,\downarrow)$ に属する電子を一つ生成/消滅させることを意味します。 $n_{i\sigma}$ は個数演算子と呼ばれる演算子で、位置 $i$、スピン $\sigma$ の電子の個数を表します。

第一項はオンサイトエネルギー項です。係数 $\epsilon_{i\sigma}$は位置 $i$ 、スピン $\sigma$ の電子がもつエネルギーを表します。今回は係数を $\epsilon_{i\sigma}=0$ とするため、この項は消滅します。

第二項はホッピング項とよばれ、隣り合う原子間での電子の飛び移りにより生じるエネルギーを表します。 $\sum$ の $\braket{i,j}$ は隣り合う2点の組み合わせを表します。今回は簡単のため、係数を $\gamma_0=1$ とします。

第三項はハバード項、あるいは相互作用項とよばれ、クーロン相互作用による電子同士の反発によるエネルギーを記述します。 $\sum_in_{i\uparrow}n_{i\downarrow}$ という和の取り方から分かるように、このモデルでは同じ位置に2つの電子があるときのみ反発すると仮定しています。係数 $U_0$ は反発の強さを表します。

まとめると、今回扱うグラフェンのハバード模型のハミルトニアンは

$$
\mathcal H = \mathcal H_0+\mathcal H_\mathrm U
$$

のように、大きく分けて2つの項で表せることがわかります。ここで、 $\mathcal H_0$ はホッピング項(隣り合う原子間の電子の飛び移りの効果)を表し、 $\mathcal H_\mathrm U$ はハバード項(同じ位置に存在する電子間の反発の効果)を表します。

量子ビットとしての表し方

グラフェン上の電子は実際には6種類の位置 $i(=0,\dots,5)$ と2種類のスピン $\sigma=(\uparrow,\downarrow)$ で指定されます。しかし本論文では、12パターンの電子を量子ビット上で表すために、代わりにスピンのない12個の粒子を考えます。すなわち量子ビットは12個用いることになり、インデックスは単なる $i(=1,\dots,12)$ のみで指定することになります。

IQPE

IQPE(Iterative Quantum Phase Estimation / 反復量子位相推定)は、QPE(Quantum Phase Estimation / 量子位相推定)を量子ビットの少ないNISQデバイスでも実行しやすい形に改良したものです。

まず、通常のQPEについて説明します。QPEはユニタリ演算子の固有値を求めるためのアルゴリズムです。ユニタリ演算子の固有値は大きさ1の複素数であることから、位相 $\phi\,(0\le\phi\le1)$ を用いて $e^{i2\pi\phi}$ と表すことができます。ここで、この位相 $\phi$ は二進数表記で $\phi=0.\phi_0\phi_1\dots\phi_{n-1\, (2)}(\phi_i=0,1)$ と表すことができます。( $n$ は精度の桁数であり、大きくするほど正確になります。) QPEは、状態 $\ket{\phi_0}\otimes\ket{\phi_0}\otimes\cdots\otimes\ket{\phi_{n-1}}$ を生成し、測定することで $\phi$ の値を求めるという仕組みのアルゴリズムです。例えば、測定して得られた状態が $\ket1\otimes\ket0\otimes\ket1$ だったとします。このとき位相は二進数表記で $\phi=0.101_{(2)}$であるため、十進数に直すことで位相が $\phi=5/8$ であることがわかります。

QPEは、エネルギー固有値を求めるのに利用することができます。系のハミルトニアンが $\mathcal H$ であるとき、時間発展演算子は $e^{-i\mathcal Ht}$であり、これはユニタリ演算子です。 $e^{-i\mathcal Ht}$ の固有値はエネルギー固有値 $E$ を用いて $e^{-iEt}$と表すことができるので、 $E=-2\pi\phi$ より $\phi$ から $E$ を求めることができます。本論文では、この考え方を利用してグラフェンの基底エネルギーを求めることになります。

IQPEの基本的な考え方はQPEと同じですが、その実装方法に違いがあります。通常のQPEは、状態 $\ket{\phi_0}\otimes\ket{\phi_0}\otimes\cdots\otimes\ket{\phi_{n-1}}$を作成する必要があるため、そのために補助ビットを桁数 $n$ の分だけ用意しなければなりません。したがって、精度を大きくするほど補助ビットが多くなり、NISQデバイスでの実装が困難になります。一方、IQPEでは反復的な操作を取り入れることにより、たった1つの補助ビットで位相を求めることができます。

IQPEの量子回路

以下に具体的なIQPEの手順を記します。

まず基底状態 $\ket{\psi_0}$と補助ビット $\ket+=\frac{1}{\sqrt2}(\ket0+\ket1)$を用意し、次の操作を $k=1$から $k=m$まで繰り返し行います。

1. 制御ユニタリゲート $\mathcal U^{2^{m-k}}$ の作用
   時間発展演算子を $\mathcal U=e^{-i\mathcal Ht}$ として、その $2^{m-k}$ 乗を、補助ビット $\ket+$ を制御ビットとした制御ゲートとして作用させます。すなわち、補助ビットが $\ket1$ の場合のみ $\ket{\psi_0}$ に $\mathcal U$ を $2^{m-k}$ 回作用させる操作を行います。 $\ket{\psi_0}$ は固有状態であることから、 $\mathcal U\ket{\psi_0}=e^{-iEt}\ket{\psi_0}=e^{i2\pi \phi}\ket{\psi_0}$ より、系全体の状態は次のようになります： $$
   \begin{align*}
   \ket+\ket{\psi_0}
   &=\frac1{\sqrt2}(\ket0+\ket1)\ket{\psi_0}\\
   &\to\frac1{\sqrt2}(\ket0+e^{i2\pi2^{m-k}\phi}\ket1)\ket{\psi_0}.
   \end{align*}
   $$ ここで、$\phi$ の二進数表記 $\phi=0.\phi_1\phi_2\dots\phi_{m\, (2)}$を利用すると $$
   \begin{align*}
   2^{m-k}\phi&=2^{m-k}\times0.\phi_1\phi_2\dots\phi_m\\
   &=\phi_1\phi_2\dots\phi_{m-k}+0.\phi_{m-k+1}+0.0\phi_{m-k+2}\dots\phi_{m}
   \end{align*}
   $$ となります。第一項 $\phi_1\phi_2\dots\phi_{m-k}$ は整数部分であることから $e^{i2\pi\phi_1\phi_2\dots\phi_{m-k}}=1$より消滅して、状態は $$
   \begin{align*}
   &\quad\frac1{\sqrt2}(\ket0+e^{i2\pi2^{m-k}\phi}\ket1)\ket{\psi_0}\\
   &=\frac1{\sqrt2}(\ket0+e^{i2\pi(0.\phi_{m-k+1}+0.0\phi_{m-k+2}\dots\phi_{m})}\ket1)\ket{\psi_0}\\
   \end{align*}
   $$ となります。
2. $Z$軸回転ゲートの作用 補助ビットに $Z$軸回転ゲート $P(-\varphi_k)$ を作用させます。ここで、$\varphi_k=0.0x_{m-k+2}\dots x_{m-1}x_m$ は $k-1$ 回目の試行までで既に求められた値であることに注意してください。この作用により $e^{i2\pi(0.\phi_{m-k+1}+0.0\phi_{m-k+2}\dots\phi_{m})}$ の $e^{i2\pi0.0\phi_{m-k+2}\dots\phi_{m}}$ の部分は打ち消されて、状態は $$
   \begin{align*}
   &\quad\frac1{\sqrt2}(\ket0+e^{i2\pi(0.\phi_{m-k+1}+0.0\phi_{m-k+2}\dots\phi_{m})}\ket1)\ket{\psi_0}\\
   &\to\frac1{\sqrt2}(\ket0+e^{i2\pi0.\phi_{m-k+1}}\ket1)\ket{\psi_0}
   \end{align*}
   $$ と変化します。
3. $H$ゲートの作用 補助ビットにアダマールゲート $H$ を作用させます。状態は $$
   \begin{align*}&\quad\frac1{\sqrt2}(\ket0+e^{i2\pi0.\phi_{m-k+1}}\ket1)\ket{\psi_0}\\&\to\frac12\big((\ket0+\ket1)+e^{i2\pi0.\phi_{m-k+1}}(\ket0-\ket1)\big)\ket{\psi_0}\\
   &=\frac12\big((1+e^{i2\pi0.\phi_{m-k+1}})\ket0+(1-e^{i2\pi0.\phi_{m-k+1}})\ket1\big)\ket{\psi_0}
   \end{align*}
   $$ と変化します。ここで、 $$
   e^{i2\pi0.\phi_{m-k+1}}=e^{i\pi\phi_{m-k+1}}=(-1)^{\phi_{m-k+1}}
   $$ であるから $$
   =\frac12\big((1+(-1)^{\phi_{m-k+1}})\ket0+(1-(-1)^{\phi_{m-k+1}})\ket1\big)\ket{\psi_0}
   $$ となります。
4. 補助ビットの測定 補助ビットを測定し、その結果が $\phi$ の第 $m-k+1$ 桁の値を求めます。測定結果がもし $\ket 0$ なら $\phi_{m-k+1}=0$ であり、 $\ket 1$ なら $\phi_{m-k+1}=1$ であることがわかります。
5. 反復 $k$を1増やして手順 2. に戻ります。 $k$ が精度の桁数 $m$ に達したらIQPEは終了です。

以上の方法により、位相 $\phi$ を最も小さい桁から順に決定していくことができます。

量子ゲートの実装

IQPEでは、時間発展演算子 $\mathcal U =e^{-i\mathcal H t}$ を利用しました。しかし、この式を量子ゲートとして量子回路上に実装するためにはいくつかの課題があります。一つは、ハミルトニアン $\mathcal H$ が $X,Y,Z$ などの量子ゲートではなく生成消滅演算子 $a, a^\dagger$ で表されていることです。もう一つは、演算子が指数関数の肩に乗っていることです。前者を解決するのがジョルダン・ウィグナー変換であり、後者を解決するのが鈴木・トロッター分解です。

ジョルダン・ウィグナー変換

ジョルダン・ウィグナー変換[6]は、生成消滅演算子 $a, a^\dagger$ をパウリ演算子 $X,Y,Z,I$ の組み合わせに置き換える変換です。この変換により、 $\mathcal H$ のような生成消滅演算子で表された演算子を量子回路上で実行可能な形に書き換えることができます。ここではジョルダン・ウィグナー変換の一般論に触れることはせず、ハミルトニアンの各項がどのように変化するかについてのみ説明します。

ホッピング項 $h_0=-\gamma_0(a_{i\sigma}^\dagger a_{j\sigma}+a_{j\sigma}^\dagger a_{i\sigma})$ は次のように変換されます：

$$
(a_i^\dagger a_j+a_j^\dagger a_i)\mapsto\frac 1 2(X_iX_j+Y_iY_j)Z_\mathrm{JW}\quad(j>i)
$$

ここで、 $Z_\mathrm{JW}=\otimes_{k=i+1}^{k=j-1}Z_k$ はJWストリングと呼ばれる演算子で、電子の反交換関係を保つために現れます。JWストリングが現れるのは1番目と6番目、7番目と12番目のように実際は隣接しているにもかかわらず番号が離れている場合で、1と2, 2と3のように番号が隣接している場合は単に

$$
(a_i^\dagger a_j+a_j^\dagger a_i)\mapsto\frac 1 2(X_iX_j+Y_iY_j)\quad(j>i)
$$

となります。

相互作用項 $h_U=U_0n_{i\sigma}n_{i\sigma’}$は次のように変換されます：

$$
n_i=a_i^\dagger a_i\mapsto\frac 1 2(I_i-Z_i),\\
n_in_j\mapsto\frac 1 4 (I_iI_j-I_iZ_j-Z_iI_j+Z_iZ_j)
$$

以上を踏まえると、ハミルトニアンは

$$
\begin{align*}\mathcal H&=-\gamma_0\sum_{\langle i,j\rangle,\sigma}(a^\dagger_{i\sigma}a_{j\sigma}+a_{j\sigma}^\dagger a_{i\sigma})+U_0\sum_in_{i\uparrow}n_{i\downarrow}\\
&\mapsto-\frac{\gamma_0}2\sum_{\langle i,j\rangle}(X_iX_j+Y_iY_j)+\frac{U_0}4\sum_{i,j}(I_iI_j-I_iZ_j-Z_iI_j+Z_iZ_j)
\end{align*}
$$

のように変換できることがわかります。

JWストリングの対処

JWストリングとは、ホッピング項に現れる量子ゲート $Z_\mathrm{JW}=\otimes_{k=i+1}^{k=j-1}Z_k$ のことです。JWストリングが現れるのは1番目と6番目、7番目と12番目のように実際は隣接しているにもかかわらず番号が離れている場合で、本来2量子ビットゲートでよい操作が6量子ビットに作用することになり、計算コストの増大につながります。本論文では、JWストリングを単なるグローバル位相 $(-1)^{N_f-1}$( $N_f$ は電子の数)として表せることが示されており、ホッピング項は

$$
(a_i^\dagger a_j+a_j^\dagger a_i)\mapsto\frac {(-1)^{N_f-1}} 2(X_iX_j+Y_iY_j)\quad(j>i)
$$

となります。導出はAppendixで行います。

鈴木・トロッター分解

ジョルダン・ウィグナー変換によりハミルトニアン $\mathcal H$ を量子ゲートで表すことに成功しました。ここで、変換後のハミルトニアンが量子ゲートの和で表せることを踏まえると、時間発展演算子 $\mathcal U = e^{-i\mathcal H t}$ を適当な行列 $A_1,\dots,A_n$ を用いて

$$
\mathcal U=e^{A_1+\cdots+A_n}
$$

と表せます。ここで、一つ問題が生じます。 $A_1$ のような単純な行列の指数関数単体 $e^{A_1}$ は、量子回路上に図2のような形で実装できることが知られています。しかし、 $A_1+\cdots+A_n$ のような一般の行列に対しては不可能です。ここで、 $e^{A_1+\cdots+A_n}$ を指数法則を用いて $e^{A_1}\cdots e^{A_n}$ と分解できればよいと考えられるかもしれませんが、行列の非可換性より一般に $e^{A_1+A_2}\neq e^{A_1}e^{A_2}$ となります。つまり、行列の和を積の形にばらして計算を行うことはできません。この問題を解決するのが鈴木・トロッター分解[7]です。

鈴木・トロッター分解は、自然数 $n$ を十分大きくとるとき、

$$
e^{A_1+A_2}\approx(e^{A_1/n}e^{A_2/n})^n=\underbrace{e^{A_1/n}e^{A_2/n}\cdots e^{A_1/n}e^{A_2/n}}_{n回繰り返し}
$$

が近似的に成立することを主張しています。つまり今回の場合、 $\mathcal H=\mathcal H_1+\mathcal H_2\cdots+\mathcal H_k$ と表すならば、

$$
e^{-i\mathcal H t}\approx(e^{-i\mathcal H_1 \Delta t}e^{-i\mathcal H_2 \Delta t}\cdots e^{-i\mathcal H_k \Delta t})^{N_\mathrm{trot}}\quad(\Delta t=t/N_\mathrm{trot})
$$

とすることで、量子回路上での時間発展演算子 $\mathcal U$ の実装が可能になります。ここで、 $N_\text{trot}$ はトロッター分解数であり、大きくするほど近似の精度が向上します。

初期状態の準備

IQPEは、与えられたユニタリ演算子の固有状態に対し、対応する固有値を求めるアルゴリズムです。今回の論文で求めたいのは基底エネルギーであるため、正確な基底状態を用意する必要があるように思われます。しかし、ごく単純な系を除き、基底状態を求めることは不可能です。ですが、実用上は正確な基底状態を求める必要はなく、基底状態に十分近い状態が用意できればよいことがわかっています。その理由を解説します。

一般の状態を $\ket\chi$ とします。ユニタリ演算子の固有状態を $\ket{\psi_0}, \ket{\psi_1},\ket{\psi_2},\dots$ とし、特に基底状態を $\ket {\psi_0}$ とします。状態 $\ket\chi$ は固有状態の線形和で表せるため、

$$
\ket\chi=c_0\ket{\psi_0}+c_1\ket{\psi_1}+\cdots
$$

と書くことができます。状態 $\ket\chi$ にIQPEを適用することを考えます。まず、各固有状態 $\ket {\psi_i}$ に対してIQPEを適用すると、補助ビット $\ket0^n$ を含め全体の状態は

$$
\ket{0}^n\ket{\psi_i}\longrightarrow_\text{IQPE}\ket{\phi_{i,0}}\ket{\phi_{i,1}}\cdots\ket{\phi_{i,n-1}}\ket{\psi_i}
$$

と変化します。IQPEはユニタリ演算子であるため、線形演算子です。ゆえに、固有状態の線形結合で表せる $\ket\chi$ にIQPEを適用すると、

$$
\ket{0}^n\ket\chi=\ket{0}^n\sum_ic_i\ket{\psi_i}\longrightarrow_\text{IQPE}\sum_ic_i\ket{\phi_{i,0}}\ket{\phi_{i,1}}\cdots\ket{\phi_{i,n-1}}\ket{\psi_i}
$$

となります。この状態を観測すると、確率 $|c_i|^2$ で固有状態 $\ket{\psi_i}$ の固有値の位相 $\phi_{i,0}\phi_{i,1}\cdots\phi_{i,n-1}$ が求まります。したがって、 $|c_0|$ が十分に大きい状態が $\ket\chi$ であれば、十分大きい確率で基底エネルギーを求めることができます。すなわち、基底状態に十分近い状態を用意して何度もIQPEと観測を実行したのち、得られた固有値の中でもっとも頻出したものを採用すれば、それが基底エネルギーとなります。

では、基底状態に十分近い状態とは具体的にどのようなものでしょうか。本論文ではスレーター行列式[8]を用いています。スレーター行列式とは、相互作用が無いフェルミ粒子系の基底状態です。今回のモデルである6原子グラフェンでは相互作用を考慮するため、この近似が正しく機能するかどうかはわかりませんが、もっとも単純な近似であるためこれを採用します。スレーター行列式はQiskit Natureに実装されており、量子コンピューター上で簡単に準備することができます。

実験

グラフェンの基底エネルギーをIQPEと厳密対角化によって求め、IQPEが基底エネルギーの探索において有効であることを検証します。まずは、量子コンピューターを古典コンピューター上で再現するQiskit Natureライブラリ[9]を用いて、ノイズのない理想的な環境下におけるシミュレーションを行います。次に、ノイズを加えたシミュレーションを行います。最後に、量子コンピューター実機上で計算を行います。

ノイズなしシミュレーション

まずは、相互作用のない場合 $U_0=0$ と相互作用がある場合 $U_0=3$ のそれぞれで、基底エネルギーを全ての電子数 $N_\text{occ}=1,\dots,11$ について計算します。トロッター分解は $N_\text{trot}=15$ で行います。IQPEはQiskit Aerシミュレーター、厳密解はQuSpinパッケージを用いて計算を行います。結果は図3のようになります。いずれの場合もIQPEの結果が厳密解と非常によく一致していることがわかります。相互作用が無い場合は $N_\text{occ}=6$ で基底エネルギーが最小、相互作用がある場合は $N_\text{occ}=4$ で基底エネルギーが最小となることがわかります。

次に、IQPEの精度の桁数 $m$ とトロッター分解数 $N_\text{trot}$ が精度に及ぼす影響を検証します。(a)では $N_\text{trot}=15$ のもとで、 $m$ を動かします。(b)では $m=5$ のもとで、 $N_\text{trot}$ を動かします。どちらも相互作用は $U_0=3$ で行います。結果は図4のようになります。大きい $m,N_\mathrm{trot}$ でエネルギーが厳密値に収束することがわかります。

次に、様々な $U_0$ で基底エネルギーを計算します。$3\le N_\mathrm{occ}\le9$ でそれぞれ検証を行います。結果は図5のようになります。いずれの場合もIQPEの結果は厳密値と一致しています。全体的に、 $U_0$ の増加に伴い基底エネルギーも増加することがわかります。

全体として、ノイズなしのシミュレーションにおいては、 $m, N_\text{trot}$ が十分大きいときIQPEの結果と厳密値が非常によく一致することがわかります。

ノイズありシミュレーション

ノイズを加えることで、より実機に近づけたシミュレーションを行います。実機ibm_strasbourgの特性に合わせ、脱分極エラー、熱緩和エラー、読み出しエラーの三種類のエラーを再現します。なお、ノイズの導入により計算量が増大する問題を回避するため、系を6原子から3原子に縮小して行います。

まずは脱分極エラーの検証を行います。脱分極エラーとは量子ゲート実行中に発生するランダムな乱れであり、量子ビットの状態が完全混合状態に置き換わるエラーのことを指します。エラー発生確率は1量子ビットゲートと2量子ビットゲートで大きく異なり、当モデルではエラー確率の基準として実機ibm_strasbourgにおけるエラー確率の中央値を採用しています。1量子ビットゲートに対しては $p_1=2.230\times10^{-4}$、2量子ビットゲートに対しては $p_2=7.986\times10^{-3}$ と設定されています。実際の検証では、基準となる確率 $p_1,p_2$ の $1/10,1/5,1/2,1,2$ 倍でそれぞれ基底エネルギーの計算を行い、エラー確率の変化が精度に及ぼす影響を調べます。

結果は図6の通りです。青が1量子ビットゲートの結果、赤が2量子ビットゲートの結果を表します。一貫して2量子ビットゲートのエラーが及ぼす影響が1量子ビットゲートよりも大きいことがわかります。また、1量子ビットゲートはエラー確率を下げていけば推定値が厳密値に近づき、ばらつきも小さくなっていくのに対し、2量子ビットゲートはエラー確率を下げてもエネルギーを高く見積もったままで、ばらつきも大きいことがわかります。

次に、熱緩和エラーの検証を行います。熱緩和エラーとは環境との相互作用により情報が失われるエラーであり、脱分極エラーと異なりゲート操作を行っていない待機時間にも生じるという特徴を持っています。熱緩和エラーは2種類あり、一つはエネルギー緩和と呼ばれる、励起状態 $\ket 1$ が基底状態 $\ket 0$ に変化するエラーです。エネルギー緩和にかかる時間 $T_1$ の基準は実機に合わせて約 $300\, \mu \text s$ と設定します。もう一つは位相緩和と呼ばれる、重ね合わせの位相差情報が失われるエラーです。位相緩和にかかる時間 $T_2$ の基準は約 $160\,\mu\text s$ と設定します。実際の検証では、実機の $1/5,1/2,1,2,5,10,50$ 倍について検証することで、熱緩和時間とIQPE精度の関係性を調べます。また、量子ゲート実行時間の基準は1量子ビットと2量子ビットの場合についてそれぞれ $60 \text{ ns},660 \text{ ns}$ とし、その $1/6,1/2,1$ 倍について検証することで、量子ゲート実行時間とIQPE精度の関係性についても調べます。

結果は図7の通りです。(b)は横軸が熱緩和時間、(c)は横軸がゲート実行時間となっています。縦軸は図6と同様に基底エネルギーの推定値です。まず(b)について説明します。(b)は、熱緩和時間を変化させていったときのIQPE推定値の変化を、青(ゲート実行時間が実機の $1/6$)、赤(実機の $1/2$)、緑(実機の値)で比較した結果を表しています。図を見てみると、熱緩和時間が短いとき、ゲート実行時間をどの値に設定しても(すなわち、青でも赤でも緑でも)基底エネルギーを厳密値より大きく見積もってしまっていることがわかります。一方で、ゲート実行時間を実機の $1/6$ 倍に設定した青線に目を向けると、熱緩和時間が長くなるにつれ基底エネルギー推定値が厳密値に収束していくこともわかります。熱緩和時間を実機の $50$ 倍という極端に理想的な条件に設定すれば、どのゲート実行時間でも(すなわち、赤線や緑線の場合でも)厳密値に収束することがわかります。

(c)は、ゲート実行時間を変化させていったときのIQPE推定値の変化を、丸(ゲート実行時間が1秒=実質的にエラーのない状態)、三角(実機の10倍の熱緩和時間)、四角(実機相当の熱緩和時間)で比較した結果です。青系の線は1量子ビットエラーを、赤系の線は2量子ビットエラーを表します。まず、熱緩和時間が極端に長い丸の場合は、厳密値によく一致しています。三角の場合においても、1量子ビットエラーでは厳密値によく一致していますが、2量子ビットエラーは途中から外れてしまうことがわかります。四角の場合は、1量子ビットでも上昇傾向が見られ、2量子ビットでは最初から大きく外れていることがわかります。

もう一つのエラーである読み出しエラーは、測定の際 $\ket0$ を $\ket 1$ に、 $\ket 1$ を $\ket 0$ に読み間違えるエラーです。生じればIQPEの結果を一桁反転させることになりますが、今回の検証では結果に大きな影響を及ぼさなかったことが本文中で述べられています。

最後に、脱分極エラーと熱緩和エラーを混ぜた場合についてこれまでと同様に検証します。結果は図8の通りです。全体として、単独のエラーのみを発生させた場合よりも、推定値が正確になっていることがわかります。一方で、データのばらつき(エラーバー)は依然大きいままであり、結果に不確実性が残ることがわかります。

全体として、脱分極エラー確率と熱緩和時間の両方を最適化することが重要であり、特に2量子ビットエラーの影響が顕著であることがわかりました。

実機での検証

二種類の量子コンピューター実機ibm_strasburgとibm_fezを用意し、 $U_0$ を動かしたときのIQPE推定値の精度の変化を調べます。IQPEの桁数 $m$ とトロッター分解数 $N_\text{trot}$ は、ibm_strasburgでは $m=4,N_\text{trot}=12$、より高性能なibm_fezでは $m=5, N_\text{trot}=15$ と設定します。それぞれ50000ショットを20回試行し、その平均値とばらつきをプロットしたものが図9です。図9には参考として、ibm_strauburgに基づいたノイズありシミュレーションの結果も掲載しています。結果としては、ibm_strauburg系は一貫した基底エネルギーを過大評価している一方で、 $m,N_\text{trot}$ を大きく設定したibm_fezでは推定値が厳密値とよく一致していることがわかります。また、全体的に $U_0$ が大きくなるにつれてエラーバーが縮小しており、相互作用が強い領域においてはエラーの影響が緩和されることがわかります。

論文のまとめ

本論文では、6原子グラフェンのハバードモデルの基底状態をIQPEによって調査しました。回路の深さを削減するため、JWストリングの簡略化を行いました。ノイズなしのシミュレーションにおいては、IQPEを用いることによって初期状態が単なるスレーター行列式でも厳密値に非常によく一致した基底エネルギーを推定できることを示しました。ノイズありシミュレーションでは、熱緩和エラーや2量子ビットゲートエラーが特に有害であり、ノイズを混合した場合は精度が向上することを示しました。量子コンピューター実機を用いた検証ではは、ibm_fezが厳密値によく一致する結果を出すことがわかりました。

再現実装

本論文に基づき、Qiskit を用いてIQPEの再現実装を行いました。まずは図3で示したものと同様に、ノイズなしの環境下において各電子数 $N_\text{occ}$ での基底エネルギーを推定しました。計算量の問題から、系を6原子から3原子に縮小して実行しました。結果は図10の通りです。相互作用がある場合もない場合も、全体的にIQPE推定値が厳密値とよく一致することが確かめられました。

また、2量子ビットゲートであるCNOTゲートのエラー率の変化が推定値に及ぼす影響も調査しました。ノイズありの環境下では結果のばらつきが大きくなるため、各エラー率について何度も計算を行ったのちその平均値とエラーバーを表示するのが理想的ですが、計算時間の都合上それぞれ一回のみ計算を行いました。結果は図11の通りです。エラー率の小さいうちは比較的厳密値に近い結果が得られましたが、エラー率を大きくするにつれ結果が大きく外れることがわかりました。

Appendix

JWストリング簡略化の証明

$(i,j)=(1,6),(7,12)$ のとき、JWストリングの簡略化の式である

$$
(a_i^\dagger a_j+a_j^\dagger a_i)\mapsto\frac {(-1)^{N_f-1}} 2(X_iX_j+Y_iY_j)
$$

が成立することを証明します。

$(i,j)=(1,6)$ の場合から証明します。まず、生成消滅演算子 $a_j,a_j^\dagger$ に対しジョルダン・ウィグナー変換を行った結果は、スピン昇降演算子 $\sigma_j=\frac12(X_j+iY_j),\sigma_j^\dagger=\frac12(X_j-iY_j)$ を用いて

$$
\begin{align*}a_j\mapsto(\otimes_{k=1}^{j-1}Z_k)\otimes\sigma_j\\
a_j^\dagger\mapsto(\otimes_{k=1}^{j-1}Z_k)\otimes\sigma_j^\dagger\end{align*}
$$

と表せます。この変換規則を用いると、 $a_6^\dagger a_1$ は

$$
\begin{align*}a_6^\dagger a_1&\mapsto(\otimes_{k=1}^5Z_k)\otimes\sigma_6^\dagger \sigma_1\\
&=(Z_1\sigma_1)\otimes(\otimes_{k=2}^5Z_k)\otimes\sigma_6^\dagger
\end{align*}
$$

と変形でき、

$$
Z_1\sigma_1=\begin{pmatrix}1&0\\0&-1\end{pmatrix}\begin{pmatrix}0&1\\0&0\end{pmatrix}=\begin{pmatrix}0&1\\0&0\end{pmatrix}=\sigma_1
$$

より

$$
\begin{align*}a_6^\dagger a_1&\mapsto\sigma_1\otimes(\otimes_{k=2}^5Z_k)\otimes\sigma_6^\dagger
\end{align*}
$$

となります。ここで、個数演算子の変換式 $n_k\mapsto\frac12(I_k-Z_k)$ を等式とみなし変形すると $Z_k=I_k-2n_k$ が得られます。これを代入すると、

$$
\begin{align*}a_6^\dagger a_1&\mapsto\sigma_1\otimes(\otimes_{k=2}^5(I_k-2n_k))\otimes\sigma_6^\dagger
\end{align*}
$$

となります。個数 $n_k$ は $0,1$ のみの値をとることを踏まえると、 $n_k=0$ の場合は $I_k-2n_k=I_k$、 $n_k=1$ の場合は $I_k-2n_k=I_k-2I_k=-I_k$ と変形できるため、

$$
\otimes_{k=2}^5(I_k-2n_k)=(-1)^{\sum_{k=2}^5n_k}\otimes_{k=2}^5I_k
$$

が得られます。1番目から6番目の粒子の個数を $N_f=\sum_{k=1}^6n_k$ と定義して $\sum_{k=2}^5n_k=N_f-n_1-n_6$ を用いると

$$
\begin{align*}a_6^\dagger a_1&\mapsto(-1)^{N_f-n_1-n_6}\sigma_1\otimes(\otimes_{k=2}^5I_k)\otimes\sigma_6^\dagger\\
&=(-1)^{N_f-n_1-n_6}\sigma_1\otimes\sigma_6^\dagger
&\end{align*}
$$

となります。ここで、ホッピング項の物理的意味を考えます。電子が一方からもう一方に飛び移るということは、隣り合う1番目と6番目のうち一方が占有されており、もう一方は占有されていないということになります。したがって指数は $N_f-n_1-n_6=N_f-1$ となり、

$$
\begin{align*}a_6^\dagger a_1&\mapsto(-1)^{N_f-1}\sigma_1\otimes\sigma_6^\dagger
&\end{align*}
$$

が得られます。更にエルミート共役を取った形を加えることで、

$$
\begin{align*}a_1^\dagger a_6+a_6^\dagger a_1&\mapsto(-1)^{N_f-1}(\sigma_1^\dagger\otimes\sigma_6+\sigma_1\otimes\sigma_6^\dagger)
&\end{align*}
$$

が得られます。ここで $\sigma_j=\frac12(X_j+iY_j),\sigma_j^\dagger=\frac12(X_j-iY_j)$ を適用すると、

$$
\begin{align*}
\sigma_1^\dagger\otimes\sigma_6+\sigma_1\otimes\sigma_6^\dagger
&=\frac12(X_1-iY_1)\otimes\frac12(X_6+iY_6)+\frac12(X_1+iY_1)\otimes\frac12(X_6-iY_6)\\
&=\frac14(X_1X_6+Y_1Y_6+iX_1Y_6-iY_1X_6)+\frac14(X_1X_6+Y_1Y_6-iX_1Y_6+iY_1X_6)\\
&=\frac12(X_1X_6+Y_1Y_6)
\end{align*}
$$

となります。以上より、証明したい式である

$$
(a_1^\dagger a_6+a_6^\dagger a_1)\mapsto\frac {(-1)^{N_f-1}} 2(X_1X_6+Y_1Y_6)
$$

が得られました。また、 $(i,j)=(7,12)$ の場合についても、

$$
\begin{align*}a_{12}^\dagger a_7&\mapsto
(\otimes_{k=1}^{11}Z_k)\otimes\sigma_{12}^\dagger (\otimes_{k=1}^6Z_k)\otimes\sigma_7\\
&=(\otimes_{k=1}^6Z_k^2)\otimes Z_7\sigma_7\otimes(\otimes_{k=8}^{11}Z_k)\otimes\sigma_{12}^\dagger\\
&=(\otimes_{k=1}^6I_k)\otimes Z_7\sigma_7\otimes(\otimes_{k=8}^{11}Z_k)\otimes\sigma_{12}^\dagger
\end{align*}
$$

のように、1番目から6番目の演算子が恒等演算子となり無視できることから、 $(i,j)=(1,6)$ の場合と同様の議論が成立し、

$$
(a_7^\dagger a_{12}+a_{12}^\dagger a_7)\mapsto\frac {(-1)^{N_f-1}} 2(X_7X_{12}+Y_{7}Y_{12})
$$

が導けます。

参考文献

[1] N. F. Mott, The basis of the electron theory of metals, with special reference to the transition metals, Proc. Phys. Soc. A 62, 416 (1949).

[2] S. Yamada, T. Imamura, and M. Machida, 16.447 TFlops and 159-Billion-dimensional exact-diagonalization for trapped Fermion-Hubbard model on the earth simulator, SC ’05: Proceedings of the 2005 ACM/IEEE Conference on Supercomputing, IEEE, (2005).

[3] M. Cerezo, A. Arrasmith, R. Babbush, S. C. Benjamin, S. Endo, K. Fujii, J. R. McClean, K. Mitarai, X. Yuan, L. Cincio, P. J. Coles, Variational quantum algorithms, Nat. Rev. Phys. 3, 625 (2021).

[4] Y. Alexeev, M. Amsler, M. A. Barroca, et al., Quantumcentric supercomputing for materials science: A perspective on challenges and future directions, Future Generation Computer Systems 160, 666 (2024).

[5] E. Farhi, J. Goldstone, and S. Gutmann, A quantum approximate optimization algorithm, arXiv:1411.4028 (2014).

[6] P. Jordan and E. Wigner, Über das Paulische Äquivalenzverbot, Z. Phys. A. 47, 631 (1928).

[7] N. Hatano and Masuo Suzuki, “Finding exponential product formulas of higher orders” In Quantum annealing and other optimization methods, (Springer, Berlin, Heidelberg, 2005) pp. 37-68.

[8] G. Ortiz, J. E. Gubernatis, E. Knill, and R. Laflamme, Quantum algorithms for fermionic simulations, Phys. Rev. A 64, 022319 (2001).

[9] Qiskit, https://quantum.cloud.ibm.com/docs/en/guides

あとがき

FTQCアルゴリズムとされている量子位相推定をもとに作られたIQPEが、現状の量子コンピューターでもうまく機能することが興味深いと感じました。本論文で扱ったのは6原子グラフェンというシンプルな系でしたが、ハバード模型を用いた他の様々な系に適用することができれば、研究の幅が大きく広がると思いました。

齊藤 孝介

東京科学大学 理学院 物理学系 学士課程