$$
\def\bra#1{\mathinner{\left\langle{#1}\right|}}
\def\ket#1{\mathinner{\left|{#1}\right\rangle}}
\def\braket#1#2{\mathinner{\left\langle{#1}\middle|#2\right\rangle}}
$$

はじめに

QuTIPという量子状態を計算するライブラリがある。
これは、量子化学計算などとは異なり、離散化された量子状態のダイナミクスを計算するのに用いるライブラリらしい。少々興味をもって遊んで見たので、簡単な使い方を書いておきたい。
なお、本記事に関するコードは、ここに置いてある。

（注）ここで計算したモデルはあくまで仮想的なものであり、あくまでQuTIPをとりあえずいじってみるための簡単なコードだと思っていただきたい。

簡単な説明

この記事では、もっとも簡単なモデルとして、（Z軸方向の）スピンが$+\frac{1}{2}$と$-\frac{1}{2}$の２状態のみをとる系を考える。この状態を、それぞれ$\ket{0}$と$\ket{1}$とおくこととしよう。
量子情報などを少しでも触ったことがある方はQbitが一つだけある状態だと思ってもらえば良い。

この系の波動関数は、次の形式の波動関数で表される。係数$a$、$b$は複素数で、その自乗が観測時にそれぞれの状態を見出す確率となる。

\psi(t) = a(t) \ket{0} + b(t) \ket{1}

ここで、

\braket{0}{0} = \braket{1}{1} = 1\\
\braket{0}{1} = \braket{1}{0} = 0\\
|a|^2 + |b|^2 = 1

である。波動関数の時間発展は、シュレーディンガー方程式
$$
i\hbar \frac{d}{dt}\psi(t) = \hat{H} \psi(t)
$$
により記述される。これを上記の波動関数で展開すると、$a,b$の時間発展は、次のように表される。

i\hbar\frac{d}{dt} \begin{pmatrix}a(t)\\b(t)\end{pmatrix} = \hat{H}\begin{pmatrix}a(t)\\b(t)\end{pmatrix}

ここに、Z軸方向に磁場をかけた場合と、X軸方向に磁場をかけた場合の系の様子を、次節でQuTIPを使って計算してみよう。

問題設定

ハミルトニアン

系のエネルギーを表す演算子である、ハミルトニアン行列($\hat{H}$)をまずは作る。
まずは、Z軸方向（正確には、スピンの量子化軸と同じ方向）に磁場がかかった状態を考えよう。
このとき、ハミルトニアン行列は下記のように表される。

\hat{H} = \mu B \sigma_z

なお、$\sigma_z$は、z軸方向に関するパウリ行列である。大きさBの磁場がかかった時に、スピン$\pm \frac{1}{2}$では感じるエネルギーの符号が逆になることを考えるとわかりやすい。

\sigma_z = \begin{pmatrix}1 && 0 \\ 0 && -1 \end{pmatrix}

以上で構成した$\hat{H}$行列は非対角成分が0なので、時間発展は行列表示をせずとも下記のように書き下した方が早い。それぞれの振幅は、独立に動いていることがわかる。

\frac{d}{dt}\begin{pmatrix}a \\ b \end{pmatrix} = \frac{i}{\hbar}\begin{pmatrix}a \\ -b \end{pmatrix}

初期状態

はじめはシンプルに、+1/2の状態のみが存在する状態$ \ket{\psi} = \ket{0} $を初期状態と置いてみよう。
実際のところは、これだとすでに状態が確定している状態と言え、その時間発展を追うことに何の意味がある、ということになるのだが、計算してみるのは自由である。

QuTIPでコード化

インストールはpipで簡単にインストールができる。

pip install qutip

今回はjupyter-notebookを利用した。グラフの可視化にはmatplotlibを用いる。必要ライブラリのインポートなどは下記のようになる。

%matplotlib inline
import matplotlib.pyplot as plt
import numpy as np
from qutip import *

次に、前節で記述したハミルトニアンと初期状態を記述する。この部分は、数式と簡単に対応がつけられると思う。
初期状態を定義するのに用いているbasis()関数は、系のとる状態を生成する関数である。第一引数は取りうる状態数を表しており、今回は$\ket{0}$のみを初期状態としているので第二引数はそれを表すために0としている。

# Hamiltonian
B = 1.0
H = B * sigmaz()

#initial state
psi0 = basis(2, 0)

最後に、計算を行う時間間隔を設定し、計算を行う。計算に置いては、mesolveというソルバを用いる。
この第５引数でエルミート演算子を配列にして渡すことで、各時間ステップでの演算子の値を取り出すことができるようになる。また、第４引数に空の配列を渡しているが、これは系の緩和を記述するものである。今回は、系の閉じた状態のみを扱っているため、気にしない。

tlist = np.linspace(0, 10, 100)

#dynamics calculation
result = mesolve(H, psi0, tlist, [], [sigmax(), sigmay(), sigmaz()])

この時の$\sigma_z, \sigma_x, \sigma_y $の演算子の期待値を表示してみよう。

fig, axes = plt.subplots(1,1)
axes.plot(tlist, result.expect[0], label = r'$\left<\sigma_x\right>$')
axes.plot(tlist, result.expect[1], label = r'$\left<\sigma_y\right>$')
axes.plot(tlist, result.expect[2], label = r'$\left<\sigma_z\right>$')
axes.legend(loc=2)

グラフが何も変化しておらず、あまり面白くない。
系の状態が$\ket{0}$の状態に確定した状態を初期状態として、しかもハミルトニアンの非対角要素が0でこれらの状態間の行き来がないため、$\sigma_z$は当然ながら1になる。

初期状態を変えてみる

次に、ハミルトニアンを同じまま、$\pm\frac{1}{2}$の状態が重ね合った状態に変えて計算してみよう。
つまり、系の初期状態の波動関数をとしてみる。

\ket{\psi} = \frac{1}{\sqrt{2}}(\ket{0} + \ket{1})

psi0 = (basis(2, 0) + basis(2,1)).unit()   #  super position : (|0> + |1>)

これで時間発展をさせた結果は下記のグラフのようになる。この時は、$\sigma_z$は決まった値を取り続けているが、$\sigma_x, \sigma_y$はそれぞれ時間とともに変化をしていることがわかる。

磁場の向きを変えてみる

今度は、考えているスピンの量子化軸とは垂直な方向に磁場をかけた状態を考えてみる。
この時のハミルトニアンを構成するコードは、X方向のパウリ行列$\sigma_x$を用いて次のようになる。

H = B * sigmax()

\sigma_x = \begin{pmatrix} 0 && 1 \\ 1 && 0 \end{pmatrix}

$\ket{0}$のみの状態を初期状態として計算を行ったときの結果のグラフを示す。

今度は、ハミルトニアンの非対角要素が0でないことから、$\sigma_z$の期待値が1と-1の間を行ったり来たりしていることがわかる。

最後に

以上で見てきたようにQuTIPでハミルトニアンを構成して、初期状態を変えたりしながら計算を行うことで、簡単に系の時間発展を追うことができるようになる。今回のケースだと、磁場の方向を変えてみたり、あるいは複数方向に同じ、または異なる大きさの磁場を変えてみることでさらに色々と試してみることで、数式だけといてもピンとこない部分などが見えてくるだろう。

なお、最後にとても言い訳がましいのだが、実際には磁場をかけている方向が量子化軸となるので、２番目の例は実験で再現するのは不可能である。（NMRなどは、多数の原子・分子の統計集団のもつ磁化をみている。）

今後は、緩和のある系などについても試してみたい。

参考文献

QuTIP (http://qutip.org/)
R.P.ファインマン、R.B.レイトン、M.サンズ(en)「ファインマン物理学 <5> 量子力学」　共立出版

ここで使った関数などに関しては、やはりQuTIPのマニュアルなどを見ていただくのが良いと思う。
二状態系に関しては量子光学の分野などで広く研究されているため、様々な本がある。私ははじめにファインマン物理学　第５巻のアンモニア・メーザーの例題で読み、非常に分かり易かった覚えがあるのでこれを参考文献とした。