量子情報理論の基本：密度演算子

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

## はじめに

密度演算子自体は量子力学の基本的な概念かもしませんが、量子情報の理論展開にはとても重要な役割を果たしています。これを出発点にして、量子情報理論の基本について、これからちょっと勉強していきたいと思います。理論の理解ができたところで、自作の量子計算シミュレータ[qlazy](https://github.com/samn33/qlazy)を使って、その概念・性質を確認してみます。

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

* [密度演算子-純粋状態・混合状態-量子コンピュータ入門-物理とか](https://whyitsso.net/physics/quantum_mechanics/QI6.html)
* [ニールセン、チャン「量子コンピュータと量子通信（２）」オーム社（2005年）](https://www.ohmsha.co.jp/book/9784274200083/)

## 純粋状態と混合状態

密度演算子の説明に入る前に、「純粋状態」と「混合状態」という概念をおさえておきます。

量子力学を学ぶと、まず状態の重ね合わせという不思議な考え方が登場し、面食らいます。例えば、1粒子の2準位系(スピンのようなもの)の場合、

```math
\ket{\psi} = \alpha \ket{0} + \beta \ket{1}
```

という具合に記述され、これは状態$\ket{0}$(例えばスピン上向き)と状態$\ket{1}$(例えばスピン下向き)の重ね合わせ(どっちでもあり得るという状態)を表現するものであって、上向きか下向きかは測定によってはじめて決定されるものであり、その出現確率は$\alpha$や$\beta$の絶対値の２乗で計算できます、というわけです。一般的には、以下のように量子状態を記述することができます。

```math
\ket{\psi} = c_0 \ket{u_0} + c_1 \ket{u_1} + \cdots
```

考えている系に応じて、シュレディンガー方程式をつくって解けば、この量子状態は完全に求まります。大事なのは、この重ね合わせの表式自体は確率的ではなく、決定論的に求まるということです(物理量は確率的にしかわからないのですが、このあたり混乱しないように)。このような量子状態のことを「純粋状態」と呼んでいます。重ね合わさり方が1種類しかなく、他のものが混じっていないという意味で「純粋」と理解すれば良いと思います。たとえて言うなら、オレ独自の香辛料の割合でつくったオレのカレールーのようなものです。その割合はオレにしかわからない、測定しないと他人にはわからない(企業秘密なので)、けれど、割合自体は決まっている、みたいな感じですかね(で、伝わる？)。

一方で、重ね合わさり方が複数あって(つまり$\\{\ket{\psi_{i}}\\}$のような集合)、その出現割合が確率的にしかわからない場合(つまり$\\{p_{i},\ket{\psi_{i}}\\}$というセットの集合)も考えることができます。ここで言う確率は量子的な確率ではなく各々の純粋状態がどのくらいの頻度で登場するかを表す統計的な確率と理解してください。このような状態＝量子状態のアンサンブルのことを「混合状態」と呼んでいます。先程のカレーのたとえでいうと、オレのカレールー以外にもAさんのカレールーやBさんのカレールーもあって、どのカレールーがお店(カレー屋さん)で出てくるか、確率的にしかわからないみたいな状況でしょうか(で、伝わる？)。

妙なたとえ話はここまでにして、「純粋状態」と「混合状態」の違いについて、もう少し具体的な物理イメージで示してみます。2準位のスピン系を例にします。例えば、状態$\ket{0}$と状態$\ket{1}$が半々に重なり合った状態は、

```math
\ket{\psi} = \frac{1}{\sqrt{2}} (\ket{0} + \ket{1})
```

と書けます。何らかの装置からスピンをもった粒子が飛び出してきて、それが上向きなのか下向きなのかわかっていない状況です(測定すると半々の確率で観測される)。これは「純粋状態」です。一方で、純粋状態で記述される複数の量子状態が確率的に飛び出してくるような状況もあり得ます。例えば、純粋状態$\ket{0}$と純粋状態$\ket{1}$が半々の割合で飛び出してくるような場合です。これは「混合状態」です。「純粋状態」と「混合状態」のこの例の場合、どっちも測定すると上向きと下向きが半々になるので、結局、「純粋」も「混合」も現象としては同じということで、わざわざ「純粋」とか「混合」とか分けなくても良いのでは、という声も聞こえてきそうです。が、違います。測定する方向を変えると、両者は決定的に違う結果を与えます。例えば{$\ket{+},\ket{-}$}方向に測定することを考えます。つまり、ブロッホ球でいうとX軸方向の測定です。とすると、前者は状態$\ket{+}$そのものなので、100%の確率でX軸上向きという結果になります。一方、後者の場合、

```math
\begin{align}
\ket{0} &= \frac{1}{\sqrt{2}} (\ket{+} + \ket{-}) \\
\ket{1} &= \frac{1}{\sqrt{2}} (\ket{+} - \ket{-})
\end{align}
```

と書き直してみればわかる通り、X軸上向きと下向きが半々の確率で出現するということになります。というわけで両者はまったく違うものと理解しておくべきです。ちなみに、カレールーのたとえ話は、この時点で破綻します。量子論的なカレールーなるものをでっち上げないと話が続きません。なので、たとえ話は話半分にしておいてください(逆に、なぜこのたとえがダメなのかを考えることで、量子力学の理解が少し進むと思います)。


## 密度演算子

### 定義

この「純粋状態」と「混合状態」を統一的に記述するための概念が「密度演算子」です。我々が相手にするような系は常に閉じていて、その中でずっと純粋状態を保ち続けるという理想的な状態にあるわけではなく、むしろ外界と相互作用をしながら純粋状態が次第に壊れていくという状況(結果として混合状態になるという状況)が普通のありがちな状態であります。そのようなリアルな現実を記述したいとするならば、「密度演算子」というのは非常に強力なツールになります。

では、その定義を示します。いま複数の(純粋)量子状態$\\{\ket{\psi_i}\\}\space(i=1,2,\cdots)$があって、各々の出現確率が$\\{p_i\\} \space (i=1,2,\cdots)$であるとすると、この量子状態のアンサンブルに対応した「密度演算子」は、

```math
\rho = \sum_{i} p_{i} \ket{\psi_{i}} \bra{\psi_{i}}  \tag{1}
```

と定義されます。純粋状態１個に対する「密度演算子」は、これの特殊な場合で、

```math
\rho = \ket{\psi} \bra{\psi}
```

と書けます。

### 時間発展

純粋状態に対する時間発展はユニタリ演算子$U$を用いたユニタリ変換で表現できます。つまり、

```math
\ket{\psi_i} \rightarrow U \ket{\psi_i}
```

でした。密度演算子の時間発展も、

```math
\begin{align}
\rho &= \sum_{i} p_i \ket{\psi_i} \bra{\psi_i} \\
&\rightarrow \sum_{i} p_i U \ket{\psi} \bra{\psi_i} U^{\dagger} = U \rho U^{\dagger} \tag{2}
\end{align}
```

ということで、$\rho$という演算子に対するユニタリ変換で表現できます。

### 測定

測定という過程はどのように密度演算子で表現できるでしょうか。いま簡単のため射影測定をイメージしておきます(別に射影測定でなくても良いのですが)。で、その測定に対応した射影演算子の集まりを$\\{M_m\\} \space (m=1,2,\cdots)$とします。初期状態$\ket{\psi_i}$に対して測定した結果、$m$が得られる確率は、

```math
p(m|i) = \bra{\psi_i} M_m^{\dagger} M_m \ket{\psi_i} = Tr(M_m^{\dagger} M_m \ket{\psi_i} \bra{\psi_i})
```

となりますので、量子状態のアンサンブルがあったときに結果$m$が得られる確率$p(m)$は、

```math
\begin{align}
p(m) &= \sum_{i} p(m|i) p_i \\
&= \sum_{i} p_i Tr(M_m^{\dagger} M_m \ket{\psi_i} \bra{\psi_i}) \\
&= Tr(M_m^{\dagger} M_m \rho)  \tag{3}
\end{align}
```

と書けます。

また、初期状態$\ket{\psi_i}$に対して測定結果$m$を得たあとの状態は、

```math
\ket{\psi_i^m} = \frac{M_m \ket{\psi_i}}{\sqrt{\bra{\psi_i} M_m^{\dagger} M_m \ket{\psi_i}}}
```

です。これに対する密度演算子$\rho_m$は、

```math
\begin{align}
\rho_m &= \sum_{i} p(i|m) \ket{\psi_i^m} \bra{\psi_i^m} \\
&= \sum_{i} p(i|m) \frac{M_m \ket{\psi_i} \bra{\psi_i} M_m^{\dagger}}{\bra{\psi_i} M_m^{\dagger} M_m \ket{\psi_i}} \\
&= \sum_{i} p_i \frac{M_m \ket{\psi_i} \bra{\psi_i} M_m^{\dagger}}{Tr(M_m^{\dagger} M_m \rho)} \\
&= \frac{M_m \rho M_m^{\dagger}}{Tr(M_m^{\dagger} M_m \rho)}  \tag{4}
\end{align}
```

と書けます。ここで、条件付き確率に対するベイズの定理、

```math
p(i|m) = \frac{p(m|i)p_i}{p(m)}
```

を使いました。

ということで、測定のプロセスを密度演算子と測定演算子を使って記述できることがわかりました。


### 数学的な性質

密度演算子の定義から、

* (1) 半正定値である
* (2) トレースが1である

ということが言えます。また、この２つの性質が成り立つ演算子であれば、それは密度演算子であるということも逆に言えます。

まず、(1)について。「半正定値」というのは、任意の$\ket{\phi}$に対して、$\bra{\phi} \rho \ket{\phi} \geq 0$が成り立つということです。密度演算子の定義から、

```math
\bra{\phi} \rho \ket{\phi} = \sum_{i} p_i \braket{\phi}{\psi_i} \braket{\psi_i}{\phi} = \sum_{i} p_i |\braket{\phi}{\psi_i}|^2 \geq 0  \tag{5}
```

ということで、半正定値であることがわかります。

次に、(2)について。密度演算子の定義から、

```math
Tr(\rho) = \sum_{i} p_i Tr(\ket{\psi_i} \bra{\psi_i}) = \sum_{i,j} p_i \braket{j}{\psi_i} \braket{\psi_i}{j} = \sum_{i,j} p_i \braket{\psi_i}{j} \braket{j}{\psi_i} = \sum_{i} p_i = 1 \tag{6}
```

ということで、確かに1になります。ここで$\ket{j}$は、いま考えている空間上での適当な正規直交基底です。上の式からわかる通り、これは確率の和が1であるという事実を反映した性質になります。

逆に「半正定値である」かつ「トレースが1である」演算子であれば、それは密度演算子である、ということも以下のように証明できます。

半正定値であれば、密度演算子$\rho$の固有値$\\{\lambda_i\\}$、固有ベクトル$\\{\ket{u_i}\\}$を使って、以下のようにスペクトル分解できます。

```math
\rho = \sum_{i} \lambda_i \ket{u_i} \bra{u_i}
```

このとき、半正定値なので、$\lambda_i$は実数で$\lambda_i \geq 0$です。また、$Tr(\rho)=1$より、$\sum_i \lambda_i = 1$です。これは密度演算子の定義にほかなりません。というわけで、証明完了です。

### 複合系

量子状態が$\ket{\Psi_{A}}, \space \ket{\Psi_{B}}, \space \ket{\Psi_{C}}, \space \cdots$という具合に複数あったとき、この全体を表す量子状態は、それらのテンソル積、

```math
\ket{\Psi} = \ket{\Psi_{A}} \otimes \ket{\Psi_{B}} \otimes \ket{\Psi_{C}}, \space, \cdots
```

で表すことができました。同様に、この複合系を密度演算子で、

```math
\rho = \rho_{A} \otimes \rho_{B} \otimes \rho_{C} \otimes \cdots
```

という具合にテンソル積で表すことができます。

ここで、複合系と混合状態を混同しないようにしてください。先程のカレーの例で言うと(しつこい？)、ラーメンの世界や餃子の世界が新たに登場して、１つの世界を形成するイメージが複合系です(謎)。冗談はさておき、真面目な物理イメージを言うと、注目している量子系Aとそれを取り囲む外界を表す量子系Bからなる全体世界を、複合系の典型例の１つと思っておけば良いと思います。

### ユニタリな自由度

密度演算子を使っていく上で１つ注意事項があります。それは、異なる量子状態アンサンブルが ~~異なる~~ 同じ密度演算子を与えることがあるということです。具体例を示します。

* 例1

```math
\ket{\phi_1} =
\begin{pmatrix}
1 \\
0
\end{pmatrix}, \space
\ket{\phi_2} =
\begin{pmatrix}
0 \\
1
\end{pmatrix}
```

という2つの量子状態の出現確率が半々$(p_1 = 1/2,\space p_2 = 1/2)$の場合

* 例2

```math
\ket{\psi_1} = \frac{1}{\sqrt{2}}
\begin{pmatrix}
1 \\
1
\end{pmatrix}, \space
\ket{\psi_2} = \frac{1}{\sqrt{2}}
\begin{pmatrix}
1 \\
-1
\end{pmatrix}
```

という2つの量子状態の出現確率が半々$(q_1 = 1/2, \space q_2 = 1/2)$の場合

例1、例2ともに、密度演算子を計算すると、

```math
\rho = \sum_{i} p_i \ket{\phi_i} \bra{\phi_i} = \sum_{i} q_i \ket{\psi_i} \bra{\psi_i} = \frac{1}{2}
\begin{pmatrix}
1 & 0 \\
0 & 1
\end{pmatrix}
```

となり、同じになります。

一般に、密度演算子の側から見ると、量子状態アンサンブルにはユニタリの自由度があると言えます。2つの量子状態アンサンブル、$\\{p_i, \ket{\phi_i}\\}, \space \\{q_i, \ket{\psi_i}\\}$があったとします。記号を簡単化するため、$\ket{\tilde{\phi_i}} = \sqrt{p_i } \ket{\phi_i}, \space \ket{\tilde{\psi_i}} = \sqrt{q_i} \ket{\psi_i}$とおきます。この2つの間に、あるユニタリ演算子$U(行列要素はu_{ij})$があって、

```math
\ket{\tilde{\phi_i}} = \sum_{j} u_{ij} \ket{\tilde{\psi_j}}

```

が成り立っているとすると、

```math
\begin{align}
\rho_{\phi} &= \sum_{i} \ket{\tilde{\phi_i}} \bra{\tilde{\phi_i}} \\
&= \sum_{i} \sum_{k,l} u_{ik} \ket{\tilde{\psi_k}} u_{li}^{*} \bra{\tilde{\psi_l}} \\
&= \sum_{k,l} \delta_{kl} \ket{\tilde{\psi_k}} \bra{\tilde{\psi_l}} \\
&= \sum_{k} \ket{\tilde{\psi_k}} \bra{\tilde{\psi_k}} = \rho_{\psi}
\end{align}
```

となり、両者の密度演算子は等しくなります。つまり、量子状態アンサンブルに対する、上のようなユニタリ変換に対して、密度演算子の値は保存されます(変化しません)。


## 純粋状態と混合状態の判別

純粋状態と混合状態を統一的に記述することができる重要な概念として密度演算子を導入しましたが、実は、これを使って、系が純粋状態なのか混合状態なのかを判別することができるということを示します。結論を先に言います。

```math
Tr(\rho^2) = 1  \tag{7}
```

の場合、純粋状態であり、

```math
Tr(\rho^2) \leq 1  \tag{8}
```

の場合、混合状態であると判別できます。証明は以下です。

【証明】

$\rho$が正定値なので、固有値$\lambda_i$と固有ベクトル(固有状態)$\ket{i}$を使って、以下のように展開できます。

```math
\rho = \sum_{i} \lambda_i \ket{i} \bra{i}
```

ここで、$\lambda_i \leq 0$で、$\rho$のトレースが1であることから、$\sum_{i} \lambda_i = 1$が成り立ちます。このとき$\rho^2$を計算すると、

```math
\begin{align}
\rho^2 &= \sum_{i,j} \lambda_i \lambda_j \ket{i} \braket{i}{j} \bra{j} \\
&= \sum_{i} \lambda_i^2 \ket{i} \bra{i}
\end{align}
```

となりますので、そのトレースは、

```math
\begin{align}
Tr(\rho^2) &= Tr(\sum_{i} \lambda_i^2 \ket{i} \bra{i}) \\
&= \sum_{i,k} \lambda_i^2 \braket{k}{i} \braket{i}{k} \\
&= \sum_{i} \lambda_i^2 \leq (\sum_{i} \lambda_i)^2 = 1
\end{align}
```

となります。

純粋状態の場合、$\lambda_i$は1つしかないので、必然的に$\sum_i \lambda_i^2 = 1$となります。つまり、純粋状態の場合のみ、$Tr(\rho^2) = 1$であり、その他の場合、$Tr(\rho^2) \leq 1$ということになります。(証明終)


## シミュレータで確認

量子計算シミュレータ[qlazy](https://github.com/samn33/qlazy)に密度演算子の機能を追加したので(v0.0.22)、その動作を確認してみます。純粋状態と混合状態の判別が密度演算子の２乗のトレースをとることでできることを、たった今示しましたので、それについて確認してみます。具体的には、$\frac{1}{\sqrt{2}} (\ket{0} + \ket{1})$という純粋状態と、$\ket{0}$と$\ket{1}$が半々で混じり合っている混合状態で試してみます。

全体のPythonコードは以下の通りです。

【2021.9.5追記】[qlazy](https://github.com/samn33/qlazy)最新版でのソースコードは[ここ](https://github.com/samn33/qlazy/blob/master/example/py/DensityOperator)に置いてあります。

```python
from qlazypy import QState,DensOp

qs_pure = QState(1).h(0)     # (|0> + |1>) / sqrt(2.0)
de_pure = DensOp(qstate=[qs_pure], prob=[1.0])

qs_pure_1 = QState(1)       # |0>
qs_pure_2 = QState(1).x(0)  # |1>
de_mixed = DensOp(qstate=[qs_pure_1,qs_pure_2], prob=[0.5,0.5])

print("== pure state ==")
de_pure.show()
print("* trace =", de_pure.trace())
print("* square trace =", de_pure.sqtrace())
print("")
print("== mixed state ==")
de_mixed.show()
print("* trace =", de_mixed.trace())
print("* square trace =", de_mixed.sqtrace())
    
qs_pure.free()
qs_pure_1.free()
qs_pure_2.free()

de_pure.free()
de_mixed.free()
```

何をやっているか、順に説明します。

```python
from qlazypy import QState,DensOp
```

通常の量子状態を表すQStateクラスと密度演算子を表すDensOpをimportします。

```python
qs_pure = QState(1).h(0)     # (|0> + |1>) / sqrt(2.0)
de_pure = DensOp(qstate=[qs_pure], prob=[1.0])
```

QStateクラスを使って、1量子ビットにアダマールをかけることで、純粋状態$\frac{1}{\sqrt{2}} (\ket{0} + \ket{1})$に対応したインスタンスを作成し、変数qs_pureに格納しています。さらに、DensOpクラスへの引数として、いま作成した量子状態の配列と確率の配列(純粋状態なので要素は1つしかないです(値は[1.0]))を与え、密度演算子のインスタンスを作成して、変数de_pureに格納しています。

```python
qs_pure_1 = QState(1)       # |0>
qs_pure_2 = QState(1).x(0)  # |1>
de_mixed = DensOp(qstate=[qs_pure_1,qs_pure_2], prob=[0.5,0.5])
```

純粋状態$\ket{0}$と$\ket{1}$を作成し、各々変数qs_pure_1,qs_pure_2に格納しています。さらに、DensOpクラスへの引数として、その配列および想定している確率[0.5,0.5]を与えて、混合状態を作成し、変数de_mixedに格納しています。

```python
print("== pure state ==")
de_pure.show()
print("* trace =", de_pure.trace())
print("* square trace =", de_pure.sqtrace())
print("")
print("== mixed state ==")
de_mixed.show()
print("* trace =", de_mixed.trace())
print("* square trace =", de_mixed.sqtrace())
```

純粋状態de_pureおよびde_mixedに対して、showメソッドで密度演算子(行列)の要素を表示します。traceメソッドでそのトレースを計算します。sqtraceで２乗のトレースを計算します。

結果は、以下の通りです。

```
== pure state ==
elm[0][0] = +0.5000+0.0000*i : 0.2500 |+++
elm[0][1] = +0.5000+0.0000*i : 0.2500 |+++
elm[1][0] = +0.5000+0.0000*i : 0.2500 |+++
elm[1][1] = +0.5000+0.0000*i : 0.2500 |+++
* trace = 1.0
* square trace = 1.0

== mixed state ==
elm[0][0] = +0.5000+0.0000*i : 0.2500 |++++
elm[0][1] = +0.0000+0.0000*i : 0.0000 |
elm[1][0] = +0.0000+0.0000*i : 0.0000 |
elm[1][1] = +0.5000+0.0000*i : 0.2500 |++++
* trace = 1.0
* square trace = 0.5
```

ということで、両方ともトレースは1で、かつ、２乗のトレースは純粋状態に対しては1、混合状態に対しては1以下(0.5)となることが確認できました(この程度の計算であればシミュレータを使わずとも簡単に手計算できるのですが、とりあえずデバッグも兼ねて)。


## おわりに

密度演算子の基本についてひとまず説明してみましたが、１つ大事なことを省いてしまいました。それは、密度演算子だけを使って、量子力学の公理を完全に言い換えることができるということです。そもそも量子力学の公理って何？ということから始めるのが、ちょっと面倒になりました。すみません。ご興味あれば、はじめに示した参考文献とか、その他適当な書籍なりをご参照くださいませ。

今後、密度演算子を起点にして、複合系で部分トレースとか、測定の一般化とか、CPTP写像とか、シュミット分解とか、いろいろ面白そうな話題が量子情報理論にはありますので、このあたりの理解も進めていきたいと思います

以上