$$
\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}}
$$
はじめに
前回の記事「量子情報理論の基本：密度演算子」で密度演算子の基本的な概念と性質について説明しました。軽くおさらいすると、量子状態には純粋状態と混合状態というものがあって、それを統一的に便利に記述できるものとして密度演算子が導入されました。時間発展や測定が密度演算子を使って書けましたし、トレースが必ず1になるとか半正定値であるという数学的な性質があったり、また、量子状態のアンサンブルに関してユニタリの自由度があったりしました。さらに、純粋状態と混合状態を判別するために密度演算子の２乗のトレースをとれば良いということもわかりました。
今回は、全体の量子系の「部分系」がどうなっているのかを知る手段として、密度演算子に「部分トレース」という手続きを施せば良いということを説明します。そして、自作の量子計算シミュレータqlazyを使って、部分トレースの動作を確認したいと思います。
参考にさせていただいたのは、以下の記事・文献です。
オブザーバブルの期待値について
部分系と部分トレースの話を進める上で、オブザーバブルの期待値を密度演算子でどう表現するかを知っておく必要があるので、はじめに説明しておきます。量子状態アンサンブル$\{p_i, \ket{\psi_i}\}$におけるオブザーバブル$A$の期待値$<A>$は、密度演算子$\rho$を使って、以下のように書くことができます。
\begin{align}
<A> &= \sum_{i} p_i \bra{\psi_i} A \ket{\psi_i} \\
&= \sum_{i,k} p_i \braket{\psi_i}{k} \bra{k} A \ket{\psi_i} \\
&= \sum_{i,k} p_i \bra{k} A \ket{\psi_i} \braket{\psi_i}{k} \\
&= \sum_{i} p_i Tr(A \ket{\psi_i} \bra{\psi_i}) \\
&= Tr(A \sum_{i} p_i \ket{\psi_i} \bra{\psi_i}) \\
&= Tr(A \rho) \tag{0}
\end{align}
部分トレースによる部分系の記述
それでは、本題に入ります。
まず、系Aで定義されるあるオブザーバブル$T$があったとき、複合系ABから見てこれと同等の効果を与えるオブザーバブルを$\tilde{T}$とします。この$T$と$\tilde{T}$は、どのように関連付けられるべきかをまず考えます。$T$の固有状態が$\ket{t}$で、複合系ABの状態が$\ket{t}\ket{\psi}$であったとします。そうすると、$T$の測定結果は100%の確率で$t$になります。系Aにおける測定に対応した射影演算子を$P_t$とすると、複合系ABでの射影演算子は$P_t \otimes I_B$です。$\tilde{T}$は、この射影演算子に固有値をかけて$t$について和をとったものに等しいので、
\tilde{T} = \sum_{t} t P_t \otimes I_B = T \otimes I_B \tag{1}
と書けます。どうでしょう。複合系ABにおけるオブザーバブル$\tilde{T}$は、部分系Aにおけるオブザーバブル$T$と部分系Bの恒等演算子$I_B$のテンソル積ということで、感覚的にはとても自然な感じがします。
では、複合系ABにおける密度演算子は、部分系Aでどのように表現されるべきでしょうか。任意のオブザーバブル$T$に対して、系Aでも複合系ABでも等しい統計的な結果を与えないといけないので、期待値は等しいはずです。つまり、
Tr(T \rho^A) = Tr(\tilde{T} \rho^{AB}) = Tr((T \otimes I_B) \rho^{AB}) \tag{2}
です。このとき、$\rho^A$が$\rho^{AB}$の「部分トレース」になっていれば、この式が成り立つことを以下で証明します。
【証明】
$\rho^{AB}$は、系Aにおける状態$\ket{a_1}, \ket{a_2}$と系Bにおける状態$\ket{b_1}, \ket{b_2}$を用いて、$\rho^{AB}=\ket{a_1} \bra{a_2} \otimes \ket{b_1} \bra{b_2}$と一般的に書くことができます。このとき、系Bについての部分トレース$Tr_B(\rho^{AB})$は、
Tr_B(\rho^{AB}) = Tr_B(\ket{a_1} \bra{a_2} \otimes \ket{b_1} \bra{b_2}) = \ket{a_1} \bra{a_2} Tr(\ket{b_1} \bra{b_2}) \tag{3}
と定義されます。
式(2)の最右辺は、
\begin{align}
Tr((T \otimes I_B) \rho^{AB}) &= Tr((T \otimes I_B) (\ket{a_1} \bra{a_2} \otimes \ket{b_1} \bra{b_2})) \\
&= Tr(T \ket{a_1} \bra{a_2} \otimes I_B \ket{b_1} \bra{b_2}) \\
&= Tr(T \ket{a_1} \bra{a_2}) \cdot Tr(\ket{b_1} \bra{b_2}) \\
&= Tr(T \ket{a_1} \bra{a_2} Tr(\ket{b_1} \bra{b_2})) \tag{4}
\end{align}
のように展開されます。式(3)の定義を式(4)に代入すると、
Tr((T \otimes I_B) \rho^{AB}) = Tr(T \cdot
Tr_B(\rho^{AB})) \tag{5}
となります。これより、式(2)が成り立つためには、
\rho^A = Tr_B(\rho^{AB}) \tag{6}
が成り立っていれば良いことがわかります。(証明終)
部分系は部分トレースで計算できることがわかったので、次に、具体例を２パターン示します。
例1：積状態に対する部分系
量子系Aの状態$\rho^A$と量子系Bの状態$\rho^B$の積状態$\rho^{AB}$というのは、
\rho^{AB} = \rho^A \otimes \rho^B
と定義されます。この全体の状態$\rho^{AB}$があったときに、部分系Aについての量子状態は、Bについての部分トレースなので、
Tr_B(\rho^{AB}) = Tr_B(\rho^A \otimes \rho^B) = \rho^A Tr(\rho^B) = \rho^A
となり、$\rho^A$という状態が得られました。同様に部分系Bについての量子状態は$\rho^B$になります。というわけで、$\rho^{AB}$が純粋状態で、上のように書ける積状態だったとすると、その部分系は純粋状態になります(かなり自明な感じですね)。
例2：量子エンタングル状態に対する部分系
では、$\rho^{AB}$が、積状態になっていない純粋状態だった場合はどうでしょうか。例えば、２量子ビットのベル状態$(\ket{00} + \ket{11})/\sqrt{2}$の場合、どうなるかということです。密度演算子は、
\begin{align}
\rho^{1,2} &= \frac{\ket{00}+\ket{11}}{\sqrt{2}} \cdot \frac{\bra{00}+\bra{11}}{\sqrt{2}} \\
&= \frac{1}{2} (\ket{00}\bra{00} + \ket{11}\bra{00} + \ket{00}\bra{11} + \ket{11}\bra{11})
\end{align}
となります。部分系2に対する部分トレースをとってみると(つまり部分系1はどうなっているかを見てみると)、
\begin{align}
\rho^1 &= Tr_{2}(\rho^{1,2}) \\
&= \frac{1}{2} (Tr_2(\ket{00}\bra{00}) + Tr_2(\ket{11}\bra{00}) + Tr_2(\ket{00}\bra{11}) + Tr_2(\ket{11}\bra{11})) \\
&= \frac{1}{2} (\ket{0}\bra{0}\braket{0}{0} + \ket{1}\bra{0}\braket{0}{1} + \ket{0}\bra{1}\braket{1}{0} + \ket{1}\bra{1}\braket{1}{1}) \\
&= \frac{1}{2} (\ket{0}\bra{0}+\ket{1}\bra{1}) \\
&= \frac{I}{2}
\end{align}
となります。２乗のトレースをとると、
Tr((\rho^1)^2) = (\frac{1}{2})^2 + (\frac{1}{2})^2 = \frac{1}{2} < 1
ということで、なんと混合状態ですね。純粋状態の部分系が混合状態になってしまうのは、ちょっと不思議な感じがしますが、量子エンタングル状態というものは、そういうものなのです。ここは本日のポイントですので、覚えておきましょう。
部分トレース計算の実装
部分トレースの定義がわかったところで、では、コンピュータプログラムとして、どうやって実装すれば良いでしょうか(シミュレータの機能として組み込みたいので)。$n$量子ビット系の密度演算子は、$2^n$行$2^n$列の複素正方行列で表現できます。これに対して、例えば、a番目とb番目とc番目の量子ビットについての部分トレースを計算したいのです。しばし悩んでしまいました。結局以下のように実装しましたので、参考まで紹介してみます(もっとスマートなアルゴリズムで、こうやればいいみたいなコードサンプルが、どこかにあるような気もするのですが、、)。
考え方は、$\rho$の行列表示、
\bra{i_1,i_2,\cdots,i_n} \rho \ket{j_1,j_2,\cdots,j_n}
のa番目とb番目とc番目の量子ビットについて、以下のような和(トレース)をとるイメージです。
\begin{align}
Tr_{a,b,c}(\rho) &= \sum_{i_a,i_b,i_c} \sum_{j_a,j_b,jc} \delta_{i_a,j_a} \delta_{i_b,j_b} \delta_{i_c,j_c} \bra{i_1,i_2,\cdots,i_n} \rho \ket{j_1,j_2,\cdots,j_n} \\
&= \sum_{i_a,i_b,i_c} \bra{i_1,\cdots,i_a,\cdots,i_b,\cdots,i_c,\cdots,i_n} \rho \ket{j_1,\cdots,i_a,\cdots,i_b,\cdots,i_c,\cdots,j_n}
\end{align}
というわけで、テキトーな疑似コードで示します。元の密度演算子(行列)はA[k][l]、部分トレース後の密度演算子(行列)はB[k][l]、元の行列の次元はNで、a番目,b番目,c番目の3つの量子ビットを部分トレースする想定です。部分トレースをとる量子ビットが3つ以外の場合については適当に読み替えてください。
B[k][l]のすべてのk,lについて0に初期化 (k,lは0から2^(n-3)-1までの整数)
for i=0 to N-1 do
k = [2進表示のiからa,b,c番目(MSBから数えて)のビットを削除して右詰め(LSB方向に)]
kk = [2進表示のiからa,b,c番目のビットを取り出して並べる]
for j=0 to N-1 do
l = [jに対して上のkと同じ処理]
ll = [jに対して上のkkと同じ処理]
if kk == ll then
B[k][l] = B[k][l] + A[i][j]
こんな感じで、qlazyに実装してみました。
シミュレータで確認
例として、4量子ビットの全体系を考えて、そのうちの2量子ビットの部分系がどうなるかを計算してみました。上で説明したように、全体系が積状態の場合、部分系は「純粋状態」になり、量子エンタングル状態の場合、部分系は「混合状態」になるはずなので、それを密度演算子の２乗のトレースをとることで確認します。
全体のPythonコードを示します。
【2021.9.5追記】qlazy最新版でのソースコードはここに置いてあります。
from qlazypy import QState, DensOp
qs = QState(4)
qs.h(0).h(1) # unitary operation for 0,1-system
qs.x(2).z(3) # unitary operation for 2,3-system
de1 = DensOp(qstate=[qs], prob=[1.0]) # product state
de1_reduced = de1.patrace(id=[0,1]) # trace-out 0,1-system
print("== partial trace of product state ==")
print(" * trace = ", de1_reduced.trace())
print(" * square trace = ", de1_reduced.sqtrace())
qs.cx(1,3).cx(0,2) # entangle between 0,1-system and 2,3-system
de2 = DensOp(qstate=[qs], prob=[1.0]) # entangled state
de2_reduced = de2.patrace(id=[0,1]) # trace-out 0,1-system
print("== partial trace of entangled state ==")
print(" * trace = ", de2_reduced.trace())
print(" * square trace = ", de2_reduced.sqtrace())
print("== partial state of entangled state ==")
qs.show(id=[2,3])
qs.free()
de1.free()
de2.free()
de1_reduced.free()
de2_reduced.free()
何をやっているか、順に説明します。
from qlazypy import QState, DensOp
量子状態クラスQStateと密度演算子クラスDensOpをimportします。
qs = QState(4)
qs.h(0).h(1) # unitary operation for 0,1-system
qs.x(2).z(3) # unitary operation for 2,3-system
4量子ビット状態を用意して、0番目と1番目の量子ビット(0,1系)に各々アダマールをかけます。2番目と3番目の量子ビット(2,3系)に各々XゲートとZゲートをかけます。これで、(0,1系)(2,3系)別個にユニタリ変換をした形になるので、純粋状態の積状態ができたことになります(qs)。量子回路で書くと、
0 --H--
1 --H--
2 --X--
3 --Z--
です。
de1 = DensOp(qstate=[qs], prob=[1.0]) # product state
de1_reduced = de1.patrace(id=[0,1]) # trace-out 0,1-system
この積状態(qs)に対する密度演算子を作成します(de1)。patraceメソッドで、(0,1系)についての部分トレースをとり、部分系(2,3系)に関する密度演算子を作成します(de1_reduced)。
print("== partial trace of product state ==")
print(" * trace = ", de1_reduced.trace())
print(" * square trace = ", de1_reduced.sqtrace())
部分系(2,3系)の密度演算子のトレースと２乗のトレースを計算して、表示します。純粋状態なので両方とも1になるはずです。
qs.cx(1,3).cx(0,2) # entangle between 0,1-system and 2,3-system
次に、現在の量子状態(qs)の部分系(0,1系)と(2,3系)をまたがる形でCNOTをかけます。先程(0,1系)にアダマールをかけていたので、これで量子エンタングル状態が発生します(するはずです)。量子回路で書くと、
0 --H------*---
1 --H---*--|---
2 --X---|--X---
3 --Z---X------
です。
de2 = DensOp(qstate=[qs], prob=[1.0]) # entangled state
de2_reduced = de2.patrace(id=[0,1]) # trace-out 0,1-system
print("== partial trace of entangled state ==")
print(" * trace = ", de2_reduced.trace())
print(" * square trace = ", de2_reduced.sqtrace())
先程と同じように、密度演算子を作成し(de2)、さらに部分トレースをとって(de2_deduced)、そのトレースと２乗のトレースを表示します。上で説明した通りであれば、混合状態になるので、トレースは1ですが、２乗のトレースは1以下の数になるはずです。
print("== partial state of entangled state ==")
qs.show(id=[2,3])
最後、量子状態の部分表示機能で、(2,3系)のみの量子状態を表示します。
結果は、以下の通りです。
== partial trace of product state ==
* trace = 1.0
* square trace = 1.0
== partial trace of entangled state ==
* trace = 1.0
* square trace = 0.25
== partial state of entangled state ==
c[00] = +0.0000+0.0000*i : 0.0000 |
c[01] = +1.0000+0.0000*i : 1.0000 |+++++++++++
c[10] = +0.0000+0.0000*i : 0.0000 |
c[11] = +0.0000+0.0000*i : 0.0000 |
量子エンタングルメント状態の部分系の密度演算子の２乗トレースが、確かに1以下(0.25)になりました。
最後の(2,3系)のみの状態表示ですが、この出力例では$\ket{01}$にピークが立っています。が、これは実行のたびに変わったりします。複数の純粋状態が確率的に出現する混合状態であるということを、まさに表しています。
おわりに
量子情報理論について、まだまだ基本の基の初歩段階です。測定の一般化とか、シュミット分解とか、いろいろ言うべき話題が残っていますので、シミュレータでの実装をできる範囲で実施しながら、理論の理解を進めていきたいです。次回は「測定の一般化」かなー、、と思っています。
以上