3.4. 【例題】量子スピン系（ハイゼンベルグ模型）¶

3.4.1. ハイゼンベルグ模型 ¶

スピン間の交換相互作用を考慮した以下のハイゼンベルグ模型は磁性の研究でよく用いられます：

(1)¶ $\mathcal{H} = \frac{1}{2} \sum_{ij} J_{ij} \bm{S}_i \cdot \bm{S}_j$

ここで、 $i, j$ はスピンを区別するラベルで1からNの整数とします。 $\bm{S}_i=(S_i^x, S_i^y, S_i^z)$ は量子力学のスピン演算子です。異なるラベルのスピン演算子は交換します。したがって、完全な交換関係は

$[S_i^{a}, S_j^{b}] = i\epsilon_{abc} S_i^{c} \delta_{ij}$

で与えられます。ここで、 $a, b=x, y, z$ で $\epsilon_{abc}$ は完全反対称テンソルです（ $\epsilon_{xyz}=+1$ ）。

ハイゼンベルグ模型では、2体相互作用で全てのスピンが繋がっているので、全系の多体状態ベクトルを考える必要があります。具体的には、 $\sigma_i=\uparrow, \downarrow$ として、

$| \sigma_1 \sigma_2 \cdots \sigma_N \rangle$

で状態ベクトルが表されます。ベクトル空間の次元が $2^N$ なので、ハミルトニアンは $2^N \times 2^N$ 行列です。 10スピン程度（ $2^{10} \approx 10^3$ ）なら対角化して全ての固有値・固有ベクトルを計算することができます。また、最大固有値のみを求めるアルゴリズム（例えばべき乗法やランチョス法がある [6]）を使えば20スピン（ $2^{20} \approx 10^6$ ）以上を扱えます。しかし、熱力学極限にはまだまだ遠いです。多体問題を正確に解くことの難しさが分かると思います。

以下で、N=2と3の場合を解いて、量子スピン系の感覚を掴みましょう。ここでは、NumPyの線形代数モジュールを使って数値的に対角化をしますが、SymPy というライブラリを使うと解析計算も可能です。興味のある人は調べてみてください。

3.4.2. 2スピン模型 ¶

まずは練習のために、2スピン系から考えます。 (1) 式で $J_{12}=J_{21}=J$ とすると

(2)¶ $\mathcal{H} = J \bm{S}_1 \cdot \bm{S}_2$

この系の状態は $|\uparrow \uparrow \rangle$ , $|\uparrow \downarrow \rangle$ , $|\downarrow \uparrow \rangle$ , $|\downarrow \downarrow \rangle$ の4つなので、ハミルトニアンは4×4行列です。 J>0のとき、スピンが反平行の状態（反強磁性状態）のエネルギーが低くなります。そのような状態は $|\uparrow \downarrow \rangle$ と $|\uparrow \downarrow \rangle$ の2つあり、基底状態はそれらの重ね合わせです。数値計算で確認しましょう。

実装 ¶

ソースコード two_spins.py

まずは1スピン演算子を表す2×2行列を作ります。

def make_spin_ops():
    ops = {}

    # S^z
    ops['Sz'] = np.diag((0.5, -0.5))

    # S^+
    sp = np.zeros((2,2))
    sp[0, 1] = 1
    ops['S+'] = sp

    # S^-
    sm = sp.transpose()
    ops['S-'] = sm

    # Sx, Sy
    ops['Sx'] = (sp + sm) / 2.0
    ops['Sy'] = (sp - sm) / 2.0j

    # S^0
    ops['I'] = np.identity(2)

    return ops

戻り値は辞書型です。 pythonでは複数の戻り値を返すことができますが、この例のようにひとつの辞書にまとめて返す方法もあります。こうすると、戻り値の順番を意識しなくてよいという点や、後から戻り値を追加した場合に、この関数を既に呼び出しているコードを変更する必要がないという点がメリットです。

注釈

辞書型を使って複数の戻り値を返す方法以外に、collections.namedtuple という名前付きのタプル（c言語の構造体のようなもの）を使う方法もあります。詳しくは Effective Python 2nd Edition [3] の "Item 19: Never unpack more than three variables when functions return multiple values" を参照してください。

10行目：numpy.diag は対角行列を作る関数です。 13行目：numpy.zeros は0で初期化された配列を生成する関数です（numpy.empty で生成すると初期化されません）。 0で初期化した後、14行目で必要な要素に値を代入します。 18行目：numpy.transpose は行列の転置を取る関数です。 26行目：numpy.identity は単位行列を生成する関数です。 NumPyの様々な配列生成法は、リファレンスマニュアルの Array creation routines にまとまっています。

次に、2スピン演算子を作ります。以下のように、クロネッカー積 numpy.kron を使えば簡単に実装できます：

def make_matrix(s1, s2):
    return np.kron(s1, s2)

2つの2×2行列を入力すると4×4行列が得られます。このように、NumPyには便利な演算が多く提供されています。リファレンスマニュアルの Linear algebra (numpy.linalg) を参照してください。

注釈

クロネッカー積については Wikipedia - クロネッカー積を参照してください。

ここまで準備ができたらあと少しです。 (2) 式のハミルトニアンを表す4×4行列を作って対角化をします。

def main():
    sp_ops = make_spin_ops()
    sz = sp_ops['Sz']
    sp = sp_ops['S+']
    sm = sp_ops['S-']
    s0 = sp_ops['I']

    hamil = make_matrix(sz, sz) + (make_matrix(sp, sm) + make_matrix(sm, sp)) / 2.0
    print("H =\n", hamil)

    eigval, eigvec = np.linalg.eigh(hamil)

    print("\nEigenvalues =\n", eigval)
    print("\nEigenvectors =")
    for i in range(4):
        print(eigvec[:, i])

45行目で対角化をして固有値と固有ベクトルを受け取っています。 numpy.linalg.eigh 関数はエルミート行列の対角化ルーチンです。固有値が必ず実数になります。 50行目：固有ベクトルは縦ベクトルとして保存されていることに注意してください。

注釈

numpy.linalg.eigh 関数は、入力された行列がエルミートかどうかはチェックしないので注意してください。エルミート行列を仮定して、左下半分の要素のみを使います（UPLO='L' オプション）。安全なコードにするなら、eighを呼び出す前に、行列がエルミートになっているかを自分で判定してください。

結果 ¶

H =
 [[ 0.25  0.    0.    0.  ]
 [ 0.   -0.25  0.5   0.  ]
 [ 0.    0.5  -0.25  0.  ]
 [ 0.    0.    0.    0.25]]

Eigenvalues =
 [-0.75  0.25  0.25  0.25]

Eigenvectors =
[ 0.          0.70710678 -0.70710678  0.        ]
[1. 0. 0. 0.]
[0.         0.70710678 0.70710678 0.        ]
[0. 0. 0. 1.]

エネルギーが-3/4のスピン一重項とエネルギーが1/4のスピン三重項に分かれていることが確認できます。

3.4.3. 3スピン模型 ¶

さて、ここまでは練習です。次に3スピン系を考えます。

(3)¶ $\mathcal{H} = J (\bm{S}_1 \cdot \bm{S}_2 + \bm{S}_2 \cdot \bm{S}_3 + \bm{S}_3 \cdot \bm{S}_1)$

状態ベクトルは8つ、ハミルトニアンは8×8行列です。 N=2の場合と異なり、全ての相互作用エネルギーを得するような反強磁性状態は作れません。例えば、 $|\uparrow \downarrow \uparrow \rangle$ 状態は、 $\bm{S}_1 \cdot \bm{S}_2$ と $\bm{S}_2 \cdot \bm{S}_3$ の相互作用エネルギーは得しますが、 $\bm{S}_3 \cdot \bm{S}_1$ で損をしてしまいます。これを 幾何学的フラストレーション と呼びます。古典的な反強磁性状態が作れないので、量子力学的な効果が顕著になります。 N=3の場合に、基底状態で非自明な縮退が残ることを見てみましょう。

実装 ¶

ソースコード three_spins.py

1スピン演算子の生成には、先ほど作った make_spin_ops 関数を再利用しましょう。

from two_spins import make_spin_ops

このようにインポートすれば使えます。ただし、2つのソースファイルは同じ場所においてください。

$N=3$ の状態空間における相互作用の行列表現は、クロネッカー積を2回使って作ることができます。

def make_matrix(s1, s2, s3):
    return np.kron(s1, np.kron(s2, s3))

$2 \times 2$ 行列を3つ入力すると、 $8 \times 8$ 行列が生成されます。例えば、 $S_1^x S_2^x$ の場合、 $S_1^x S_2^x = S_1^x S_2^x I_3$ と考えて（ $I$ は単位行列）、make_matrix(sx, sx, s0) とします。

メイン部分は、2スピンの場合からハミルトニアンの部分だけ書き換えれば出来上がりです。

def main():
    sp_ops = make_spin_ops()
    sz = sp_ops['Sz']
    sp = sp_ops['S+']
    sm = sp_ops['S-']
    s0 = sp_ops['I']

    hamil = make_matrix(sz, sz, s0) + (make_matrix(sp, sm, s0) + make_matrix(sm, sp, s0)) / 2.0 \
          + make_matrix(s0, sz, sz) + (make_matrix(s0, sp, sm) + make_matrix(s0, sm, sp)) / 2.0 \
          + make_matrix(sz, s0, sz) + (make_matrix(sm, s0, sp) + make_matrix(sp, s0, sm)) / 2.0
    print("H =\n", hamil)

    eigval, eigvec = np.linalg.eigh(hamil)

    print("\nEigenvalues =\n", eigval)
    print("\nEigenvectors =")
    for i in range(8):
        print(eigvec[:, i])

結果 ¶

H =
 [[ 0.75  0.    0.    0.    0.    0.    0.    0.  ]
 [ 0.   -0.25  0.5   0.    0.5   0.    0.    0.  ]
 [ 0.    0.5  -0.25  0.    0.5   0.    0.    0.  ]
 [ 0.    0.    0.   -0.25  0.    0.5   0.5   0.  ]
 [ 0.    0.5   0.5   0.   -0.25  0.    0.    0.  ]
 [ 0.    0.    0.    0.5   0.   -0.25  0.5   0.  ]
 [ 0.    0.    0.    0.5   0.    0.5  -0.25  0.  ]
 [ 0.    0.    0.    0.    0.    0.    0.    0.75]]

Eigenvalues =
 [-0.75 -0.75 -0.75 -0.75  0.75  0.75  0.75  0.75]

Eigenvectors =
[-0.         -0.81622341  0.42640143  0.          0.38982197 -0.
 -0.         -0.        ]
[ 0.          0.          0.         -0.81622341  0.          0.42640143
  0.38982197  0.        ]
[ 0.          0.02111916  0.69631062  0.         -0.71742978  0.
  0.          0.        ]
[ 0.          0.          0.         -0.02111916  0.         -0.69631062
  0.71742978  0.        ]
[-0.         -0.57735027 -0.57735027  0.         -0.57735027 -0.
 -0.         -0.        ]
[-0.         -0.         -0.          0.57735027  0.          0.57735027
  0.57735027 -0.        ]
[1. 0. 0. 0. 0. 0. 0. 0.]
[0. 0. 0. 0. 0. 0. 0. 1.]

4重縮退が2つ（エネルギー+3/4と-3/4）という結果になりました。

スピン3つの合成状態は、全スピンS=3/2, S=1/2, S=1/2の3種類で、縮退度はそれぞれ4重、2重、2重です。 S=3/2は強磁性状態で、高エネルギー状態（エネルギー+3/4）が対応しています。基底状態（エネルギー-3/4）の4重縮退は、S=1/2が2つ縮退していることを意味しています。この結果は、全スピン $\bm{S}^2$ と方位量子数 $S_z$ に加えて、もう一つハミルトニアンと可換な演算子（量子数）が存在すること示唆します。

実は、以下の演算子がハミルトニアンと可換です：

(4)¶ $C=\bm{S}_1 \cdot (\bm{S}_2 \times \bm{S}_3)$

これはスカラースピンカイラリティと呼ばれます。この演算子の固有値の符号から、3つのスピンの相対的な向き（右手・左手）が区別されます。

注釈

(4) 式の右辺はスカラー三重積と呼ばれる形です。その絶対値が3つのベクトルで張られる平行六面体（平行四辺形の3次元版）の体積を与え、符号は3つのベクトルが右手系の関係にあればプラス、左手系の関係にあればマイナスです。また、どれか2つのベクトルが平行または反平行だとゼロになります。

【発展】交換関係 ¶

ソースコード three_spins_chiral.py

演算子 $C$ がハミルトニアンと交換可能であることを数値計算で確かめます。

まずは、全スピン演算子で交換関係の計算をして、正しく計算できることを確認します。答えの分からない問題を解く前に、答えの分かっている問題でテストをすることは数値計算ではとても重要です。全スピン演算子 $\bm{S}^2, S^x, S^y, S^z$ を以下のように作ります：

    # total spin: Sx, Sy, Sz, S^2
    tot_sx = make_matrix(sx, s0, s0) + make_matrix(s0, sx, s0) + make_matrix(s0, s0, sx)
    tot_sy = make_matrix(sy, s0, s0) + make_matrix(s0, sy, s0) + make_matrix(s0, s0, sy)
    tot_sz = make_matrix(sz, s0, s0) + make_matrix(s0, sz, s0) + make_matrix(s0, s0, sz)
    tot_s = tot_sx @ tot_sx + tot_sy @ tot_sy + tot_sz @ tot_sz
    print("S^x =\n", tot_sx)
    print("S^2 =\n", tot_s)

ここでも make_matrix 関数を使って、1スピン演算子を表す2×2行列から3スピン系の8×8行列を作ります。 32行目：行列積には @ 演算子を使います（numpy.matmul 関数と等価です）。

次に、交換関係を判定する関数を用意します。

def does_commute(m1, m2):
    """Return True if m1 and m2 commute"""
    return np.allclose(m1 @ m2, m2 @ m1)

この関数は2つの行列 m1, m2 が交換可能な場合に True を返します。配列の比較には numpy.allclose を使います。

注釈

numpy.allclose 関数は、2つの配列の全ての要素が ある精度の範囲で 等しい場合に True を返します。精度は相対誤差 rtol と絶対誤差 atol があり、後者のデフォルトは atol=1e-08 です。コンピュータの実数は 丸め誤差 を含むので、実数同士の比較は a==b とせずに、必ず abs(a-b)<atol として誤差を考慮に入れます（abs は絶対値を返す関数）。

NumPy変数に対して様々な判定を行う関数は NumPy Reference の Logic functions にまとまっています。一度目を通しておくとよいでしょう。

この関数を使って交換関係を計算してみると

    # commutation relations
    print("[H, Sz]  :", does_commute(hamil, tot_sz))
    print("[H, S^2] :", does_commute(hamil, tot_s))
    print("[S^2, Sz]:", does_commute(tot_s, tot_sz))
    print("[Sz, Sx] :", does_commute(tot_sz, tot_sx))

結果：

[H, Sz]  : True
[H, S^2] : True
[S^2, Sz]: True
[Sz, Sx] : False

期待通りの結果が得られました。計算の確認のため、交換しない例もひとつ入れました。

さて、それでは目的のスカラースピンカイラリティ演算子 (4) を調べます。 $C$ は3つのスピンを含む演算子ですが、これまでと同様に make_matrix 関数を使って作れます。

    # scalar spin chirality, C
    scalar_chiral_op = make_matrix(sx, sy, sz) - make_matrix(sx, sz, sy) \
                     + make_matrix(sy, sz, sx) - make_matrix(sy, sx, sz) \
                     + make_matrix(sz, sx, sy) - make_matrix(sz, sy, sx)
    print("\nC =\n", scalar_chiral_op)

    # commutation relations
    print("[H, C]   :", does_commute(hamil, scalar_chiral_op))
    print("[C, Sz]  :", does_commute(scalar_chiral_op, tot_sz))
    print("[C, S^2] :", does_commute(scalar_chiral_op, tot_s))

結果：

C =
 [[0.+0.j   0.+0.j   0.+0.j   0.+0.j   0.+0.j   0.+0.j   0.+0.j   0.+0.j  ]
 [0.+0.j   0.+0.j   0.+0.25j 0.+0.j   0.-0.25j 0.+0.j   0.+0.j   0.+0.j  ]
 [0.+0.j   0.-0.25j 0.+0.j   0.+0.j   0.+0.25j 0.+0.j   0.+0.j   0.+0.j  ]
 [0.+0.j   0.+0.j   0.+0.j   0.+0.j   0.+0.j   0.-0.25j 0.+0.25j 0.+0.j  ]
 [0.+0.j   0.+0.25j 0.-0.25j 0.+0.j   0.+0.j   0.+0.j   0.+0.j   0.+0.j  ]
 [0.+0.j   0.+0.j   0.+0.j   0.+0.25j 0.+0.j   0.+0.j   0.-0.25j 0.+0.j  ]
 [0.+0.j   0.+0.j   0.+0.j   0.-0.25j 0.+0.j   0.+0.25j 0.+0.j   0.+0.j  ]
 [0.+0.j   0.+0.j   0.+0.j   0.+0.j   0.+0.j   0.+0.j   0.+0.j   0.+0.j  ]]
[H, C]   : True
[C, Sz]  : True
[C, S^2] : True

確かに、 $C$ がハミルトニアンと交換することが確認できました。また、 $C$ は $\bm{S}^2, S^z$ とも交換します。このことから、ハミルトニアン (3) が3つの演算子 $\bm{S}^2, S^z, C$ と同時対角化可能であることがわかりました。

注釈

ここでは同時対角化可能であることだけ示しました。ハミルトニアンの固有状態を実際に $\bm{S}^2, S^z, C$ の量子数で分類するにはここで示した方法だけでは不十分です。どのようにしたらできるか考えてみてください。

3.4. 【例題】量子スピン系（ハイゼンベルグ模型）¶

3.4.1. ハイゼンベルグ模型 ¶

3.4.2. 2スピン模型 ¶

実装 ¶

結果 ¶

3.4.3. 3スピン模型 ¶

実装 ¶

結果 ¶

【発展】交換関係 ¶

目次

前のトピックへ

次のトピックへ