光連続量計算と調和振動子から位相平面プロットまで

光連続量計算について学び始めました。永井君の社内資料やQiitaを中心に学んでいます。凄いわかりやすいので、ほぼ98%永井君の資料をトレースしてみます。

また、わからないところの捕捉はこちらがよかったです。

参考：
https://eman-physics.net/quantum/creat_op.html

固有エネルギーの離散化

どうやら古典の世界では、電磁波は電場と磁場のエネルギーを持ってるらしいです。高校でこんなの学んだかどうか完全に忘れてしまいました。位置qと運動量pを用いて、
$$
E = \frac{1}{2}(p^2 + \omega^2 q^2)
$$
このようにかけるようです。これは調和振動子の形になっているので、ばねとして扱うことができそうです。
$$
E = \frac{p^2}{2m} + \frac{1}{2}m\omega^2 x^2
$$
これを量子で見たときに、調和振動子のハミルトニアンは、シュレディンガー方程式より、
$$
\hat{H} = \frac{\hat{p}^2}{2m} + \frac{1}{2}m\omega^2 \hat{x}^2
$$
となり、波動関数を用いて、
$$
\hat{H}\psi(x) = E\psi(x)
$$

を満たすとのことで、さっそくいい感じですね。
次に消滅生成演算子というのが出てきますが、wikipediaによると、「消滅演算子は、状態の粒子の数を1だけ減らす演算子である。生成演算子は、状態の粒子の数を1だけ増やす演算子で、消滅演算子のエルミート共役をとったものである。」とあります。

ということで、消滅演算子 $\hat{a}$ と生成演算子 $\hat{a}^{\dagger}$ を導入します。
$$
\hat{a} = \frac{1}{(2m \hbar \omega)^{1/2}}(m\omega \hat{x} + i\hat{p}_x) \
\hat{a}^{\dagger} = \frac{1}{(2m \hbar \omega)^{1/2}}(m\omega \hat{x} - i\hat{p}_x)
$$
これらを使い、先ほどのハミルトニアンは、
$$
\hat{H} = \hbar \omega \big(\hat{a}^{\dagger}\hat{a} + \frac{1}{2}\big)
$$
と計算できます。手計算していてわかりましたが、 $\hat{x}$ と $\hat{p}$ の計算では交換関係といってかける順番で計算が変わるので計算結果はきちんと計算する必要がありました。
$$
[\hat{a}, \hat{a}^{\dagger}] = \hat{a}\hat{a}^{\dagger} - \hat{a}^{\dagger}\hat{a} = 1
$$
先ほどのハミルトニアンから、固有値方程式 $H\psi_n = E_n \psi_n$ に消滅演算子と生成演算子をそれぞれ作用させると、

\hat{H} \hat{a}^{\dagger}\psi_n = (\hbar \omega + E_n)\hat{a}^{\dagger}\psi_n \\ \hat{H} \hat{a}\psi_n = (-\hbar \omega + E_n)\hat{a}\psi_n

演算子によって、エネルギー固有値がそれぞれ減少・増加しており、エネルギーが量子化されているのが分かりました。これは光子を一つ消滅・生成することと対応して考えられるようで、光子数をベースに光量子計算をするのが光連続量計算のようです。

光子数基底

上記から、光子数を量子状態として量子計算ができそうです。量子状態を光子数の増減で表現して、 $|0>,|1>,|2>,....$ のような状態をつかい、光子数0からスタートし、先ほどの消滅演算子と生成演算子を使って、のようにできます。

\hat{a}^{\dagger}|n\rangle = (n+1)^{1/2}| n+1\rangle \\ \hat{a}|n\rangle = n^{1/2}|n-1\rangle \\ \hat{a}|0\rangle = 0

永井氏の資料から、上記消滅・生成演算子（昇降演算子）を次元を決めて下記のように書いてみることで、

def _downMat(dim, order=1):
    if order == 0:
        A = np.eye(dim)
        return A
    else:
        A = np.zeros([dim, dim])
        for i in np.arange(order, dim):
            A[i - order, i] = np.prod(np.sqrt(np.arange(i, i - order, -1)))
        return A

def _upMat(dim, order=1):
    if order == 0:
        A = np.eye(dim)
        return A
    else:
        A = np.zeros([dim, dim])
        for i in np.arange(order, dim):
            A[i, i - order] = np.prod(np.sqrt(np.arange(i, i - order, -1)))
        return A

import numpy as np

#次元数
d = 4

#消滅演算子
a = _downMat(d)
print("消滅演算子 (a)")
print(a)

#初期状態をたとえば|2>としてみる
state = np.eye(d)[2]
print("\n|2> 状態")
print(state)

#消滅演算子を使って量子状態を変化
print("\na |2> = sqrt(2) |1>")
print(a @ state)

消滅演算子 (a)
[[0.         1.         0.         0.        ]
 [0.         0.         1.41421356 0.        ]
 [0.         0.         0.         1.73205081]
 [0.         0.         0.         0.        ]]

|2> 状態
[0. 0. 1. 0.]

a |2> = sqrt(2) |1>
[0.         1.41421356 0.         0.        ]

消滅演算子によって、光子数基底の|2>状態が|1>状態に遷移しており、係数に光子数の平方根が表れています。
同様に生成演算子も、

#生成演算子、次元は先ほどと同じ4
adag = _upMat(d)
print("生成演算子 (adag)")
print(adag)

#初期状態を|1>としてみる
state = np.eye(d)[1]
print("\n|1> 状態")
print(state)


#生成演算子を使って量子状態を変化
print("\nadag |1> = sqrt(2) |2>")
print(adag @ state)

生成演算子 (adag)
[[0.         0.         0.         0.        ]
 [1.         0.         0.         0.        ]
 [0.         1.41421356 0.         0.        ]
 [0.         0.         1.73205081 0.        ]]

|1> 状態
[0. 1. 0. 0.]

adag |1> = sqrt(2) |2>
[0.         0.         1.41421356 0.        ]

また、
$$
\hat{a}^{\dagger} \hat{a} |n\rangle = n |n\rangle
$$
が成り立つことより、

#初期状態をたとえば|2>としてみる
state = np.eye(d)[2]
print("\n|2> 状態")
print(state)

#演算子をかけてみる
print("\nadag a |2>")
print(adag @ a @ state)


|2> 状態
[0. 0. 1. 0.]

adag a |2>
[0. 0. 2. 0.]

わかりやすいですね。

adag @ a

array([[0., 0., 0., 0.],
       [0., 1., 0., 0.],
       [0., 0., 2., 0.],
       [0., 0., 0., 3.]])

演算子

演算子はユニタリ演算氏となり、光子数の存在確率の総和が1で保存されるように計算をするようです。演算子自体をユニタリにするため（？）、エルミート行列を用いた行列指数関数の形にして状態ベクトルに作用させるようです。

変位ゲート

D(\alpha) = \exp(\alpha \hat{a}^{\dagger} - \alpha^{*}a^{\dagger})

スクイージングゲート

S(\epsilon) = \exp{(\frac{1}{2}\epsilon^{*}\hat{a}^{2}} - \frac{1}{2}\epsilon\hat{a}^{\dagger 2})

\epsilon = re^{2i\phi}

ビームスプリッターゲート

\hat{H} = i\hbar \theta (\hat{a}_{1}^{\dagger}\hat{a}_{2} -\hat{a}_{1} \hat{a}_{2}^{\dagger})

BSgate = \exp(\frac{i}{\hbar}\hat{H})

カーゲート

H = \hbar \frac{\chi}{2} (\hat{a}^{\dagger})^2 \hat{a}^2

K = \exp{(\frac{i}{\hbar}H)} = \exp{(i\frac{\chi}{2} (\hat{a}^{\dagger})^2 \hat{a}^2)} = \exp{(i\frac{\chi}{2} \hat{a}^{\dagger} \hat{n} \hat{a})}

位相平面とウィグナー関数

光子数状態を表現するのに、位置と運動量を軸に取る位相平面を利用し、ウィグナー関数を使って表現できます。
真空状態|0>を位相平面上にウィグナー関数を利用してプロットします。次元は15にし、量子モードは1にします。

from blueqat import photonqat as pq

F = pq.Fock(1, cutoff = 15).run()
(x, p, W) = F.Wigner(0, plot = 'y', xrange = 5.0, prange = 5.0) # plot Wigner function

<Figure size 432x288 with 1 Axes>

明るいほうが存在確率が高くなる表記になっています。スクイーズド状態は、

r = .5
phi = np.pi/2
eps = r * np.exp(2j * phi)

F = pq.Fock(1, cutoff = 15).S(0, eps).run()
(x, p, W) = F.Wigner(0, plot = 'y', xrange = 5.0, prange = 5.0) # plot Wigner function

<Figure size 432x288 with 1 Axes>

ゲートのうちガウシアンゲートは古典的なガウス分布で状態を表せられるようです。ガウス状態を効率的に表現するガウシアンシミュレータを使って同様の計算をしてみます。

G = pq.Gaussian(1).S(0, eps).run()
(x, p, W) = G.Wigner(0, plot = 'y', xrange = 5.0, prange = 5.0) # plot Wigner function

<Figure size 432x288 with 1 Axes>

また、ウィグナー関数へのプロットだけでなく、特定のqumodeの光子数状態を状態ベクトルとして取り出せます。

F = pq.Fock(1, cutoff = 15) # (qumode数, cutoff 数)
F.n_photon(0, 0) # (qumode, photon number state)
F.run()

(x, p, W) = F.Wigner(0, plot = 'y', xrange = 5.0, prange = 5.0) # plot Wigner function

print("Photon number state:\n", F.state)

<Figure size 432x288 with 1 Axes>

Photon number state:
 [1.+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.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j]

次に古典的な状態で表せない非線形なゲートを見てみます。変位ゲートを作用させてからカーゲートを作用させると、

alpha = 2
chi = .2

F = pq.Fock(1, cutoff = 20)
F.D(0, alpha)
F.Kerr(0, chi)
F.run()

x, p, W = F.Wigner(0) # plot

<Figure size 432x288 with 1 Axes>

なんかひん曲がった分布になりました。こうした状態はガウシアンシミュレータでは表現できないのですが、実際に量子コンピュータの実機としても作りづらいようです。

例題：シュレディンガーの猫状態

二つの反対の状態によって構成される量子状態である猫状態を見てみます。
$$
|\mathrm{cat}\rangle \propto |\alpha \rangle \pm |-\alpha\rangle
$$

alpha = 1
mode = 0
F = pq.Fock(1, cutoff = 10)
F.cat(0, alpha, 'e').run()

(x, p, W) = F.Wigner(mode, plot = 'y', xrange = 5.0, prange = 5.0) # plot

<Figure size 432x288 with 1 Axes>

alpha = 1
mode = 0
F = pq.Fock(1, cutoff = 15)
F.cat(0, alpha, 'o').run()

(x, p, W) = F.Wigner(0, plot = 'y', xrange = 5.0, prange = 5.0) # plot

<Figure size 432x288 with 1 Axes>

よくわかりません。そのうち勉強します。
今回は調和振動子から位相平面プロットまでを見ました。光連続量は準位系が多く、量子ゲートが基本行列指数関数の形になっているのでユニタリ行列になっていて、量子ビットみたいな単純な軸の回転みたいなところから入らないので理解が難しいと感じる人もいそうだなという感想を持ちました。しかし、基本は同じなので意外と共通テクが通用したり、両方覚えておいても損はないなと感じました。GBSや連続量QAOAなど光連続量を用いた計算はたびたび開発されてきているので、興味ある方々は、ぜひ皆さんも勉強してみましょう。
以上です。