Zenn
Open7

密度行列・部分トレース

derwindderwind

量子状態 \ket{\psi} の密度行列は \ket{\psi}\bra{\psi} で与えられる。\ket{\psi} = \ket{\psi_1} \otimes \ket{\psi_0} というテンソル積で書かれているとして、\ket{\psi_1} = \begin{bmatrix} a \\ b \end{bmatrix}, \ket{\psi_0} = \begin{bmatrix} c \\ d \end{bmatrix} とする。\ket{\psi}^\dagger を計算したい。

\begin{align*} \ket{\psi_1} \otimes \ket{\psi_0} &= \begin{bmatrix} ac \\ ad \\ bc \\ bd \end{bmatrix} \\ &\xrightarrow{\dagger} [a^*c^*\ \ a^*d^*\ \ b^*c^*\ \ b^*d^*] = [a^*\ \ b^*] \otimes [c^*\ \ d^*] = \bra{\psi_1} \otimes \bra{\psi_0} \end{align*}

となる。よって、\ket{\psi}^\dagger = \bra{\psi} である。

上記から、\ket{\psi}\bra{\psi} = (\ket{\psi_1} \otimes \ket{\psi_0}) (\bra{\psi_1} \otimes \bra{\psi_0}) = (\ket{\psi_1}\bra{\psi_1}) \otimes (\ket{\psi_0}\bra{\psi_0}) が従う。

derwindderwind

一般には、 \rho = \begin{bmatrix}\rho_{00} & \rho_{01} \\ \rho_{10} & \rho_{11}\end{bmatrix} が密度行列とは、\operatorname{tr}(\rho) = 1, \rho \geq 0, \rho: エルミートの時を言う。

この時、スペクトル分解定理から、スペクトルを \{ p_i \} として、直交基底 \{ \ket{\psi_i} \} により \rho = \sum_i p_i \ket{\psi_i}\bra{\psi_i} と書け、仮定より p_i \geq 0、\sum p_i = 1 が従う。

つまり、量子システムが多くの状態 \ket{\psi_i} の可能性があるとしてその確率を p_i とする時に純粋状態の線型結合で記述されている形になる。

どの状態 \ket{\psi} にあるかが正確に分かっている量子システムを純粋状態にあるという。つまり \rho = \ket{\psi}\bra{\psi} の時である。

例) \rho = \frac{1}{2} \begin{bmatrix}1 & 1 \\ 1 & 1\end{bmatrix} = \ket{+}\bra{+} = (\frac{\ket{0} + \ket{1}}{\sqrt{2}})(\frac{\bra{0} + \bra{1}}{\sqrt{2}}) は純粋状態である。

純粋状態でないものを混合状態と呼ぶ。

例) \rho = \frac{1}{2} \begin{bmatrix}1 & 0 \\ 0 & 1\end{bmatrix} = \frac{1}{2} \ket{0}\bra{0} + \frac{1}{2} \ket{1}\bra{1} は混合状態である。

derwindderwind

\ket{\psi} = \cos \frac{\theta}{2} \ket{0} + e^{i\phi} \cos \frac{\theta}{2} \ket{1} とする。この時、（この純粋状態に対する）密度行列 \rho は、

\begin{align*} \rho = \ket{\psi} \bra{\psi} &= \begin{bmatrix} \cos \frac{\theta}{2} \\ e^{i\phi} \sin \frac{\theta}{2} \end{bmatrix} \begin{bmatrix} \cos \frac{\theta}{2}\quad e^{-i\phi} \sin \frac{\theta}{2} \end{bmatrix} \\ &= \begin{bmatrix} \cos^2 \frac{\theta}{2} & e^{-i\phi} \sin \frac{\theta}{2} \cos \frac{\theta}{2} \\ e^{i\phi} \sin \frac{\theta}{2} \cos \frac{\theta}{2} & \sin^2 \frac{\theta}{2} \end{bmatrix} \\ &= \begin{bmatrix} \cos^2 \frac{\theta}{2} & \frac{1}{2} \sin \theta [\cos \phi - i \sin \phi] \\ \frac{1}{2} \sin \theta [\cos \phi + i \sin \phi] & \sin^2 \frac{\theta}{2} \end{bmatrix} \tag{1} \end{align*}

となる。

\begin{align*} X = \begin{bmatrix} 0 & 1 \\ 1 & 0 \end{bmatrix},\ Y = \begin{bmatrix} 0 & -i \\ i & 0 \end{bmatrix},\ Z = \begin{bmatrix} 1 & 0 \\ 0 & -1 \end{bmatrix} \end{align*}

に注意して、\frac{2 \cos^2 \alpha - 1}{2} = \frac{1}{2}(\cos^2 \alpha - \sin^2 \alpha) = \frac{1}{2} \cos (2\alpha) を使うと (1) は

\begin{align*} \rho &= \frac{1}{2} \! \begin{bmatrix} 1 & 0 \\ 0 & 1 \end{bmatrix} \!+\! \begin{bmatrix} \cos^2 \! \frac{\theta}{2} \!-\! \frac{1}{2} & 0 \\ 0& - \left(\cos^2 \! \frac{\theta}{2} \!-\! \frac{1}{2} \right) \end{bmatrix} \!+\! \begin{bmatrix} 0 & \frac{1}{2} \! \sin \theta [\cos \! \phi \!-\! i \sin \! \phi] \\ \frac{1}{2} \! \sin \theta [\cos \! \phi \!-\! i \sin \! \phi] & 0 \end{bmatrix} \\ &= \frac{1}{2} I + \frac{1}{2} \! \cos \! \theta Z + \begin{bmatrix} 0 & \frac{1}{2} \! \sin \! \theta [\cos \! \phi - i \sin \! \phi] \\ \frac{1}{2} \sin \theta [\cos \phi - i \sin \phi] & 0 \end{bmatrix} \\ &= \frac{1}{2} I + \frac{1}{2} \cos \theta Z + \frac{1}{2} \sin \theta \cos \phi X + \frac{1}{2} \sin \theta \sin \phi Y \\ &= \frac{1}{2} \left( I + xX + yY + zZ \right) \end{align*}

と書ける。ここで、\operatorname{tr}(X) = \operatorname{tr}(Y) = \operatorname{tr}(Z) = \operatorname{tr}(XY) = \operatorname{tr}(YZ) = \operatorname{tr}(ZX) = 0 に注意して、

\begin{align*} \begin{cases} x = \sin \theta \cos \phi = \operatorname{tr}(X\rho)\\ y = \sin \theta \sin \phi = \operatorname{tr}(Y\rho) \\ z = \cos \theta = \operatorname{tr}(Z\rho) \\ \end{cases} \end{align*}

である。よって ブロッホベクトル を思い出すと、密度行列を成分分解した係数と 1 対 1 対応していることが分かる。

また、正規直交基底 \{ \ket{i} \} を使うと、

\begin{align*} \braket{X} &= \braket{\psi | X | \psi} \\ &= \bra{\psi} \sum_i \braket{i | X | \psi} \ket{i} \\ &= \sum_i \braket{i | X | \psi} \braket{\psi | i} = \operatorname{tr}(X \ket{\psi} \bra{\psi}) = \operatorname{tr}(X \rho) \end{align*}

であるので、x = \braket{X} である。同様に、y = \braket{Y}, z = \braket{Z} が分かる。

この辺の話は、量子コンピュータと量子通信I 演習 2.72 や 量子コンピュータと量子通信III 8.4.2 量子プロセストモグラフィーにも似たような解説が記載されている。

derwindderwind
qc = QuantumCircuit(1)
qc.x(0)

という量子回路がある場合に、X ゲートにノイズが乗っていて 確率 p \in (0, 1) で状態 \ket{0} のままとする。この場合、終状態での密度行列は確率 1-p で \ket{1}\bra{1} で確率 p で \ket{0}\bra{0} である。つまり、

\begin{align*} \rho &= (1-p) \ket{1}\bra{1} + p \ket{0}\bra{0} \\ &= \begin{bmatrix} p & 0 \\ 0 & 1 - p \end{bmatrix} \end{align*}

という Bloch 球の内部に対応する混合状態になる（p = \frac{1}{2} の時については上で触れた）。x = \operatorname{tr}(X \rho) = 0, y = \operatorname{tr}(Y \rho) = 0, z = \operatorname{tr}(Z \rho) = 2p - 1 に注意すると、\braket{Z} = 2p - 1 となる。よって、\braket{X}^2 + \braket{Y}^2 + \braket{Z}^2 = (2p - 1)^2 < 1 となる。
純粋状態 \ket{\psi} = \alpha \ket{0} + \beta \ket{1},\ |\alpha|^2 + |\beta|^2 = 1 の場合には、

\begin{align*} &\ \braket{X}^2 + \braket{Y}^2 + \braket{Z}^2 \\ =&\ \braket{\psi | X | \psi}^2 + \braket{\psi | Y | \psi}^2 + \braket{\psi | Z | \psi}^2 \\ =&\ (\alpha^* \beta + \beta^* \alpha)^2 + (-i \alpha^* \beta + i \beta^* \alpha)^2 + (\alpha^* \alpha - \beta^* \beta)^2 \\ =&\ 2 \alpha \alpha^* \beta \beta^* + |\alpha|^4 + |\beta|^4 \\ = & (\alpha^* \alpha + \beta^* \beta)^2 = (|\alpha|^2 + |\beta|^2)^2 = 1 \end{align*}

となることに注意したい。

derwindderwind

X ゲートにノイズが乗っていて 確率 p \in (0, 1) で状態 \ket{0} のままとする。

\rho = (1-p) \ket{1}\bra{1} + p \ket{0}\bra{0} について、確率 p で余計な flip が発生すると、回路のユニタリは実質 U_B = XX になっている。確率 1-p で正常に動作するほうは回路のユニタリは U_A = X である。よって、\rho = (1-p) U_A\ket{0}\bra{0}U_A^\dagger + p U_B\ket{0}\bra{0}U_B^\dagger と書き直すことができる。

より一般的に書くと、

\begin{align*} \rho = \sum_K p_K U_K \ket{0} \bra{0} U_K^\dagger,\quad K \in \{A, B\} \end{align*}

となる。これを踏まえてエラーチャネルを

\begin{align*} \mathcal{E}(\rho_\text{in}) = \sum_K p_K U_K\rho_\text{in} U_K^\dagger \end{align*}

と書く。純粋状態を始状態とする密度行列を \rho_\text{in} とする場合、既に見たようにこれは \rho_\text{in} = \frac{1}{2}(I + x_\text{in} X + y_\text{in} Y + z_\text{in} Z) と書ける。X \rho_\text{in} X = \frac{1}{2}(I + x_\text{in} X - y_\text{in} Y - z_\text{in} Z) を直接計算で得るので、今回のノイズ入り X ゲートの回路のケースでは、

\begin{align*} \rho_\text{out} &= \mathcal{E}(\rho_\text{in}) \\ &= (1-p) \frac{1}{2}(I + x_\text{in} X + y_\text{in} Y + z_\text{in} Z) + p \frac{1}{2}(I + x_\text{in} X - y_\text{in} Y - z_\text{in} Z) \\ &= \frac{1}{2} (I + x_\text{in} X + (1-2p) y_\text{in} Y + (1-2p) z_\text{in} Z) \end{align*}

となる。上で見た密度行列と ブロッホベクトル との対応付けを思い出すと、

\begin{align*} \begin{bmatrix} x_\text{in} \\ (1-2p) y_\text{in} \\ (1-2p) z_\text{in} \end{bmatrix} \end{align*}

に対応することになるので、純粋状態から始めた結果、Bloch 球が Y 軸と Z 軸方向に潰れたような楕円球にうつる形 (“shrink our Bloch sphere”) になる。

derwindderwind

X \frac{1}{2} (I + xX + yY + zZ) X = \frac{1}{2} (I + xX - yY - zZ) であったので、上で見た \rho_\text{out} = \mathcal{E}(\rho_\text{in}) の式を合わせると、ノイズ入りの X ゲートを d 回適用、つまり

qc = QuantumCircuit(1)
for _ in range(d):
    qc.x(0)

すると、終状態の密度行列は \rho_\text{final} = \frac{1}{2} (I + x_\text{in} X + (2p-1)^d y_\text{in} Y + (2p-1)^d z_\text{in} Z) となることが分かる。
特に、\ket{\psi} = \ket{0} を始状態とする時、\rho_\text{in} = \frac{1}{2}(I + Z) であり、終状態の密度行列は \frac{1}{2}(I + (2p-1)^d Z) となる。この時、\braket{Z} = \operatorname{tr}(Z \rho_\text{final}) = (2p-1)^d である。p \ll 1 とすると、Taylor 展開の近似で、(2p-1)^d = (-1)^d (1-2p)^d \approx (-1)^d (1-2dp) となる。理想的な X ゲートでは終状態は \frac{1}{2}(I + (-1)^d Z) であり、\braket{Z}_\text{ideal} = (-1)^d であるので、\braket{Z}_\text{ideal} - \braket{Z} \approx 2 (-1)^d dp というように、p について線形なエラーが乗ることが分かる。（インコヒーレントエラー）

derwindderwind

部分トレース

合成系 \ket{\psi} = \ket{\psi_A}\ket{\psi_B} を考える。さらに、系の密度行列 \rho = \ket{\psi}\!\bra{\psi} = \ket{\psi_A}\ket{\psi_B}\!\bra{\psi_A}\bra{\psi_B} を考える。

この時、\rho から系 A の密度行列は以下のように求まる:

\begin{align*} \rho_A = \operatorname{tr}_B (\rho) &= \operatorname{tr}_B (\ket{\psi_A}\ket{\psi_B}\!\bra{\psi_A}\bra{\psi_B}) \\ &= \operatorname{tr}_B (\ket{\psi_A}\!\bra{\psi_A} \otimes \ket{\psi_B}\!\bra{\psi_B}) \\ &= \ket{\psi_A}\!\bra{\psi_A} \operatorname{tr}(\ket{\psi_B}\!\bra{\psi_B}) = \ket{\psi_A}\!\bra{\psi_A} \end{align*}

ここで、系 B の正規直交基底を \{ \ket{\phi_j} \} とすると、

\begin{align*} \operatorname{tr}_B (\ket{\psi_A}\!\bra{\psi_A} \otimes \ket{\psi_B}\!\bra{\psi_B}) &= \sum_j I_A \otimes \bra{\phi_j} (\ket{\psi_A}\!\bra{\psi_A} \otimes \ket{\psi_B}\!\bra{\psi_B}) I_A \otimes \ket{\phi_j} \\ &= \ket{\psi_A}\!\bra{\psi_A} \otimes \sum_j \braket{\phi_j | \psi_B}\!\braket{\psi_B | \phi_j} \\ &= \ket{\psi_A}\!\bra{\psi_A} \sum_j |\braket{\phi_j | \psi_B}|^2 \\ &= \ket{\psi_A}\!\bra{\psi_A} |\ket{\psi_B}|^2 = \ket{\psi_A}\!\bra{\psi_A} \end{align*}

となる。要するに、

\begin{align*} \rho_A = \ket{\psi_A}\!\bra{\psi_A} = \sum_j (I_A \otimes \bra{\phi_j}) \rho (I_A \otimes \ket{\phi_j}) \end{align*}

である。純粋状態による合成系のみ記述したが、混合状態でも同様の計算である。

Python

import numpy as np


I = np.eye(2)
ZERO = np.array([1, 0], dtype=complex).reshape(-1, 1)
ONE = np.array([0, 1], dtype=complex).reshape(-1, 1)


def kron(*args):
    if len(args) < 2:
        raise ValueError("At least two arrays are required for Kronecker product")
    
    return reduce(np.kron, args)


_ = kron
state = _((ZERO+ONE)/np.sqrt(2), (ZERO-ONE)/np.sqrt(2))
rho = state @ state.conj().T

subsystem_rho1 = 0
for psi in [ZERO, ONE]:
    subsystem_rho1 += _(I, psi).conj().T @ rho @ _(I, psi)
print(subsystem_rho1)
evals, evecs = np.linalg.eigh(subsystem_rho1)
subsystem_state1 = evecs[:, np.argmax(evals)]
print(subsystem_state1)

subsystem_rho2 = 0
for psi in [ZERO, ONE]:
    subsystem_rho2 += _(psi, I).conj().T @ rho @ _(psi, I)
print(subsystem_rho2)
evals, evecs = np.linalg.eigh(subsystem_rho2)
subsystem_state2 = evecs[:, np.argmax(evals)]
print(subsystem_state2)