[Qutip] pythonで量子状態の虚時間発展をシミュレートする

はじめに

虚時間発展法のシミュレートの記事がなかったのでメモとして残します．
今回は，Qutipによる素直な実装と，Qiskitに実装されている変分法による近似であるVarQITEの実装を行います.
-> まずはQutipでの実装のみで，VarQITEはまた追記します．

Qutipによる実装

qutip == 5.0.4

ハミルトニアンを定義します．ここでは，横磁場イジング模型のハミルトニアンを使用します．問題のサイズをNとします．

from qutip import sigmaz, sigmax, qeye, tensor

# 問題のサイズ = 量子ビット数
N = 2

# 演算子のリスト
sz = [tensor([sigmaz() if i == j else qeye(2) for j in range(N)]) for i in range(N)]
sx = [tensor([sigmax() if i ==  j else qeye(2) for j in range(N)]) for i in range(N)]

問題ハミルトニアンのインスタンスを作成します．横磁場イジング模型は，以下の形をしていますので，　素直に実装します．今回は，相互作用係数$J_{ij}$および局所磁場$h_i$はランダム$(0〜1)$にします．
$$
H = \sum_{i<j} J_{ij} \sigma_i^z \sigma_j^z + \sum_{i} h_i \sigma_i^x
$$

import random

H = 0
for i in range(N):
    for j in range(N):
        H += random.random() * sz[i] * sz[j] 

for i in range(N):
    H += random.random() * sx[i]

H

上のコードを実行すると，ランダムな係数で定義されたQobjクラスのインスタンスが作成されるはずです．

（実行例）

それでは，虚時間発展を行う関数を定義します．今回は素直に，虚時間シュレディンガー方程式
$$
\dfrac{\partial \ket{\psi(\tau)}} {\partial \tau} = - \hat{H} \ket{\psi(\tau)},
$$
において，ハミルトニアンが時間依存しない場合の解
$$
\ket{\psi(\tau)} = A(\tau) e^{-H \tau} \ket{\psi (0)} = A(\tau) \sum_{i} c_i e^{-E_i \tau } \ket{\phi_i}
$$
を用います．ここで$A(\tau)$は規格化定数です．

def SimulateImaginaryTimeEvolution(H, psi, T):
    """
    量子状態psiの虚時間発展をシミュレートする関数
    """
    num_steps = 100
    dt = T / num_steps
    psi_history = [psi] # psi:初期状態
    psi_evol =  psi.copy()
    
    for _ in range(num_steps):
        U =  (-1 * H * dt).expm() # e^{-H*dt}
        psi_evol = (U * psi_evol).unit() # 虚時間発展演算子はnon-unitaryなので規格化
        psi_history.append(psi_evol)
    return psi_evol, psi_history # 終状態とhistoryをreturn

では実行してみましょう．適当な初期状態と発展時間Tを用意してあげます．

from qutip import basis

ket0, ket1 = basis(2,0), basis(2,1) # |0〉 と |1〉
init_state = tensor([(ket0+ket1).unit() for _ in range(N)]) # |++..+〉

T = 10
final_state, history = SimulateImaginaryTimeEvolution(H, init_state, T)

プロットしてみましょう．ついでに，H.eigenenergies()を用いて求めた基底エネルギーの厳密解もプロットします．

from matplotlib import pyplot as plt
from qutip import expect

tlist = [t for t in range(int(100)+1)]
plt.plot(tlist, [expect(H, psi) for psi in history], label='QITE') # 各状態におけるHの期待値を計算

plt.axhline(y=H.eigenenergies()[0], color='r', linewidth=0.7, label='Ground state energy') # 基底状態のエネルギーを赤の横線でプロット

plt.title('Imaginary time evolution')
plt.xlabel('imaginary time')
plt.ylabel('energy expectation value')
plt.legend()
plt.show()

できました．基底状態に収束していそうです．
念の為，基底状態$\ket{\phi_0}$とのフィデリティもプロットします．なおフィデリティは以下の式で計算されます．
$$
fidelity = | \braket{\phi_0 |\psi (\tau)}|^2
$$

ground_state =  H.eigenstates()[1][0]
plt.plot(tlist, [abs(psi.dag()*ground_state)**2 for psi in history], label='QITE')

plt.title('fidelity with ground state')
plt.ylabel('fidelity')
plt.xlabel('imaginary time')
plt.legend()
plt.show()

できました．