ハイゼンベルグの運動方程式をPythonを使って解けないか？

QuTiPを使えという話ですが，NumPyやSympyをつかって第二量子化したハミルトニアンのハイゼンベルグの運動方程式を立式できないだろうかというのが気になったのでやってみる。
sympyの演算子などのドキュメントが英語しかなかったので，それも含めてちょうどいいかなと思った。（ただpythonを使って第二量子化した演算子を取り扱おうとしていて，日本語しか読めないなんて人はまずいないだろう）

anacondaにデフォルトで入ってる，sympyモジュールにはブラケットを定義，ブラケットの内積や外積，交換関係を計算してくれるモジュールがあったので，テストしてみた。

計算モデル

ただの２準位系を想定して計算をする。ハミルトニアンを
$$\mathcal{H}=\hbar\omega|1\rangle\langle 1|$$
と定義して，例えば$|0\rangle\langle1|$の運動方程式が導出できるかやってみる。ここで，$|0\rangle$, $|1\rangle$はそれぞれ２準位系の基底状態と励起状態で，$\hbar,~\omega$は換算プランク定数と２準位系の共鳴周波数。

sympyのモジュールを利用したハミルトニアンの記述

$|0\rangle$, $|1\rangle$を定義してみる。

from sympy.physics.quantum import Ket

up=Ket('1')
low=Ket('0')

print(up,low)

を実行すると，期待通りに

|1> |0>

が返ってきたので，外積モジュールを利用してハミルトニアンを作成する。

from sympy.physics.quantum import Dagger,Ket, OuterProduct
from sympy import Symbol
up=Ket('1')
low=Ket('0')
freq=Symbol('omega')
plank_const=Symbol('hbar')
H=plank_const*freq*OuterProduct(up,Dagger(up))

print(H)

を実行して，

hbar*omega*|1><1|

ここで利用したsympyモジュールDaggerは複素共役を計算してくれるもので，OuterProductでブラとケットの外積の計算をしてくれる。

ここまででハミルトニアンをpython上で書くことができた。

ハイゼンベルグの運動方程式を解いてみる

本題

Commutator($\mathcal{O}$,$\mathcal{H}$)で計算を行ってくれることを期待した。

    from sympy.physics.quantum import Commutator, Dagger
    from sympy.physics.quantum import Ket, OuterProduct
    from sympy import Symbol
    up=Ket('1')
    low=Ket('0')
    freq=Symbol('omega')
    plank_const=Symbol('hbar')


    upper=OuterProduct(up,Dagger(low))
    H=plank_const*freq*OuterProduct(up,Dagger(up))
    result=-1/plank_const*Commutator(upper,H)

    print(result)

ここでupperは上昇（生成？）演算子の$|1\rangle\langle 0|$で，最終行は単なるハイゼンベルグの運動方程式（のつもり）だ。結果は

-omega*[|1><0|,|1><1|]

となって，計算が出来ていない。sympyのドキュメントを読んでいると

.expand(commutator=True).doit()

をつければよいと書いていたので，実行してみる。

-omega*(|1><0|*|1><1| - |1><1|*|1><0|)

無事に途中まで計算ができた。状態は直行しているので，一つ目の項が0になるはずだが，sympyのブラケットの定義だけでは自動でクロネッカーのデルタに書き換えてくれないみたい。

何やらいい方法があれば教えてほしい。QuTipを使えばいいんだろうなあ