[QuTiP]ブロッホ球のアニメーションの計算方法

概要

前回の記事で演算子の定義の方法を解説したので，今回はブロッホ球のアニメーションを描写してみます。

今回のコードはすべてこちらにありますので，ご参考に。

ブロッホ球とは

wikipediaの項目を読んでもらった方が正しい記述になると思いますが，私なりの言葉で記載しておきます。
ブロッホ球は，スピンなどの2つの量子状態しかとらない対象を球面上に表す「表記方法」です。
2つの量子状態の時間発展ダイナミクスを球状の動きとして記述することで理解（想像？）しやすいことが利点です。

このブロッホ球で描かれる対象は色々あります。
私が初めてブロッホ球に触れたのは医療現場で画像診断に使われているNMRを理解するためでした。
NMRはスピンエコーを利用して体のスキャン画像を作成しています。
このスピンのダイナミクスを理解するためにブロッホ球を利用していました。

また，量子コンピュータの量子ビットは古典ビットと同じくon-offの二値しかとらないので，その時間発展ダイナミクスはブロッホ球面上の動作として描写できます。
microsoftのQ#のチュートリアルのpauli operatorの項目にもブロッホ球が出てきます。
ここではゲート操作として出てくるので，知っておくべき内容かなと思います。

基本の描写方法

QuTiPのドキュメントによれば，ブロッホ球のクラスがQuTiPには実装されています。

from qutip import Bloch
b = Bloch()
b.show()

と，このようにBloch球が描写できます。保存方法は

b.save(filename)

です。Z方向が取ることができる二つの量子状態になっています。
また，適切な操作をすることでx方向やy方向にも初期状態を作ることができます。

適当な操作を行って任意の初期状態を作ったとします。初期状態を描写するためにこのブロッホ球面上にベクトルを書いてみます。|0>方向のベクトルを書きます。

コード

b = Bloch()
pnt = [0,0,1]
b.add_vectors(pnt)
b.show()
b.save("/Bloch_zvector.png")

pntでx, y, zの位置の情報を与えて，add_vectorsで与えて位置へ中心からベクトル
を描くことができます。
当然ですが，pntをの要素を変えれば描かれるベクトルも変わります。例えば，pnt = [-1/2, -1/2, 1/2]であれば，

となります。

脇道にそれました。この初期状態|0>のベクトルから次第に変化していく様子を点で描写したいので，次は点の描写のやり方をみてみます。これも簡単で，add_pointsメソッドを使えば良いです。

コード

b = Bloch()
pnt = [0,0,1]
b.add_vectors(pnt)
pnt = [0,1/2,1/2]
b.add_points(pnt)
b.show()
b.save(data_dir+"/Bloch_vec_poi.png")

青色の点が描写されているのがわかります。
これらを利用して計算を行った結果を可視化してみます。

計算結果

いきなりですが，計算結果のアニメーションは次のようになります。

このアニメーションは，先のドキュメントのexampleのコードを実行したものになります。

初期状態$|0\rangle$が励起状態で，$|1\rangle$への緩和の様子です。
ここではベクトルの方向に沿って，系の運動が行われます。¹

このあたりの話を真面目に書き出すと趣旨から大きく脱線するので，

• 全体的な外観が知りたければScullyのQuantum Optics
• 緩和の導出や古典統計物理との対応をなどの丁寧な記述であれば，CarmichaelのStatistical Methods in Quantum Optics 1: Master Equations and Fokker-Planck Equations
• 緩和の話を中心に知りたければ，GardinerとZollerのQuntum Noise
  私はCarmichaelの本が好きでした。

本筋に戻って，ここで解いているハミルトニアンは

\mathcal{H} = \omega (\cos\theta \sigma_z + \sin\theta \sigma_x)

であり，緩和を表すマスター方程式は

\dot{\rho} = -i[\mathcal{H}, \rho] + \Gamma(n+1)(2\sigma_+\rho\sigma_- [\sigma_+\sigma_-, \rho])\\
+\Gamma(n+1)(2\sigma_-\rho\sigma_+ -[\sigma_-\sigma_+, \rho])\\
+\gamma(\sigma_z\rho\sigma_z -2\rho)

です。この方程式を解いて，上のアニメーションのようにブロッホ球面上のダイナミクスを計算しています。

これらのダイナミクスを解くためのコードは，

def qubit_integrate(w, theta, gamma1, gamma2, psi0, tlist):
    # 演算子の定義
    sx = sigmax(); sy = sigmay(); sz = sigmaz(); sm = sigmam()
    # ハミルトニアン
    H = w * (cos(theta) * sz + sin(theta) * sx)
    # マスター方程式
    c_op_list = []
    n_th = 0.5 # temperature
    rate = gamma1 * (n_th + 1)
    if rate > 0.0: c_op_list.append(sqrt(rate) * sm)
    rate = gamma1 * n_th
    if rate > 0.0: c_op_list.append(sqrt(rate) * sm.dag())
    rate = gamma2
    if rate > 0.0: c_op_list.append(sqrt(rate) * sz)


    # evolve and calculate expectation values
    output = mesolve(H, psi0, tlist, c_op_list, [sx, sy, sz])  
    return output.expect[0], output.expect[1], output.expect[2]

に当たります。外部から，共鳴周波数$\omega$と$\theta$，緩和レート$\Gamma$, $\gamma$，初期状態$\phi_0$を受け取っています。

演算子の定義

    # 演算子の定義
    sx = sigmax(); sy = sigmay(); sz = sigmaz(); sm = sigmam()

は前回の記事のままですね。
ハミルトニアン

    # ハミルトニアン
    H = w * (cos(theta) * sz + sin(theta) * sx)

は，先程記したものをそのまま打ち込んだものになります。

マスター方程式は

    n_th = 0.5 # temperature
    rate = gamma1 * (n_th + 1)
    if rate > 0.0: c_op_list.append(sqrt(rate) * sm)

で，マスター方程式の1つ目の緩和を表すことができています。
ここまで定義すれば，解きたい時間と解きたい期待値がどれかmsolveメソッドの入れてあげることで計算ができます。

    output = mesolve(H, psi0, tlist, c_op_list, [sx, sy, sz])  

これの第一引数がハミルトニアン，第二引数が初期状態，第三引数が解きたい時間，第四引数が緩和，第五引数が解きたい期待値がどれか，になります。
第五引数の順番に答えが帰ってくるので，

    return output.expect[0], output.expect[1], output.expect[2]

で計算結果を返してくれます。私はずっと手でマスター方程式を解いていたので，なんて簡単なんだ！と感動しました。

def animate(i):
    sphere.clear()
    sphere.add_vectors([np.sin(theta),0,np.cos(theta)])
    sphere.add_points([sx[:i+1],sy[:i+1],sz[:i+1]])
    sphere.make_sphere()
    return ax

結果を一点一点アニメーションにしているのが次のコードになります。
ダイナミクスの方向をadd_vectorsメソッドで定義し，時間ごとに少しずつadd_pointsメソッドで点を増加させています。

これをmatplotlibのanimation.FuncAnimationで動画化しています。
最後にmp4で保存して，今回の動画を作りました。

ani.save('bloch_sphere.mp4', fps=20)

ただ，このままでは，Qiitaに貼れなかったので，ffmpegを使ってmp4をgifに変換しました。²

#　まとめ

• ブロッホ球について簡単に説明しました。
• QuTiPのブロッホ球の描写方法を記載しました。
• ドキュメントどおりにアニメーションを作成しました。

1. 多分。天下り的に計算をしていたので検証はしていません。問題があれば指摘してほしいです。
   また，$|0\rangle$→$|1\rangle$への緩和を表すために，量子マスター方程式を使ってダイナミクスを解いています。 ↩

2. 例えば　https://qiita.com/wMETAw/items/fdb754022aec1da88e6e ↩