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

## はじめに

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

## Qutipによる実装
 ```
 qutip == 5.0.4
 ```
ハミルトニアンを定義します．ここでは，横磁場イジング模型のハミルトニアンを使用します．問題のサイズをNとします．
```python
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
$$

```python
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クラス`のインスタンスが作成されるはずです．

（実行例）
![スクリーンショット 2024-10-15 23.59.00.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3469107/0cc65c3f-d851-0aac-9340-2147f168a069.png)

それでは，虚時間発展を行う関数を定義します．今回は素直に，虚時間シュレディンガー方程式
$$
        \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)$は規格化定数です．

```python
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を用意してあげます．
```python
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()`を用いて求めた基底エネルギーの厳密解もプロットします．
```python
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()
```

![output.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3469107/2b010553-cf40-c57b-33e1-fbc6323c5c70.png)

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

```python
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()
```
![output2.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3469107/863992dd-b1e2-17dc-3c4e-2da723e0572c.png)

できました．