QuTiPでフロケ系のシミュレーション

# 概要
量子フロケ系とは、時間依存ハミルトニアン$H(t)$が時間$T$の周期を有している量子系のことです。
```math
H(t + T) = H(t)
```

量子系のシミュレーションパッケージである[QuTiP](https://qutip.org/)には、フロケ系のシミュレーションに必要なツールが揃っているので、それを使ってみます。

# 理論的背景
時間依存するハミルトニアン系のシュレディンガー方程式は次のように書けます。
```math
i \hbar \frac{d}{dt} |\psi(t)\rangle = H(t) |\psi(t)\rangle
```
この解は時間発展演算子$U(t_2, t_1)$を用いて
```math
|\psi(t)\rangle = U(t, 0)|\psi(0)\rangle
```
と書けます。ただし、
```math 
U(t_2, t_1):= \mathcal{T}\left[ \exp \left(-\frac{i}{\hbar} \int_{t_1}^{t_2} ds H(s) \right)\right] 
```
であり、$\mathcal{T}$は[時間順序積](https://ja.wikipedia.org/wiki/%E6%99%82%E9%96%93%E9%A0%86%E5%BA%8F%E7%A9%8D)を表します。

## 時間発展演算子の性質
$t_1$から$t_2$までの時間発展と、$t_2$から$t_3$までの時間発展を合わせると$t_1$から$t_3$までの時間発展となるので、次の式が成り立ちます。
```math 
U(t_3, t_2)U(t_2, t_1)= U(t_3, t_1)
```
また$H$の周期性から、$U(t_2, t_1)$の$T$の時間並進に対する不変性がわかります。
```math
\begin{align} 
U(t_2+T, t_1+T) &= \mathcal{T}\left[ \exp \left(-\frac{i}{\hbar} \int_{t_1+T}^{t_2+T} ds H(s) \right)\right] \\
& = \mathcal{T}\left[ \exp \left(-\frac{i}{\hbar} \int_{t_1}^{t_2} ds H(s + T) \right)\right] \\
& = \mathcal{T}\left[ \exp \left(-\frac{i}{\hbar} \int_{t_1}^{t_2} ds H(s ) \right)\right] \\
& = U(t_2, t_1)
\end{align}
```

## フロケの定理
$U(T, 0)$はユニタリ演算子なので、あるエルミート演算子$H_F$を用いて、
```math 
U(T, 0) = \exp(-i H_F T/\hbar)
```
と表すことができます。また、$U(t,0)\exp(i H_F t/\hbar)$もユニタリ演算子なので、時刻$t$に依存するあるエルミート演算子$K_F(t)$を用いて、
```math 
U(t,0)\exp(i H_F t/\hbar) = \exp(-i K_F(t)/\hbar)
```
と表すことができます。以上より、
```math 
U(t,0) = \exp(-i K_F(t)/\hbar)\exp(-i H_F t/\hbar)
```
が成り立ちます。
ここで、
```math 
\begin{align}
U(t+T,0) & = U(t+T, T)U(T, 0)\\
& = U(t, 0)U(T, 0)\\
& = \exp(-i K_F(t)/\hbar)\exp(-i H_F t/\hbar)\exp(-i H_F T/\hbar)\\
& = \exp(-i K_F(t)/\hbar)\exp(-i H_F\times (t+T)/\hbar)
\end{align}
```
であり、左辺は$\exp(-i K_F(t+T)/\hbar)\exp(-i H_F\times (t+T)/\hbar)$であることに注意すると、$K_F(t)$の周期性
```math 
K_F(t+T) = K_F(t)
```
がわかります。
$H_F$の固有値を{$\epsilon_n$}、固有状態を{$|n\rangle$}とし、初期状態$|\psi(0)\rangle$を{$|n\rangle$}で展開したときの係数を{$a_n$}と置いてみます。
```math 
|\psi(0)\rangle = \sum_n a_n |n\rangle
```
時刻$t$の状態$|\psi(t)\rangle$は、
```math 
\begin{align}
|\psi(t)\rangle & = U(t, 0)|\psi(0)\rangle \\
&= \sum_n a_n U(t,0)|n\rangle \\
&= \sum_n a_n  \exp(-i K_F(t)/\hbar)\exp(-i H_F t/\hbar)|n\rangle \\
&= \sum_n a_n  \exp(-i K_F(t)/\hbar)\exp(-i \epsilon_n t/\hbar)|n\rangle \\
&= \sum_n a_n \exp(-i \epsilon_n t/\hbar) |u_n(t)\rangle \\
\end{align}
```
となります。ただし、
```math 
|u_n(t)\rangle　:= \exp(-i K_F(t)/\hbar)|n\rangle 
```
であり、$K_F(t)$の周期性から、$|u_n(t)\rangle$も周期性を持ちます。
```math 
|u_n(t+T)\rangle　=|u_n(t)\rangle
```
このように$|\psi(t)\rangle$が周期$T$を持つ$|u_n(t)\rangle$を用いて$\exp(-i \epsilon_n t/\hbar) |u_n(t)\rangle $の重ね合わせで書けることを、**フロケの定理**と言います。

以上より、数値計算において$|\psi(t)\rangle$は次の手順で求めることができます。

1. $U(T,0)$を計算し、対角化することで{$\epsilon_n$}、{$|n\rangle$}を求める。
2. $t\in [0, T]$に対して$|u_n(t)\rangle = \exp(i \epsilon_n t/\hbar) U(t, 0)|n\rangle$ を求める。
3. $|\psi(0)\rangle$の{$|n\rangle$}による展開係数{$a_n$}を求め、$|\psi(t)\rangle = \sum_n a_n \exp(- i \epsilon_n t/\hbar) |u_n(t')\rangle$を計算する。ただし、$t'\in [0, T]$は$t$を$T$で割った余り。

# QuTiPでの計算
次のハミルトニアンを取り扱います。
```math 
H(t)=− \frac{\Delta}{2} \sigma_x − \frac{\epsilon_0}{2}   \sigma_z + \frac{A}{2} \sin(\omega t)\sigma_z.
```

まず{$|n\rangle = |u_n(0)\rangle$}と{$\epsilon_n$}を求めてみましょう。2次元なので2つの固有状態と固有値が生じます。

```python
import numpy as np
from qutip import *
import matplotlib.pyplot as plt

delta = 0.2 * 2*np.pi
eps0 = 1.0 * 2*np.pi
A = 2.5 * 2*np.pi
omega = 1.0 * 2*np.pi
T = 2*np.pi / omega

H0 = - delta/2.0 * sigmax() - eps0/2.0 * sigmaz()
H1 = A/2.0 * sigmaz()
H = [H0, [H1, 'sin(w * t)']] # H0 + H1*sin(w*t)の意味。時間依存項を表すQuTipの記法。
args = {'w': omega} # これをfloquet_modesの引数に渡すと、wにomegaが代入される。
f_modes_0, f_energies = floquet_modes(H, T, args)
print(f_energies)
print(f_modes_0[0].data)
print(f_modes_0[1].data)
```

実行結果
```
[-2.83131212  2.83131212]
  (0, 0)	(0.72964231187078+0j)
  (1, 0)	(-0.39993745864107016+0.5546820043084018j)
  (0, 0)	(0.39993745864107005+0.5546820043084018j)
  (1, 0)	(0.72964231187078+0j)
```

{$\epsilon_n$}は定義上$2\pi \hbar /T$の整数倍の任意性がありますが、QuTiPでは$-\pi\hbar /T \le \epsilon_n \le \pi\hbar /T$に制限されているようです。`f_modes_0`は[Qobj](https://qutip.org/docs/latest/apidoc/classes.html#qobj)クラスのリストです。`data`属性で量子状態のスパース行列を取得できます。

次に時間発展を計算してみましょう。
初期状態を
```math 
|\psi(0)\rangle = 
\begin{pmatrix}
1\\
0
\end{pmatrix}
```
として、展開係数{$a_n$}を計算してみます。

```python
psi0 = basis(2,0) 
f_coeff = floquet_state_decomposition(f_modes_0, f_energies, psi0)
print(f_coeff)
```
実行結果
```
[(0.72964231187078+0j), (0.39993745864107005-0.5546820043084018j)]
```

引数の`f_energies`は別に必要ないし、フロケとも関係ない関数ですが、QuTiPで与えられた基底で展開係数求めるのはこの関数しかないようです…

状態$|\psi(t)\rangle $は次のように計算できます。

```python
psi_t = floquet_wavefunction_t(f_modes_0, f_energies, f_coeff, 1, H, T, args) #時刻1での状態
print(psi_t.data)
```

実行結果
```
  (0, 0)	(-0.952247951725807+0.019771617832092026j)
  (1, 0)	(-0.2471425782149468-0.17819502685421923j)
```

最後に数演算子$N := \sigma_+ \sigma_-$の期待値を$t \in [0, 10T]$の範囲でプロットしてみます。

```python
tlist  = np.linspace(0.0, 10 * T, 101)
p_ex = np.zeros(len(tlist))
for n, t in enumerate(tlist):
    psi_t = floquet_wavefunction_t(f_modes_0, f_energies, f_coeff, t, H, T, args)
    p_ex[n] = expect(num(2), psi_t)

plt.plot(tlist, real(p_ex), 'o')
plt.xlabel('Time')
plt.ylabel('Occupation probability')
plt.legend(("Floquet $P_1$"))
plt.show()
```

![ダウンロード.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/271560/b08f9856-fce2-db39-f9be-fdb337388b76.png)

# 参考リンク
QuTiP公式チュートリアル
https://qutip.org/docs/4.1/guide/dynamics/dynamics-floquet.html

Konrad Viebahn, Introduction to Floquet theory (2020).
https://boulderschool.yale.edu/sites/default/files/files/floquetlecture_Konrad%20Viebahn.pdf