Table of contents
第4章 量子ダイナミクスシミュレーション
シュレディンガー方程式とは、量子ダイナミクスとは
ダイナミクスのシミュレーション
物理系のシミュレーションのゴールは、系(着目している対象)の初期状態(最初の状態)と、その系のダイナミクス＝時間発展(系が時間の経過に伴ってどう変化するか)のルールが与えられたときに、一定時間後に系がどういう状態にあるかを知ること。 系の時間発展のルールは、一般に微分方程式で表される。古典力学系ならニュートンの運動方程式、電磁気学系ならマクスウェル方程式、流体力学系ならナビエ・ストークス方程式、といった具合である。したがって、物理系のシミュレーションの本質は、微分方程式を数値的に解くことに帰着される。そして、当然ながら量子系のシミュレーション、すなわち量子力学に従う系のシミュレーションもその中に含まれる。
量子系のシミュレーション
量子系を古典コンピュータでシミュレーションする場合、系のサイズ（粒子やスピンの数）が大きくなると、計算にかかる時間が指数関数的に増大するため、すぐに太刀打ちできなくなってしまう。例えば、n個の量子ビットのある系をシミュレーションするためにはおよそ\(2^n\)個の方程式を解かなければならない(参照：Nielsen-Chuang 4.7 Simulation of quantum systems)。指数関数的な速度向上が期待される量子コンピュータを用いて量子系のシミュレーションを行うことで、上述の問題を解決しようというアイデアが生まれた。
用語整理
シュレディンガー方程式
\[i\hbar \frac{\partial \left| \psi (t)\right>}{\partial t} = H \left| \psi (t) \right>\]ハミルトニアン
| シュレディンガー方程式に出てくる演算子であり、系の時間発展を司る。量子コンピュータでシミュレーションする範囲においては、ただの大きな行列だと思っておいても良い。第1章で学んだ量子操作と異なる点は、ハミルトニアンはエルミート\(H^\dagger = H\)だということ。そのため、一般にハミルトニアンを作用させた状態$$H \left | \psi (t) \right>$$を量子コンピュータ上で作成することはできない(ユニタリかつエルミートであるパウリ演算子X, Y, Zを除く)。 |
ハミルトニアンは、系のエネルギーに対応するオブザーバブルでもある。ハミルトニアンには、粒子同士に働く力や、外部からの力（電場・磁場など）といった系の全ての情報が詰まっている。このハミルトニアンをシュレディンガー方程式に入れ、それを解くことで、系のエネルギーの値が得られる。
オブザーバブル(物理量)
量子力学では、物理量(エネルギー、磁化=系の示す磁力の大きさなど)は状態に作用するエルミート演算子として記述される。演算子Aがエルミートであるとは\(A^\dagger = A\)を満たすこと。 Aを対角化した時の固有値を\(\{ a_i\}_i\)、固有ベクトルを\(\{ \left| a_i \right>\}_i\)とする。状態\(\left| \psi \right>\)について、物理量Aを観測(射影測定)したとき、観測値\(a_i\)が
\[\left| \psi \right> = \sum_i c_i \left| a_i \right>\]| の係数$$ | c_i | ^2$$に比例した確率で得られる(1章の射影測定はZというオブザーバブルに関する射影測定)。 |
エネルギー固有状態(基底状態・励起状態)
ハミルトニアン\(H\)は系のエネルギーを表すオブザーバブルなので、その固有値がわかれば系が取り得るエネルギーの値(固有エネルギー)がわかる。系のエネルギーがわかると「実際にどの状態が実現しやすいか」が分かる。与えられたハミルトニアンのもとでシュレディンガー方程式を解くと、一般に複数の解(状態とエネルギーの組)が得られる。それらの中で、最も低いエネルギーを持つ、すなわち最も安定な状態は基底状態と呼ばれ、自然の中で最も実現しやすい状態である(系の温度が上がると、より高エネルギーな状態も出現しやすくなる)。
ダイナミクス(時間発展)
ハミルトニアンが時間に依存しない系のシュレディンガー方程式を解くと
\[\left| \psi (t) \right> = e^{-iHt/\hbar} \left| \psi (0) \right>\]| となる。量子ダイナミクスシミュレーションとは、つまるところこの右辺をいかに計算するかという問題である。ハミルトニアンHは量子ビットの数nに対して\(2^n\)次元と指数的に大きくなるので、古典コンピュータでは非常に大変。よくある問題設定は、とある初期状態$$\left | \psi (0) \right>\(からスタートし、時間\)t$$だけ経過した後にあんらかの物理量Aの期待値 |
を計算する、というもの。
von-Neumann方程式
混合状態と呼ばれる状態についてのダイナミクスを記述する式。シュレディンガー方程式に似ている。
トロッター分解を用いた量子シミュレーション
この手法はダイナミクスのシミュレーションとしてはある意味最も愚直な手法であり、NISQデバイスでも一定の規模の問題までは動作する。ここでは1次元横磁場イジングモデルという物理系の量子ダイナミクスをシミュレートしてみよう。
(リー・)トロッター分解とは
正方行列A, Bの和の指数関数を、それぞれの指数関数の積に近似する以下の公式のこと。
\[e^{\delta (A+B)} = e^{\delta A} \cdot e^{\delta B} + \mathcal{O} (\delta^2)\](A, Bは行列なので\(e^{A+B} \neq e^A \cdot e^B\)に注意)
トロッター分解を用いた量子シミュレーションの仕組み
量子状態に時間発展はシュレディンガー方程式
\[i\frac{\partial }{\partial t} \left| \psi (t) \right> = H \left| \psi (t) \right>\]に従う(\(\hbar = 1\))とした。特に、ハミルトニアンHが時間に依存しない場合は、これを解くことで\(\left| \psi (t) \right> = e^{-iHt} \left| \psi (0) \right>\)となる。すなわち、量子系のダイナミクスのシミュレーションでは\(e^{-iHt}\)という演算子が計算できれば良い。しかしn量子ビット系だとHは\(2^n\)次元となり、非常に大きな計算となる。ゆえに古典コンピュータでこの指数関数を計算するのは困難。そして量子コンピュータでもあらゆる\(2^n\)次元の行列Hの指数関数を解くことは難しいと考えられている。
しかし、ハミルトニアンが特定の構造を持つ場合、量子コンピュータは\(e^{-iHt}\)を効率よく計算できる。それに用いるのがトロッター分解。
具体例
n量子ビットからなるイジングモデル
\[H = \sum_{i=1}^{n-1} Z_i Z_{i+1}\]を考える。このHは\(2^n\)次元の行列だが、トロッター分解を用いると
\[e^{-iHt} = e^{-i(\sum_{t=1}^{n-1} Z_i Z_{i+1})t} = \left( e^{-i(\sum_{t=1}^{n-1} Z_i Z_{i+1})t \frac{1}{M}} \right)^M \simeq \left( e^{-i Z_1 Z_2 \frac{t}{M}} \cdot e^{-i Z_2 Z_3 \frac{t}{M}} \cdots \right)^M \tag{1}\]ここでMは\(e^{-iHt}\)の作用を分割する総数であり、近似の精度\(\mathcal{O}((t/M)^2)\)が十分に小さくなるように選ぶ。
最右辺に着目すると小さい行列\(e^{-i Z_i Z_{i+1}\frac{t}{M}}\)のnM個の積で書かれている。よって指数的に大きな行列の作用を小さな行列(2量子ビットゲート)の多項式個の作用の積で書き表すことができた。
このようにして、ハミルトニアンが少数の項の和でかける場合には、トロッター分解を用いて量子コンピュータで高速に計算を行うことが可能。
物理・量子化学の分野で興味あるハミルトニアンは、大体は効率的なトロッター分解ができる形をしている。
量子ダイナミクスの実装(1): イジングモデル
物質が磁場に反応する性質を磁性と呼ぶ。ミクロな世界では電子や原子核が非常に小さな磁性を持ち、スピンと呼ばれる。大雑把にはスピンとは小さな磁石のことであると思って良い。私たちが日常で目にする磁石は、物質中の電子(原子核)のスピンが示す小さな磁性が綺麗に揃った結果、巨大な磁性が現れたものである。
イジングモデルはそのような物質中のスピンの振る舞いを記述するモデルで、磁石の本質を取り出したモデルとして考案された。定義は
である。ここでJはモデルのパラメータ、nはスピンを持つ粒子の数、\(Z_i\)はi番目の粒子のスピンのz軸方向の成分。前節での記述より、物理量はオブザーバブルというエルミート演算子で表現されるので、スピンが示す磁化(磁力の大きさ)も演算子で表現される。空間は3次元だから、磁化の独立な成分としてx, y, z軸の3方向の磁化が考えられ、それぞれに対応する演算子はX, Y, Zである(パウリ演算子X, Y, Zの名前の由来はここからきている)。大まかにいううとi番目の粒子=i番目の量子ビットである。
\[Z \left| 0\right\rangle = \left| 0\right\rangle, \ Z \left| 1\right\rangle = -\left| 1\right\rangle, \\]| より、$$\left | 0\right\rangle\(はスピンが+z方向を向いている状態、\)\left | 1\right\rangle$$はスピンが-z方向を向いている状態を表している。 | |
| このモデルを構成する\(Z_i Z_{i+1}\)は隣り合う量子ビット間の相互作用を表している。J>0のとき、\((Z_i, Z_{i+1}) = (1, 1) {\rm or} (-1, -1)\)で\(Z_i Z_{i+1} = 1\)、そして\((Z_i, Z_{i+1}) = (1, -1) {\rm or} (-1, 1)\)で\(Z_i Z_{i+1} = -1\)となる。よってエネルギーが低いのはスピンが互い違いになる$$\left | 010101 \cdots \right\rangle\(のような場合である。逆にJ<0ではスピンが揃った場合\)\left | 000 \cdots \right\rangle, \left | 111 \cdots \right\rangle$$の方が、エネルギーが低く安定になり、系は巨大な磁化を持った磁石として振る舞う。 |
実はこのモデルそのものは量子コンピュータを使わなくても簡単に解けてしまう。しかしここではこれをトロッター分解してその基本を実装してみよう。
\(e^{-iHt}\)のトロッター分解は(1)式で導いている。必要になるのは\(e^{-i\delta Z_i Z_{i+1}}\)というゲートをどのように実装するかである。これは次のようにすれば良い(証明は後述)。
今回はz方向の全磁化
\[\left\langle \psi \left| \sum_{i=1}^n Z_i \right| \psi \right\rangle\]| という物理量に着目してその時間変化を追いかけてみよう。初期状態はスピンが+z方向に揃った$$\left | 000 \cdots \right\rangle\(とする。また簡単のため、系が周期的になっているとして\)Z_{i+1} = Z_1$$とする。結果は以下のようになる。 |
z方向の全磁化は時間に依存しない。イジングモデルではz方向の相互作用しかないので、z方向の全磁化も保存する。
次はx軸方向に磁場を加えた横磁場イジングモデルの時間発展を考えてみよう。
量子ダイナミクスの実装(2): 横磁場イジングモデル
イジングモデルにx軸方向の一様な磁場をかけたモデル
\[H = \sum_{i=1}^n Z_i Z_{i+1} + h\sum_{i=1}^n X_i\]hは横磁場の強さを表す係数、そして\(X_i\)はi番目の粒子のx方向の磁化を表すパウリ演算子(オブザーバブル)。このモデルでトロッター分解を行うと、先ほどのイジングモデルに加えて\(e^{-ihX_1\frac{t}{M}} \cdots e^{-ihX_n\frac{t}{M}}\)という項が加わる。実際に解いた結果が下図。
h=0のイジングモデルの場合は全磁化が一定だったのに対して、横磁場イジングモデルは全磁化の値が負になることもある。これは横磁場(x方向の磁場)の影響でz方向に揃っていたスピンが乱れて乱雑になり、結果として全磁化の値が変化していると考えれば良い。
量子ダイナミクスの実装(3): 厳密解との比較
トロッター分解には誤差がある。上で計算したダイナミクスがどれほどの制度なのか、\(e^{-iHt}\)を直接計算した厳密なダイナミクスと比べてみる。
途中式の証明
| 具体的に4つの状態、すなわち$$\left | 00 \right\rangle, \left | 01 \right\rangle, \left | 10 \right\rangle, \left | 11 \right\rangle$$に作用させたときに両辺が同じものになれば良い。 |
| ここで$$Z \left | 0 \right\rangle = \left | 0 \right\rangle, Z \left | 1 \right\rangle = - \left | 1 \right\rangle$$より |
一方右辺からは
\[e^{ - i\delta \overbrace{Z_i Z_{i+1}}^{1}} \left| 00 \right\rangle = e^{-i\delta } \left| 00 \right\rangle\]となる。これを他のものにも計算してやれば良い。すると
\[{\rm CNOT}_{i, i+1} \cdot ({\bf 1} \otimes e^{-i\delta Z_{i+1}}) \cdot {\rm CNOT}_{i, i+1} \left| 01 \right\rangle = {\rm CNOT}_{i, i+1} \cdot ({\bf 1} \otimes e^{-i\delta Z_{i+1}}) \left| 01 \right\rangle = {\rm CNOT}_{i, i+1} e^{i\delta} \left| 01 \right\rangle = e^{i\delta} \left| 01 \right\rangle\] \[{\rm CNOT}_{i, i+1} \cdot ({\bf 1} \otimes e^{-i\delta Z_{i+1}}) \cdot {\rm CNOT}_{i, i+1} \left| 10 \right\rangle = {\rm CNOT}_{i, i+1} \cdot ({\bf 1} \otimes e^{-i\delta Z_{i+1}}) \left| 11 \right\rangle = {\rm CNOT}_{i, i+1} e^{i\delta} \left| 11 \right\rangle = e^{i\delta} \left| 10 \right\rangle\] \[{\rm CNOT}_{i, i+1} \cdot ({\bf 1} \otimes e^{-i\delta Z_{i+1}}) \cdot {\rm CNOT}_{i, i+1} \left| 11 \right\rangle = {\rm CNOT}_{i, i+1} \cdot ({\bf 1} \otimes e^{-i\delta Z_{i+1}}) \left| 10 \right\rangle = {\rm CNOT}_{i, i+1} e^{-i\delta} \left| 10 \right\rangle = e^{-i\delta} \left| 11 \right\rangle\]が、全て右辺の結果と一致することがわかる。よって等式が成り立つ。
参考文献(というよりこれを見て勉強中)