Strawberry Fieldsで光量子計算をする(その2) インストール〜HOM干渉計

$$
\def\bra#1{\mathinner{\left\langle{#1}\right|}}
\def\ket#1{\mathinner{\left|{#1}\right\rangle}}
\def\braket#1#2{\mathinner{\left\langle{#1}\middle|#2\right\rangle}}
$$

今回はついにStrawberry Fieldsをインストールします。
次に、光の量子的性質がシンプルに確認できる、HOM干渉計について少し解説します。
最後に、Strawberry Fields でHOM干渉計を実装します。

ところで、Strawberry Fields は GUIパッケージも用意されてるようです(僕はまだノータッチですが)。
とりあえず、GUIは使いません。

余談1：光量子計算用のコードには Blackbird quantum programming language という名前が別に付いてるので、最初は Strawberry Fields は GUI環境を指すと思ってました。ただしStrawberry Fields の Tutorial で普通にコーディングをしているので、Strawberry Fields と Black Birds で名前が別れてる理由がよくわからない。

余談2：改めて言うほどの事じゃないかもですけど、XANADU はソフトウェア郡の名前にBeatlesの曲名を使っててクールですね。

Strawberry Fieldsのインストール

Strawberry Fields はPython 上で動くので、基本的には Python を入れて公式Tutorialで挙がってるパッケージ郡を pip install すればOKです。
Tensorflow のバージョンだけ気をつけましょう。

公式ではAnaconda環境が推奨されています。
執筆時点で私は以下の環境です。
- MacBook Pro 2017 (MacOS Mojave 10.14.1)
- Python 3.6.5
- Anaconda3-5.1.0
- numpy 1.14.2
- SciPy 1.1.0
- NetworkX 2.2
- tensorflow 1.6.0
- StrawberryFields 0.9.0

また、TensorFlow のGPU版をインストールしていれば strawberryfields-gpu を入れることでGPUのサポートが受けられるようです。

HOM干渉計

HOM干渉計はとってもシンプルです。

必要なコンポーネントは以下の通り。
- 互いに量子相関可能な2つの単一光子(光子数状態：$\ket{1}$)
- Beam Splitter (BS)
- 光子1個分のエネルギーを判別可能な光子検出器

図で書くと以下のような感じです。

BSに入力1と入力2から光が入射し、それぞれ一定の割合で出力1, 出力2へ反射or透過します。
反射と透過の割合は$\eta$で決まり、$\eta=1/2$ とすると50:50です。
ここでは入射光を生成演算子$\hat{a_{1}}^{\dagger}, \hat{a_{2}}^{\dagger}$で表し、BSが生成演算子に対しユニタリ変換$U_{BS}$として作用します。

U_{BS}=
\begin{bmatrix}
\sqrt{1-\eta} & \sqrt{\eta}
\\ \sqrt{\eta} & -\sqrt{1-\eta}
\end{bmatrix}

それぞれの出力から出てきた光は光子検出器で光子数をカウントされます。
古典的に考えると、仮にBSの反射率が1/2なら出力する光子数は(出力1,出力2)=(2,0), (1,1), (0,2)の3通りで確率はそれぞれ1/4, 1/2, 1/4となるでしょう。

量子力学ではどうなるか。
入力1, 2の状態$\ket{\psi_{in}}_{12}$をもう少し詳しく見ます。(この記事のおかげでケットがきれいに書けました！)
入力1, 2の真空場(光子が存在しない)状態を$\ket{0}_{12}$として表すと、

$$
\ket{\psi_{in}}_{12} = \hat{a_{1}}^{\dagger} \hat{a_{2}}^{\dagger} \ket{0}_{12}
$$

となります。
ここで生成演算子の性質として、光子数状態に作用させると$\hat{a}^{\dagger}\ket{n}=\sqrt{n+1}\ket{n+1}$ であり、$n=0$(真空場)の場合 $\hat{a}^{\dagger}\ket{0}=\ket{1}$ であることを認識すると、$\ket{\psi_{in}}_{12}$ は入力1, 2それぞれが1 光子状態 $\ket{1}$ となっていることがイメージできると思います。

BSを通ることで$\hat{a_{1}}^{\dagger}, \hat{a_{2}}^{\dagger}$ が$U_{BS}$によって変換されるので、簡単のため $\eta=1/2$ とすると出力状態 $\ket{\psi_{out}}_{12}$は

\begin{align}
\ket{\psi_{out}}_{12} &= U_{BS}\ (\hat{a_{1}}^{\dagger} \hat{a_{2}}^{\dagger} \ket{0}_{12}) \\
&=\sqrt{\frac{1}{2}}( \hat{a_{1}}^{\dagger}+\hat{a_{2}}^{\dagger}) \sqrt{\frac{1}{2}}( \hat{a_{1}}^{\dagger}-\hat{a_{2}}^{\dagger}) \ket{0}_{12} \\
&=\frac{1}{2}( \hat{a_{1}}^{\dagger}\hat{a_{1}}^{\dagger}-\hat{a_{1}}^{\dagger}\hat{a_{2}}^{\dagger}+ \hat{a_{1}}^{\dagger}\hat{a_{2}}^{\dagger} - \hat{a_{2}}^{\dagger}\hat{a_{2}}^{\dagger})\ \ket{0}_{12}
\end{align}

ここで、2つの入力光子 $\hat{a_{1}}^{\dagger}$ と $\hat{a_{2}}^{\dagger}$ が区別不可能(;偏光や波長, その他測定可能な量が等しい)な場合、右辺の第二項と第三項が打ち消し合い、以下の式のようになります。

\begin{align}
\ket{\psi_{out}}_{12}&=\frac{1}{2}( \hat{a_{1}}^{\dagger}\hat{a_{1}}^{\dagger}  - \hat{a_{2}}^{\dagger}\hat{a_{2}}^{\dagger})\ \ket{0}_{12} \\
&=\frac{1}{\sqrt{2}}(\ket{2}_1 - \ket{2}_2)
\end{align}

これが意味するのは、「必ず出力1, 出力2のいずれかに光子が2個存在する」です。
出力1, 出力2に光子が1個ずつ存在する解は打ち消されて無くなりました。
実験的には、光子検出器1と光子検出器2で同時に光子が検出されないことを確認します。

以上が、光において量子力学が古典の直感に反するシンプルな例です。
Strawberry Fieldsで再現できるか確かめてみましょう。

HOM干渉計実験のStrawberry Fieldsによる実装

Jupiter Notebookからそのまま持ってきました。

import strawberryfields as sf
from strawberryfields.ops import *

# Engineを初期化する, 引数は用意するqmodesの数
eng, q = sf.Engine(2)

# Engineで動かす量子回路を定義
with eng:
    # 入力状態を定義、qumodesを初期化
    Fock(1)  | q[0]
    Fock(1)  | q[1]

    # Beam Splitterを通す
    BSgate() | (q[0], q[1])

    # 光子数基底でqumodesを測定
    Measure  | q[0]
    Measure  | q[1]

# 量子回路を実行し、計算後の状態をstateに出力
state = eng.run("fock", cutoff_dim=6)

# "Measure"で測定した出力光子数を変数に出力
n_ = q[0].val
m_ = q[1].val

print("Output1 photon number : ", n_)
print("Output2 photon number : ", m_)

Tensor Flowでmodelを定義してから実行するのと似ていますね。
まずEngineを初期化します。qumodes(返り値"q")は2個(要素数2)とします。

eng, q = sf.Engine(2)

次に、入力状態を用意します。Fock は Fock state に由来していて、Photon number stateの別名です。ここでは2つのqumodesを1光子数状態に初期化しています。

    Fock(1)  | q[0]
    Fock(1)  | q[1]

次にBeam Splitterを通し、光子数を測定します。
光子数状態でqumodesを初期化したため、測定は光子数基底で行います。

    # Beam Splitterを通す
    BSgate() | (q[0], q[1])

    # 光子数基底でqumodesを測定
    Measure  | q[0]
    Measure  | q[1]

ここで実行です。
1つ目の引数は計算のバックエンドで、"fock", "gaussian", "tf"から選べます。
"fock"バックエンド上では全ての連続量状態や操作を行えるようです。ただし計算が遅く、誤差が出やすいケースがあるとのこと。
"gaussian"バックエンドは光子数状態やその他非ガウス操作・状態を扱えないですが、"fock"よりも処理が大分軽いです。ガウス操作・状態のみ扱う場合は"gaussian"。
"tf"はまだ使っていないのでよくわかりません...が、Tensorflowを使えそうですね。

2つ目の引数cutoff_dim は、計算に使用する最大の光子数状態です。
一般に完全な連続量量子計算を実現するには無限大の光子数状態が必要なのですが、計算面でも物理実装面でも不可能なので有限個で近似します。
今回の場合は2より大きければ何でも良いかと。

state = eng.run("fock", cutoff_dim=6)

実行すると以下のようになると思います。
必ずOutput1かOutput2のどちらかが光子数2となるはずです。

Output1 photon number :  2
Output2 photon number :  0

確認用に、上記の処理を任意の回数(ここでは1000回)繰り返すようにしてみました。

import strawberryfields as sf
from strawberryfields.ops import *

case1 = 0
case2 = 0

for i in range(1000):
    eng, q = sf.Engine(2)
    with eng:
        Fock(1)  | q[0]
        Fock(1)  | q[1]

        BSgate() | (q[0], q[1])
        Measure  | q[0]
        Measure  | q[1]

    state = eng.run("fock", cutoff_dim=6)
    n_ = q[0].val
    m_ = q[1].val
    if n_ == 2:
        case1 += 1
    elif m_ == 2:
        case2 += 1

print('The case that...')
print('mode 1 has 2 photon : ', case1)
print('mode 2 has 2 photon : ', case2)

print('Sum : ', case1 + case2)

実行結果の一例です。
確率1で、どちらかのモードが光子を2つ持っています。

The case that...
mode 1 has 2 photon :  485
mode 2 has 2 photon :  515
Sum :  1000

まとめ

以上で、駆け足になりましたがHOM干渉計で見られる光子の量子的な振る舞いをStrawberry Fields実装によって確認しました。
各関数や使用可能な状態、ゲートの詳細は公式Tutorialで解説されていますのでぜひ参照ください。
次は量子テレポーテーションかな...。
ご拝読ありがとうございました！

参考にした文献

・arXiv:1711.00080 [quant-ph] "Hong-Ou-Mandel Interference"
(なぜか2017年にとっても気合の入ったHOM干渉計のまとめが出ていた)