量子情報理論の基本：量子チャネル

$$
\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}}
$$

はじめに

周囲の環境と相互作用する開いた量子系は、環境からの雑音（ノイズ）の影響を受けながら、一般に非ユニタリで時間発展していきます。前回の記事で説明した通り、この状態変化は「全体系に対する部分トレース」「Kraus表現」「CPTPマップ」といった形で理論的に記述できるのでした。現実の量子系への応用を考える際には、単なる抽象理論だけでなく、雑音の影響で量子系自体がどのように壊れていくのかをパターン化・モデル化して理解していくことが重要です。通信技術の分野で入力と出力の間をつなぐ通信路＝雑音過程のことを「チャネル」と表現したりしますが、そのアナロジーで量子系のこの状態変化のことを「量子チャネル」と言うようです。今回の記事では、その「量子チャネル」の典型的なパターンをいくつか取り上げて、その特徴を記述してみたいと思います。さらに、量子計算シミュレータqlazyで、そのうちの一つを再現してみます。

参考にさせていただいたのは、以下の文献です。

1. ニールセン、チャン「量子コンピュータと量子通信（２）」オーム社（2005年）
2. ニールセン、チャン「量子コンピュータと量子通信（３）」オーム社（2005年）
3. 石坂、小川、河内、木村、林「量子情報科学入門」共立出版（2012年）
4. 富田「量子情報工学」森北出版（2017年）

ブロッホ表示

量子チャネルの特徴を理解するために、密度演算子の「ブロッホ表示」を知っておく必要があるので、その説明からはじめます。１量子ビットの量子状態は、３次元空間上の単位球＝ブロッホ球の球面上の１点で表現することができます。球面の北極を$\ket{0}$、南極を$\ket{1}$とし、Z軸とのなす角を$\theta$、XY平面への射影とX軸とのなす角を$\phi$とすると、任意の量子状態（純粋状態）は、

\ket{\psi} = \cos\frac{\theta}{2} \ket{0} + e^{i\phi} sin\frac{\theta}{2} \ket{1}  \tag{1}

と表すことができます。これに対応した密度演算子は、

\rho = \frac{1}{2} (\sigma_0 + u_1 \sigma_1 + u_2 \sigma_2 + u_3 \sigma_3)  \tag{2}

と書くことができ、この表現のことを密度演算子の「ブロッホ表示」と呼び、ベクトル$(u_1,u_2,u_3)$を「ブロッホ・ベクトル」と呼びます。ここで、$\sigma_0,\sigma_1,\sigma_2,\sigma_3$はパウリ行列（$\sigma_0$は単位行列）を表しており、$u_1,u_2,u_3$は、

\begin{align}
u_1 &= \sin\theta \cos\phi \\
u_2 &= \sin\theta \sin\phi \\
u_3 &= \cos\theta  \tag{3}
\end{align}

です¹。つまり、純粋状態のブロッホ・ベクトル$(u_1,u_2,u_3)$は、ブロッホ球面上の点であり、これは、量子状態をブロッホ球上で表現したときの座標値に等しいです。

ブロッホ表示は、混合状態にも拡張できます。$\{\sigma_i\}$には、$Tr(\sigma_{i} \sigma_{j})=\delta_{ij}$という関係があるので、$u_1,u_2,u_3$は$\rho$から、

u_i = Tr(\sigma_{i} \rho)  \tag{4}

のように求めることができます。混合状態は純粋状態の線形和（正確には係数の和が１になるので「凸結合」）として表現されるので、式(4)で計算されたブロッホ・ベクトルは、もはや球面上には存在しません。この後、見ていきますが、量子チャンネルによる密度演算子の変化は、このブロッホ・ベクトルの変化として視覚化して理解することができます。

量子チャネルの典型例

それでは本題に入ります。量子チャネルの典型的な例として、

• ビット反転
• 位相反転
• 分極解消
• 振幅ダンピング
• 位相ダンピング

を取り上げ、ブロッホ・ベクトルの変化のパターンや物理的なイメージについて見ていきます。

ビット反転

密度演算子が、確率$p$で以下のように変化する場合を考えます。

\Gamma(\rho) = (1-p) \rho + p \sigma_1 \rho \sigma_1   \tag{5}

右辺第２項の$\sigma_1 \rho \sigma_1$は、パウリのX行列を密度演算子に適用する（量子回路でいうと、１量子ビットにXゲートを演算する）ことに相当しますので、全体として、確率$p$で「ビット反転(bit flip)」する量子チャネルを表しています。例えば、量子コンピュータの１つの量子ビットが、何らかの外部からの雑音によって、確率的にビット反転するようなことをイメージしておけば良いです。

式(5)を見ながらKraus表現を思い出していただくと、この量子チャネルに対するKraus演算子が見えてきます。書き出すと、

\begin{align}
M_0 &= \sqrt{1-p} I \\
M_1 &= \sqrt{p} \sigma_1  \tag{6}
\end{align}

です。

さて、先程説明したブロッホ・ベクトルは、この量子チャネルでどう変化するでしょうか。式(4)を使えば、簡単に計算できます。変化する前後のブロッホ・ベクトルを各々$(u_1,u_2,u_3),\space (u_1^{\prime},u_2^{\prime},u_3^{\prime})$とすると、

\begin{align}
u_{1}^{\prime} &= Tr((1-p) \sigma_{1} \rho + p \rho \sigma_{1}) \\
&= Tr((1-p) \sigma_{1} \rho + p \sigma_{1} \rho) \\
&= Tr(\sigma_{1} \rho) \\
&= u_1 \\
u_{2}^{\prime} &= Tr((1-p) \sigma_{2} \rho + p \sigma_{2} \sigma_{1} \rho \sigma_{1}) \\
&= Tr((1-p) \sigma_{2} \rho - p \sigma_{1} \sigma_{2} \rho \sigma_{1}) \\
&= (1-p) Tr(\sigma_{2} \rho) - p Tr(\sigma_{1} \sigma_{2} \rho \sigma_{1}) \\
&= (1-p) u_2 - p u_2 \\
&= (1-2p) u_2 \\
u_{3}^{\prime} &= (1-2p) u_3  \tag{7}
\end{align}

となります²。ブロッホ・ベクトルのX成分は変わらず、Y,Z成分が$(1-2p)$倍になります。ということは、ブロッホ球面上にあった純粋状態は、ビット反転の量子チャネルによって、Y,Z方向に縮小した楕円体上に移動することになります。

位相反転

密度演算子が、確率$p$で以下のように変化する場合を考えます。

\Gamma(\rho) = (1-p) \rho + p \sigma_3 \rho \sigma_3   \tag{8}

先程の$\sigma_1$が$\sigma_3$に変わりました。$\sigma_3$はパウリのZ行列なので、全体として、確率$p$で「位相反転(phase flip)」する量子チャネルを表しています。先程と同様に、量子コンピュータの１つの量子ビットがノイズの影響で、確率的に位相を変えるというイメージをもっていただいて良いです。

式(8)より、Kraus演算子は、

\begin{align}
M_0 &= \sqrt{1-p} I \\
M_1 &= \sqrt{p} \sigma_3  \tag{9}
\end{align}

です。ブロッホ・ベクトルの変化は、結果だけを示すと（先程と同様の計算なので省略します）、

\begin{align}
u_1^{\prime} &= (1-2p) u_1 \\
u_2^{\prime} &= (1-2p) u_2 \\
u_3^{\prime} &= u_3  \tag{10}
\end{align}

となります。この量子チャネルによって、ブロッホ球面上にあった純粋状態は、X,Y方向に縮小した楕円体上に移動することになります。

分極解消

密度演算子が、確率$p$で以下のように変化する場合を考えます。

\Gamma(\rho) = p \frac{I}{2}+ (1-p) \rho   \tag{11}

これは確率$p$で完全な混合状態である$I/2$に変化することに相当しています³。この量子チャネルは「分極解消(depolarizing)」と呼ばれています。なぜ「分極」の「解消」かというと、ブロッホ球上で特定の方向を向いている純粋状態というのは、例えば、すべての電子のスピンが特定の方向を向いた状態です（当たり前ですが）。つまり、この物理系を外から見ると「分極」しているように見えます。これが完全な混合状態になるということは、全体としてスピンがどっちを向いているかわからない状態（＝完全にランダムにあっち向いたり、こっち向いたりしている状態）です。つまり、「分極」が「解消」している状態になります。

さて、改めて、式(10)の形を見ていただくと、ビット反転や位相反転の場合と違って、Kraus表現の格好になっていません。が、密度演算子に対する恒等式⁴、

\frac{1}{2} (\rho + \sigma_{1} \rho \sigma_{1} + \sigma_{2} \rho \sigma_{2} + \sigma_{3} \rho \sigma_{3}) = I \tag{12}

を使うことで、以下のように式変形できます（$q=3p/4$とおきました）。

\Gamma(\rho) = (1-q) \rho + \frac{q}{3} (\sigma_{1} \rho \sigma_{1} + \sigma_{2} \rho \sigma_{2} + \sigma_{3} \rho \sigma_{3})   \tag{13}

ということで、Kraus演算子が見えてきました。

M_0 = \sqrt{1-q} I, \space M_1 = \sqrt{\frac{q}{3}} \sigma_{1}, \space M_2 = \sqrt{\frac{q}{3}} \sigma_{2}, \space M_3 = \sqrt{\frac{q}{3}} \sigma_{3}  \tag{14}

です。ブロッホ・ベクトルの変化は、これも結果だけ示すと、

\begin{align}
u_{1}^{\prime} &= (1-p) u_{1} \\
u_{2}^{\prime} &= (1-p) u_{2} \\
u_{3}^{\prime} &= (1-p) u_{3}  \tag{15}
\end{align}

となります。したがって、この量子チャネルによって、ブロッホ球面上にあった純粋状態は、全体に縮小した球面上に移動することになります。

振幅ダンピング

次に、もう少し、現実の物理現象として想起できそうな量子チャネルについて見ていきます。まず、「振幅ダンピング(amplitude dumping)」です。

例えば、フォトンを自然放出する原子や、スピン多体系が環境と相互作用しながら平衡状態に近づく現象や、干渉計や共振器を通して散乱や減衰を受けるフォトンなど、何らかの相互作用によって、考えている量子系からエネルギーが散逸する過程を具体的にイメージすれば良いです。

一例として、フォトンが何かに当たって確率的に散乱するプロセスを考えてみます。具体的には「ビームスプリッタ」と呼ばれる光学素子にフォトンが入射することを考えます。

ということで、振幅ダンピングの説明の前に、ちょっとだけ「ビームスプリッタ」とは何か、について説明しておきます。２つのプリズムを合わせて、その境界面をハーフミラーにしたような光学素子がビームスプリッタです。境界面に斜め45度から入射した光(A)は、ある確率で斜め45度の角度で反射するか、あるいはそのまま透過します。一方、境界面に対して(A)とは反対の斜め45度から入射する別の光(B)を考えると、同じようにある確率で反射したり透過したりします。入射する２つの光の強度を各々$A,B$とすると、その光たちは、境界面の反射率（透過率）に応じて各々混じり合って出射することになります。その強度を$A^{\prime},B^{\prime}$とすると、以下のような関係が成り立ちます（$\theta$は、境界面の反射率（透過率）から決まるパラメータです）。要するに、Aが減った分だけBが増え、あるいはAが増えた分だけBが減る関係の直交変換（A,Bをベクトル成分とした回転）です。

\begin{align}
A^{\prime} &= A \cos\theta + B \sin\theta \\
B^{\prime} &= -A \sin\theta + B \cos\theta   \tag{16}
\end{align}

ここまでは、古典光学におけるビームスプリッタの説明です。が、ここは量子チャネルを説明するところなので、大量の光が一度に入射するのではなく、量子効果が現れるように、ゼロ個または１個のフォトンが入射するようなプロセスを考えてみます。したがって、ここで光の量子化についての知識が必要になるのですが、一から説明すると長くなってしまうので、天下ります。光の量子化は、生成演算子$a^{\dagger}$、消滅演算子$a$によって定義される数演算子$a^{\dagger} a$の固有状態$\ket{n}$を使って、記述することができます。結果、以下が成り立ちます。

\begin{align}
a^{\dagger} a \ket{n} &= n \ket{n} \\
a^{\dagger} \ket{n} &= \sqrt{n+1} \ket{n+1} \\
a \ket{n} &= \sqrt{n} \ket{n-1} \\
[a,a^{\dagger}] &= 1  \tag{17}
\end{align}

ここで、$\ket{n}$の中の$n$はフォトンの数を表していると思ってください（電子スピンの上向き、下向きのイメージではないのでご注意）。とすると、確かに$a^{\dagger} a$は数演算子の名にふさわしい動きをしていますし、生成演算子も消滅演算子も、確かにその名の通り、フォトンを増やしたり減らしたりする動きをしていることがわかります。

それでは、前置きが長くなりましたが、振幅ダンピングの説明に移ります。先程のビームスプリッタでは光が通る経路としてAとBを考えて、両方から光が入射するという想定でした。いまゼロ個あるいは１個のフォトンを考え、そのエネルギー損失を評価してみたいので、経路Aを環境系(S)、経路Bを注目系(S)として、この２つの経路にゼロ個または１個のフォトンが入射する想定にします。そうすると、入射のパターンは、$\ket{0}_E \ket{0}_S, \ket{1}_E \ket{0}_S, \ket{0}_E \ket{1}_S$の３パターンになります（$\ket{1}_E \ket{1}_S$はフォトンが２個のパターンになるので考えません）。古典光学のビームスプリッタの変換式から、この３パターンは各々、以下のように変化します。

\begin{align}
\ket{0}_E \ket{0}_S &\rightarrow \ket{0}_E \ket{0}_S \\
\ket{1}_E \ket{0}_S &\rightarrow \cos\theta \ket{1}_E \ket{0}_S - \sin\theta \ket{0}_E \ket{1}_S \\
\ket{0}_E \ket{1}_S &\rightarrow \sin\theta \ket{1}_E \ket{0}_S + \cos\theta \ket{0}_E \ket{1}_S   \tag{18}
\end{align}

この変換を演算子$U$と表すと、

U = \ket{0}_E \ket{0}_S \bra{0}_E \bra{0}_S + (\cos\theta \ket{1}_E \ket{0}_S - \sin\theta \ket{0}_E \ket{1}_S) \bra{1}_E \bra{0}_S + (\sin\theta \ket{1}_E \ket{0}_S + \cos\theta \ket{0}_E \ket{1}_S) \bra{0}_E \bra{1}_S  \tag{19}

と書けます。計算すればわかるように、これはユニタリ演算子です⁵。ということで、全体系に対するユニタリ演算子がわかったので、Kraus演算子は、定義から、以下のように求めることができます。

\begin{align}
M_0 &= \bra{0}_E U \ket{0}_E = \ket{0}_S \bra{0}_S + \cos\theta \ket{1}_S \bra{1}_S \\
M_1 &= \bra{1}_E U \ket{0}_E = \sin\theta \ket{0}_S \bra{1}_S  \tag{20}
\end{align}

ここで、$\gamma = \sin^{2} \theta$とおくと、$\gamma$はフォトンを失う確率を表すパラメータとなり、上のKraus演算子は、以下のように書き直せます。

\begin{align}
M_0 &= \ket{0}_S \bra{0}_S + \sqrt{1-\gamma} \ket{1}_S \bra{1}_S \\
M_1 &= \sqrt{\gamma} \ket{0}_S \bra{1}_S  \tag{21}
\end{align}

これより、Kraus表現は、

\begin{align}
\Gamma(\rho) &= M_{0} \rho M_{0}^{\dagger} + M_{1} \rho M_{1}^{\dagger} \\
&= (\ket{0} \bra{0} + \sqrt{1-\gamma} \ket{1} \bra{1}) \rho (\ket{0} \bra{0} + \sqrt{1-\gamma} \ket{1} \bra{1}) + \sqrt{\gamma} \ket{0} \bra{1} \rho \ket{1} \bra{0}  \tag{22}
\end{align}

となります（面倒なので注目系を表す添字Sは省略しました）。式(21)(22)をじっと見ていただくと、なんとなくわかってきますが、$M_0$は、$\ket{0}$が来たらそのまま出して、$\ket{1}$が来たらその確率を減らす効果を表していて、一方、$M_1$は、$\ket{1}$が来たら$\ket{0}$に遷移させる効果を表しています。全体として、入ってきたフォトンがそのまま出ていく確率を減らす、つまり減衰させる効果を表す量子チャネルになっています。

ブロッホ・ベクトルの変化を見てみると、

\begin{align}
u_1 &= \sqrt{1-\gamma} u_1 \\
u_2 &= \sqrt{1-\gamma} u_2 \\
u_3 &= \gamma + \sqrt{1-\gamma} u_3 \tag{23}
\end{align}

と計算できます。これより、ブロッホ球は、$\gamma$が大きくなるに従い、X軸、Y軸方向に縮小していき、同時にZ軸成分は1に近づいていきます。すなわち、$\ket{0}$＝北極に向かって、どんどん縮小していくパンケーキのイメージですね。

位相ダンピング

現実の物理現象として想起できそうな量子チャネルのもう一つの例は、「位相ダンピング(phase dumping)」です。環境との相互作用によって量子状態の位相が変化する過程なので、これは量子系固有の雑音過程と言えます。一般に、量子状態をエネルギー固有状態の重ね合わせで表現したとき、各々の固有状態はその固有エネルギーの値に比例する形で位相変化していきます。雑音によってその固有状態間の相対位相がずれていくというのが位相ダンピングのイメージです。

先程の振幅ダンピングと同様に、ビームスプリッタに入射するフォトンの散乱過程として位相ダンピングを記述してみます。ただし、今回はフォトンは減衰しません。つまり、散乱によってフォトンがなくなることはないので、ハミルトニアンに

\ket{1}_S \bra{1}_S

という因子がないといけません。また、この散乱によって環境系でフォトンが生成される過程を表す

\ket{1}_E \bra{0}_E

という因子もありえます⁶。そうだとすると、さらに、全体をエルミートにするために

\ket{0}_E \bra{1}_S

という因子（逆に消滅を表す因子）もハミルトニアンには必要です。この条件を満たす、最も簡単なものは、

H = \chi \ket{1}_S \bra{1}_S (\ket{1}_E \bra{0}_E + \ket{0}_E \bra{1}_E)  \tag{24}

です⁷。

このハミルトニアンで表される系の時間発展は、以下のユニタリ演算子で記述できます。

U = exp(iHt) = exp(i \chi t \ket{1}_S \bra{1}_S (\ket{1}_E \bra{0}_E + \ket{0}_E \bra{1}_E))  \tag{25}

指数関数のテイラー展開を使いながら、地道に計算すると、

U = \ket{0}_S \bra{0}_S (\ket{0}_E \bra{0}_E + \ket{1}_E \bra{1}_E) + \cos(\chi t) \ket{1}_S \bra{1}_S (\ket{1}_E \bra{0}_E + \ket{0}_E \bra{1}_E) + i \sin(\chi t) \ket{1}_S \bra{1}_S (\ket{1}_E \ket{0}_E + \ket{0}_E \bra{1}_E)  \tag{26}

となります⁸。$\gamma = \sin^{2}(\chi t)$とおくと、

U = \ket{0}_S \bra{0}_S (\ket{0}_E \bra{0}_E + \ket{1}_E \bra{1}_E) + \sqrt{1-\gamma} \ket{1}_S \bra{1}_S (\ket{1}_E \bra{0}_E + \ket{0}_E \bra{1}_E) + i \sqrt{\gamma} \ket{1}_S \bra{1}_S (\ket{1}_E \ket{0}_E + \ket{0}_E \bra{1}_E)  \tag{27}

となり、Kraus演算子は、

\begin{align}
M_0 &= \bra{0}_E U \ket{0}_E = \ket{0} \bra{0} + \sqrt{1-\gamma} \ket{1} \bra{1} \\
M_1 &= \bra{1}_E U \ket{0}_E = i \sqrt{\gamma} \ket{1} \bra{1}  \tag{28}
\end{align}

となります（面倒なので添字Sは省略しました）。したがって、位相ダンピングの量子チャネルは、

\Gamma(\rho) = (\ket{0} \bra{0} + \sqrt{1-\gamma} \ket{1} \bra{1}) \rho (\ket{0} \bra{0} + \sqrt{1-\gamma} \ket{1} \bra{1}) + \gamma \ket{1} \bra{1} \rho \ket{1} \bra{1} \tag{29}

と書くことができます。

さて、ここで面白い補足をしておきます。実は、この量子チャネルは、位相反転とまったく同じ効果を表しています。位相反転の式を再度掲載します。

\Gamma(\rho) = (1-p) \rho + p \sigma_3 \rho \sigma_3   \tag{7}

これと、式(29)は同じです。一見違って見えますが、ちょっと式変形してみます。まず、式(29)です。

\begin{align}
\Gamma(\rho) &= \ket{0}\bra{0}\rho\ket{0}\bra{0} + \sqrt{1-\gamma} \ket{0}\bra{0}\rho\ket{1}\bra{1} + \sqrt{1-\gamma} \ket{1}\bra{1}\rho\ket{0}\bra{0} + (1-\gamma) \ket{1}\bra{1}\rho\ket{1}\bra{1} + \gamma \ket{1}\bra{1}\rho\ket{1}\bra{1} \\
&= \ket{0}\bra{0}\rho\ket{0}\bra{0} + \sqrt{1-\gamma} \ket{0}\bra{0}\rho\ket{1}\bra{1} + \sqrt{1-\gamma} \ket{1}\bra{1}\rho\ket{0}\bra{0} + \ket{1}\bra{1}\rho\ket{1}\bra{1}  \tag{30}
\end{align}

次に、式(7)です。途中で、$\sum_{i} \ket{i}\bra{i} = 1$を使います。

\begin{align}
\Gamma(\rho) &= (1-p) \rho + p (\ket{0}\bra{0} - \ket{1}\bra{1}) \rho (\ket{0}\bra{0} - \ket{1}\bra{1}) \\
&= (1-p) (\ket{0}\bra{0} + \ket{1}\bra{1}) \rho (\ket{0}\bra{0} + \ket{1}\bra{1}) + p (\ket{0}\bra{0} - \ket{1}\bra{1}) \rho (\ket{0}\bra{0} - \ket{1}\bra{1}) \\
&= \ket{0}\bra{0} \rho \ket{0}\bra{0} + (1-2p) \ket{0}\bra{0} \rho \ket{1}\bra{1} + (1-2p) \ket{1}\bra{1} \rho \ket{0}\bra{0} + \ket{1}\bra{1} \rho \ket{1}\bra{1}  \tag{31}
\end{align}

式(30)と式(31)を比較すると、

p = \frac{(1-\sqrt{1-\gamma})}{2}  \tag{32}

とおけば一致することがわかります。したがって、ブロッホ球の変化も位相反転の場合とまったく同じになります。

歴史的に位相ダンピングは、何らかの散乱過程による連続的な位相誤差の結果として生じると理解されていたので、確率的に起きる離散的な位相反転の結果と同等というこの知見は、ちょっとした発見だったようで、これによって始めて量子誤り訂正が発展しました（とニールセン、チャンに書いてありました）⁹。

シミュレータによる量子チャネルの再現

それでは、いままで説明した中から、振幅ダンピングを取り上げて、その効果をqlazyで確認してみます。

振幅ダンピング

最初の純粋状態を$(u_1,u_2,u_3)=(1/\sqrt{3},1/\sqrt{3},1/\sqrt{3})$、減衰率を表すパラメータ$\gamma$をとりあえず0.5として、この量子チャネルの前後で、密度演算子と密度演算子の2乗トレース（純粋状態か混合状態かの確認）とブロッホ・ベクトルの値を出力するようにしました。また、$a^{\dagger}a$の期待値、つまりエネルギーの期待値（正確にはフォトン数の期待値ですが光の周波数を固定すればエネルギーに比例しているので代用）も合わせて計算して表示するようにしました。

全体のPythonコードを以下に示します。

【2021.9.5追記】qlazy最新版でのソースコードはここに置いてあります。

import math
import numpy as np
from qlazypy import DensOp

# Pauli Matrix

Sigma_0 = np.eye(2)
Sigma_1 = np.array([[0,1],[1,0]])
Sigma_2 = np.array([[0,-1j],[1j,0]])
Sigma_3 = np.array([[1,0],[0,-1]])

def get_coordinate(densop=None):

    u_1 = densop.expect(matrix=Sigma_1)
    u_2 = densop.expect(matrix=Sigma_2)
    u_3 = densop.expect(matrix=Sigma_3)

    return (u_1,u_2,u_3)
    
def make_densop_matrix(u_1,u_2,u_3):

    matrix = (Sigma_0+u_1*Sigma_1+u_2*Sigma_2+u_3*Sigma_3) / 2.0

    return matrix

def make_kraus(gamma=0.0):

    transmit = math.sqrt(1.0-gamma)
    reflect = math.sqrt(gamma)
    
    kraus = []
    kraus.append(np.array([[1,0],[0,0]])+transmit*np.array([[0,0],[0,1]]))
    kraus.append(reflect*np.array([[0,1],[0,0]]))

    return kraus

def make_hamiltonian():

    return np.array([[0,0],[0,1]])
    
if __name__ == '__main__':

    print("== parameter ==")

    gamma = 0.5
    H = make_hamiltonian()
    print("gamma =", gamma)
    
    print("== initial density operator ==")
    
    u_1 = math.sqrt(1/3)
    u_2 = math.sqrt(1/3)
    u_3 = math.sqrt(1/3)
    D = make_densop_matrix(u_1,u_2,u_3)
    de = DensOp(matrix=D)

    de.show()
    print("square trace =",de.sqtrace())
    print("(u_1,u_2,u_3) = ({0:.3f},{1:.3f},{2:.3f})".format(u_1,u_2,u_3))
    print("expect value of energy =", de.expect(matrix=H))

    [M_0,M_1] = make_kraus(gamma=gamma)

    print("== finail density operator ==")

    de.instrument(kraus=[M_0,M_1])

    de.show()
    print("square trace =", de.sqtrace())
    (u_1,u_2,u_3) = get_coordinate(densop=de)
    print("(u_1,u_2,u_3) = ({0:.3f},{1:.3f},{2:.3f})".format(u_1,u_2,u_3))
    print("expect value of energy =", de.expect(matrix=H))

    de.free()

何をやっているか、ざっくり説明します。

print("== parameter ==")

gamma = 0.5
H = make_hamiltonian()
print("gamma =", gamma)

$\gamma$の値を設定して、フォトン数演算子の行列表現を作っています（ハミルトニアンの代用品として）。この演算子の固有状態$\{\ket{n}\} \space (n=0,1)$の基底で$a,a^{\dagger},a^{\dagger}a$を行列表現すると、

a =
\begin{pmatrix}
0 & 1 \\
0 & 0
\end{pmatrix}, \space
a^{\dagger} =
\begin{pmatrix}
0 & 0 \\
1 & 0
\end{pmatrix}, \space
a^{\dagger} a =
\begin{pmatrix}
0 & 0 \\
0 & 1
\end{pmatrix}  \tag{33}

です。make_hamiltonian関数で、$a^{\dagger} a$の行列表現を設定しています。

print("== initial density operator ==")
    
u_1 = math.sqrt(1/3)
u_2 = math.sqrt(1/3)
u_3 = math.sqrt(1/3)
D = make_densop_matrix(u_1,u_2,u_3)
de = DensOp(matrix=D)

ブロッホ・ベクトルの成分を$(1/\sqrt{3},1/\sqrt{3},1/\sqrt{3})$に設定します（何でも良かったのですが適当です）。make_densop_matrix関数でそれに対応した密度演算子の行列表現を計算して、変数Dに格納します。Dから密度演算子クラスのインスタンスdeを作成します。

de.show()
print("square trace =",de.sqtrace())
print("(u_1,u_2,u_3) = ({0:.3f},{1:.3f},{2:.3f})".format(u_1,u_2,u_3))
print("expect value of energy =", de.expect(matrix=H))

「密度演算子の各要素」「２乗トレース」「ブロッホ・ベクトルの成分」「ハミルトニアンの期待値」の初期値を表示します。

[M_0,M_1] = make_kraus(gamma=gamma)

print("== finail density operator ==")

de.instrument(kraus=[M_0,M_1])

make_kraus関数で、設定した$\gamma$に対するKraus演算子をリスト[M_0,M_1]に格納します。そして、insutrumentメソッド¹⁰でKraus表現にします。つまり、量子チャネルに通します。

de.show()
print("square trace =", de.sqtrace())
(u_1,u_2,u_3) = get_coordinate(densop=de)
print("(u_1,u_2,u_3) = ({0:.3f},{1:.3f},{2:.3f})".format(u_1,u_2,u_3))
print("expect value of energy =", de.expect(matrix=H))

量子チャネル通過後の「密度演算子の各要素」「２乗トレース」「ブロッホ・ベクトルの成分」「ハミルトニアンの期待値」を表示します。

実行結果は以下です。まず$\gamma=0.5$の場合です。

== parameter ==
gamma = 0.5
== initial density operator ==
elm[0][0] = +0.7887+0.0000*i : 0.6220 |+++++++
elm[0][1] = +0.2887-0.2887*i : 0.1667 |+++
elm[1][0] = +0.2887+0.2887*i : 0.1667 |+++
elm[1][1] = +0.2113+0.0000*i : 0.0447 |+
square trace = 1.0
(u_1,u_2,u_3) = (0.577,0.577,0.577)
expect value of energy = 0.21132487
== finail density operator ==
elm[0][0] = +0.8943+0.0000*i : 0.7998 |+++++++++
elm[0][1] = +0.2041-0.2041*i : 0.0833 |++
elm[1][0] = +0.2041+0.2041*i : 0.0833 |++
elm[1][1] = +0.1057+0.0000*i : 0.0112 |+
square trace = 0.9776709
(u_1,u_2,u_3) = (0.408,0.408,0.789)
expect value of energy = 0.10566243

ブロッホ・ベクトルを見ると、北極に少し近づいていることがわかります。また、エネルギー期待値は半分になっています。次に、$\gamma=0.8$としてみると、以下のようになりました（最終状態のみ切り出して表示）。

== finail density operator ==
elm[0][0] = +0.9577+0.0000*i : 0.9173 |++++++++++
elm[0][1] = +0.1291-0.1291*i : 0.0333 |+
elm[1][0] = +0.1291+0.1291*i : 0.0333 |+
elm[1][1] = +0.0423+0.0000*i : 0.0018 |+
square trace = 0.98570938
(u_1,u_2,u_3) = (0.258,0.258,0.915)
expect value of energy = 0.04226497

さらに、北極に近づき、エネルギーも減衰しました。最後に$\gamma=1.0$としてみると、以下のようになりました（最終状態のみ切り出して表示）。

== finail density operator ==
elm[0][0] = +1.0000+0.0000*i : 1.0000 |+++++++++++
elm[0][1] = +0.0000+0.0000*i : 0.0000 |
elm[1][0] = +0.0000+0.0000*i : 0.0000 |
elm[1][1] = +0.0000+0.0000*i : 0.0000 |
square trace = 1.0
(u_1,u_2,u_3) = (0.000,0.000,1.000)
expect value of energy = 0.0

というわけで、全エネルギーを使い切ってしまいましたが、ついに北極点に到達しました。おめでとうございます！

おわりに

量子チャネルの物理的イメージも含めて、基本的な理解ができました。この勢いで「量子誤り訂正」に進みたい気がしてきましたが、その前におさえておくべき基本事項（「シュミット分解」とか「純粋化」とか「エントロピー」とか諸々）がいくつか残っているので、次回以降、順に見ていきたいと思います。

以上

1. $\sigma_0=\ket{0}\bra{0}+\ket{1}\bra{1},\sigma_1=\ket{0}\bra{1}+\ket{1}\bra{0},\sigma_2=-i\ket{0}\bra{1}+i\ket{1}\bra{0},\sigma_3=\ket{0}\bra{0}-\ket{1}\bra{1}$に注意しながら、三角関数の半角公式、倍角公式を使えば、式(1)から式(2)は導けます。 ↩

2. パウリ行列の性質：$\sigma_{i} \sigma_{i} = I, \space \{\sigma_{i},\sigma_{j}\}=\sigma_{i} \sigma_{j} + \sigma_{j} \sigma_{i} = 0$、およびトレースの性質：$Tr(aA+bB)=aTr(A)+bTr(B), \space Tr(ABC) = Tr(BCA) = Tr(CAB)$を使いました。 ↩

3. この場合の完全な混合状態は、$\ket{0}$と$\ket{1}$が半々に混ざりあった状態なので、密度演算子の定義から、$I/2$になるのは簡単にわかりますね。 ↩

4. パウリ行列の定義と密度演算子のトレースが１であることを使って、左辺を地道に計算すれば、わかります。 ↩

5. 参考文献4にもユニタリです、と書いてありますが、$\ket{1}_E \ket{1}_S$を入れていないので、ユニタリの次元は３ということになります。２のべきではない次元は、量子情報的にはあまり出てこないと思うのですが、こんな場合もあって良いのでしょうか。わかりませんが、今の場合、特に問題なさそうです。 ↩

6. 無から有が生まれる過程ですね。生成消滅演算子を導入して第２量子化するというのは、要はこういう状況も記述する必要に迫られて編み出された方法と言えます。もともとのシュレーディンガー方程式（第１量子化）では粒子の生成消滅は記述できなかったので、、 ↩

7. 振幅ダンピングの説明のときは、古典光学の現象と対応させる形で散乱過程を行列表現しましたが、位相ダンピングの場合、ハミルトニアンをいきなり人為的に作ってしまっています。位相雑音は量子系特有の現象なので、うまく対応づけられるような古典光学の現象がないからだと思います。なんとなく釈然としないですが、ここは我慢です。 ↩

8. 長くなるので計算過程はバッサリ省略しました。どうしても自分で導出したいという方は、参考文献4を見てください。もう少し式変形のヒントが書いてあります。 ↩

9. ニールセン、チャンの本では、Kraus演算子のユニタリな自由度に基づき両者相互に変換できるということで同等、という説明がなされていました。どんなユニタリ変換？という具体的な式変形が示されおらず、自力導出もちょっと面倒な気がしたので、ここでは、$\Gamma(\rho)$が両者同じになるよね、という説明の仕方にしました。 ↩

10. この演算のことを「CP-instrument」と言ったりするので、それに従いメソッド名としました。 ↩