スピンのラーモア歳差運動を量子計算でシミュレーション

$$
\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}}
$$

はじめに

シュレーディンガー方程式の時間発展は、量子回路で表現することができます(前回記事)。その応用として、スピンのラーモア歳差運動をとりあげ、自作の量子計算シミュレータqlazyを使って、時間発展の様子を確認してみました。

ラーモア歳差運動については、以下を参考にさせていただきました。

• EMANの物理学:量子力学:スピンのラーモア歳差運動

ラーモア歳差運動

回転するコマが斜めになりながら、倒れずに鉛直軸周りにくるりくるり回る運動のことを「歳差運動」といいます。静磁場中に置かれた円運動する電荷(=磁気モーメント)にも同じ現象が起きます。原理は両方とも同じで、重力や静磁場によって引き倒されようとするのですが、逆に倒されまいとする力が同時に働くことによって、重力や静磁場方向を軸に回転する運動が発生します。磁気モーメントの場合の回転運動を「ラーモア歳差運動」というようです。

スピンは磁気モーメントに似ていますが、ちょっと違います。が、スピンの場合も「ラーモア歳差運動」が発生することがわかっています。それを量子計算で再現・確認するのが、今回のテーマです。

シミュレーション(ブロッホ球上での量子状態の軌跡)

いま、X軸プラスの方向に静磁場をかけ、そこに１個のスピンが置かれているとします。で、スピンの初期状態は$\ket{0}$、つまりZ軸プラス方向としておきます。このとき、シュレーディンガー方程式のハミルトニアンは単純にパウリ行列$X$そのものです。この前提で、波動関数$\ket{\psi(t)}$の時間発展を、量子計算シミュレータqlazyで計算していきます。

以下にコードを示します。

【2021.9.5追記】qlazy最新版でのソースコードはここに置いてあります。

import math
import cmath
from qlazypy import QState,Observable

def bloch(alpha,beta):

    # eliminate phase factor
    beta = beta / alpha
    alpha = 1.0

    # normalize
    norm = math.sqrt(1.0 + abs(beta)**2)

    alpha /= norm
    beta /= norm

    # (alpha,beta) => (theta,phi)
    theta = (2.0 * math.acos(alpha))
    if abs(alpha)<0.00001 or abs(beta)<0.00001: # north pole or south pole
        phi = 0.0
    else:
        phi = -cmath.phase(math.sin(theta/2.0) / beta)

    # unit => pi*radian
    theta /= math.pi
    phi /= math.pi
    
    return theta, phi
    
def main():

    hm = Observable("x_0")

    for i in range(21):
        t = i*0.05
        qs = QState(1)

        # time evolution
        qs.evolve(observable=hm, time=t, iter=100)

        # quantum stat in bloch spere
        theta, phi = bloch(qs.amp[0],qs.amp[1])

        print("time = {0:.2f}, theta = {1:.2f}*PI, phi = {2:.2f}*PI"
              .format(t,theta,phi))

        del qs
    
    del hm
    
if __name__ == '__main__':
    main()

何をやっているか、少し説明します。

from qlazypy import QState,Observable

量子状態を格納しておく「QStateクラス」に加えて、ハミルトニアンなどのオブザーバブルを格納するための「Observableクラス」をimportしておきます(v0.0.9で追加しました)。

main()関数で、

hm = Observable("x_0")

としていますが、これで今回のハミルトニアンXを定義しています。引数の"x_0"は「0番目の粒子に対応したパウリ行列X」という意味です。もしハミルトニアンが$-2.0 X_{0} \otimes Y_{1} + Z_{2}$だったら、引数文字列は"-2.0x_0y_1+z_2"と指定します。

量子状態の初期値は、

qs = QState(1)

で$\ket{0}$になります。時間発展は、evolveメソッドを使って、

qs.evolve(observable=hm, time=t, iter=100)

のように指定します。ここでtimeで時間、iterで近似の繰り返し回数を指定します。リー・トロッターの積公式の指数に相当する整数値で、近似の精度を上げるためには、時間timeよりも十分大きなiterを設定するのが良いです。今回、時間を0.0から1.0まで変化させてみました。

theta, phi = bloch(qs.amp[0],qs.amp[1])

で、時間発展させたqsが、ブロッホ球上でどの位置にあるのかを計算しています。つまり、

\begin{align}
\ket{\psi} &= \alpha \ket{0} + \beta \ket{1} \\
&= e^{i \gamma} (\cos \frac{\theta}{2} \ket{0} + e^{i \phi} \sin \frac{\theta}{2} \ket{1})
\end{align}

で定義された$(\theta,\phi)$の組を求めています。引数として指定したqs.amp[0],qs.amp[1]は、上の式で$\alpha,\beta$と書かれているものです。bloch()関数はコードの上の方に定義されている通りです(説明省略)。

後は、この$(\theta,\phi)$を順に表示しているだけです。

さて、結果はどうなるかというと、

time = 0.00, theta = 0.00*PI, phi = 0.00*PI
time = 0.05, theta = 0.10*PI, phi = 0.50*PI
time = 0.10, theta = 0.20*PI, phi = 0.50*PI
time = 0.15, theta = 0.30*PI, phi = 0.50*PI
time = 0.20, theta = 0.40*PI, phi = 0.50*PI
time = 0.25, theta = 0.50*PI, phi = 0.50*PI
time = 0.30, theta = 0.60*PI, phi = 0.50*PI
time = 0.35, theta = 0.70*PI, phi = 0.50*PI
time = 0.40, theta = 0.80*PI, phi = 0.50*PI
time = 0.45, theta = 0.90*PI, phi = 0.50*PI
time = 0.50, theta = 1.00*PI, phi = 0.00*PI
time = 0.55, theta = 0.90*PI, phi = -0.50*PI
time = 0.60, theta = 0.80*PI, phi = -0.50*PI
time = 0.65, theta = 0.70*PI, phi = -0.50*PI
time = 0.70, theta = 0.60*PI, phi = -0.50*PI
time = 0.75, theta = 0.50*PI, phi = -0.50*PI
time = 0.80, theta = 0.40*PI, phi = -0.50*PI
time = 0.85, theta = 0.30*PI, phi = -0.50*PI
time = 0.90, theta = 0.20*PI, phi = -0.50*PI
time = 0.95, theta = 0.10*PI, phi = -0.50*PI
time = 1.00, theta = 0.00*PI, phi = 0.00*PI

となります。ちょっと、見づらくてすみませんが、主にthetaの値に注目していただけると、きちんと回転(歳差運動)していることがわかります。phiは最初0.5で、thetaは0.0(北極)から次第に大きくなり、theta=1.0(南極)を境にphiが-0.5になり、逆にthetaの値が減少していきます。つまり、YZ平面上を１回転しています。説明するまでもないですが、北極と南極のphiは単純に計算すると発散するので0.0にしています。それと、角度の単位は$\pi$ラジアンです。なので、phi=0.5は$0.5\pi$ラジアンを表します。

シミュレーション(期待値)

ついでに、期待値も計算してみました。量子状態がYZ平面上を動くので、YとZの期待値も同様の動きをするはずです。ということの確認です。

コードは以下です。

【2021.9.5追記】qlazy最新版でのソースコードはここに置いてあります。

import math
from qlazypy import QState,Observable

def main():

    hm = Observable("x_0")
    ob_y = Observable("y_0")
    ob_z = Observable("z_0")

    for i in range(6):
        t = i*0.05
        qs = QState(1)

        # time evolution
        qs.evolve(observable=hm, time=t, iter=100)

        # expectation values of Observable Y,Z
        exp_y = qs.expect(observable=ob_y)
        exp_z = qs.expect(observable=ob_z)

        # argument
        if abs(exp_z)<0.00001: # exp_theta = 0.5*PI, if exp_z = 0.0 
            exp_theta = 0.5
        else:
            exp_theta = math.atan(exp_y/exp_z) / math.pi

        print("time = {0:.2f}, <y> = {1:.2f}, <z> = {2:.2f}, <theta> = {3:.2f}*PI"
              .format(t, exp_y, exp_z, exp_theta))

        del qs
    
    del hm
    
if __name__ == '__main__':
    main()

何をやっているか、少し説明します。

ob_y = Observable("y_0")
ob_z = Observable("z_0")

ハミルトニアンに加え、YとZをオブザーバブルとして定義しておきます。

時間発展は先程と同様で、そのqsに対して、expectメソッドで期待値を計算します。引数にはYとZのオブザーバブルを指定します。

exp_y = qs.expect(observable=ob_y)
exp_z = qs.expect(observable=ob_z)

あとは、この期待値から偏角を計算して、表示しています。

さて、結果はどうなるかというと、

time = 0.00, <y> = 0.00, <z> = 1.00, <theta> = 0.00*PI
time = 0.05, <y> = 0.31, <z> = 0.95, <theta> = 0.10*PI
time = 0.10, <y> = 0.59, <z> = 0.81, <theta> = 0.20*PI
time = 0.15, <y> = 0.81, <z> = 0.59, <theta> = 0.30*PI
time = 0.20, <y> = 0.95, <z> = 0.31, <theta> = 0.40*PI
time = 0.25, <y> = 1.00, <z> = -0.00, <theta> = 0.50*PI

となります(定義域の処理など面倒だったので最初の６サンプルだけにしました、汗)。ここで、<y>,<z>はY,Zに対応した期待値の時間変化です。<theta>は、それから割り出した偏角(y/zのarctan)です。先程のブロッホ球上の動きと一致していることがわかります。

ということで、確認終了です。めでたしめでたし。

おわりに

スピンが１個しかないとても単純な系ですが、とりあえず量子計算で量子状態の時間発展を正しく再現することができました。もう少し複雑な系で何か試してみたい気がしてきますが、今後の課題としておきます。

それにしても、古典力学的な法則が全く関与しない、純粋にシュレーディンガー方程式だけで、こんな運動が導けるのは、少し不思議な気がします。

以上