Skip to content
Quantum Mechanics

Simulating Quantum Channels

A quantum-channel calculation should answer two questions: what happens to the state, and is the proposed map physically allowed even when the system is entangled with a reference? This example computes both. It uses small complex matrices so that you can compare every numerical result with an analytic formula.

Helpful background. Kraus operators describe channel action; the Choi matrix tests complete positivity. The code needs Python 3 and NumPy; paste the complete block into a notebook cell or a .py file.

Channels as linear maps on matrices

Section titled “Channels as linear maps on matrices”

Use the basis (∣0⟩,∣1⟩)(|0\rangle,|1\rangle), with ∣0⟩|0\rangle the ground state for amplitude damping. A Kraus map acts on any matrix AA by

Φ(A)=∑rKrAKr†.\Phi(A)=\sum_r K_r A K_r^\dagger.

It preserves trace when ∑rKr†Kr=I\sum_rK_r^\dagger K_r=I. For damping probability pp, the operators are

K0=(1001−p),K1=(0p00).K_0=\begin{pmatrix}1&0\\0&\sqrt{1-p}\end{pmatrix}, \qquad K_1=\begin{pmatrix}0&\sqrt p\\0&0\end{pmatrix}.

The distinction between a general matrix and a density matrix matters when building a Choi matrix. For example, the linear depolarizing map is

Φλ(A)=λA+(1−λ)Tr⁡(A)I2.\Phi_\lambda(A)=\lambda A+(1-\lambda)\operatorname{Tr}(A)\frac I2.

The trace factor cannot be omitted: off-diagonal basis matrices have trace zero. Substituting I/2I/2 indiscriminately would produce the wrong Choi matrix.

We use output first, reference second:

JΦ=∑i,j=01Φ(Eij)⊗Eij,Eij=∣i⟩⟨j∣.J_\Phi=\sum_{i,j=0}^1\Phi(E_{ij})\otimes E_{ij}, \qquad E_{ij}=|i\rangle\langle j|.

For a channel, JΦ≥0J_\Phi\ge0 and Tr⁡outJΦ=I\operatorname{Tr}_{\mathrm{out}}J_\Phi=I. Its trace is 22, not 11; JΦ/2J_\Phi/2 is the associated normalized bipartite state.

Compute an output state and its Choi matrix

Section titled “Compute an output state and its Choi matrix”

The input ∣+y⟩=(∣0⟩+i∣1⟩)/2|+y\rangle=(|0\rangle+i|1\rangle)/\sqrt2 tests complex conjugation as well as populations. At p=0.36p=0.36 the expected output is

Φ(ρy)=(0.68−0.40i0.40i0.32).\Phi(\rho_y)=\begin{pmatrix}0.68&-0.40i\\0.40i&0.32\end{pmatrix}.

Here @ means matrix multiplication, .conj().T means the adjoint, and np.kron forms a tensor product.

import numpy as np


I = np.eye(2, dtype=complex)
tol = 1e-12
p = 0.36
Ks = [np.diag([1, np.sqrt(1-p)]),
      np.array([[0, np.sqrt(p)], [0, 0]])]


def apply(Ks, A):
    return sum((K @ A @ K.conj().T for K in Ks),
               np.zeros_like(A, dtype=complex))


def choi(phi):
    J = np.zeros((4, 4), dtype=complex)
    for i in range(2):
        for j in range(2):
            E = np.zeros((2, 2), dtype=complex)
            E[i, j] = 1
            J += np.kron(phi(E), E)
    return J


def trace_output(J):
    # Indices are output-row, input-row, output-column, input-column.
    return np.einsum('aiaj->ij', J.reshape(2, 2, 2, 2))


rho_y = np.array([[0.5, -0.5j], [0.5j, 0.5]])
out = apply(Ks, rho_y)
J = choi(lambda A: apply(Ks, A))
assert np.linalg.norm(sum(K.conj().T @ K for K in Ks)-I) < tol
assert np.linalg.norm(J-J.conj().T) < tol
assert np.linalg.eigvalsh(J).min() >= -tol
assert np.linalg.norm(trace_output(J)-I) < tol
assert np.allclose(out, [[0.68, -0.4j], [0.4j, 0.32]])
assert np.allclose(np.linalg.eigvalsh(J), [0, 0, p, 2-p])
print('output state:', out, sep='\n')
print('Choi eigenvalues:', np.linalg.eigvalsh(J))


lam = -0.4
bad_J = choi(lambda A: lam*A + (1-lam)*np.trace(A)*I/2)
assert np.allclose(np.linalg.eigvalsh(bad_J), [-0.1, 0.7, 0.7, 0.7])
print('invalid depolarizing Choi:', np.linalg.eigvalsh(bad_J))


# Apply a Hadamard rotation first, then damping.
H = np.array([[1, 1], [1, -1]]) / np.sqrt(2)
composed = apply([K @ H for K in Ks], rho_y)
sequential = apply(Ks, H @ rho_y @ H.conj().T)
assert np.allclose(composed, sequential)
assert not np.allclose(composed, H @ out @ H.conj().T)

The valid Choi eigenvalues are 0,0,0.36,1.640,0,0.36,1.64 up to roundoff. The deliberately invalid depolarizing map has eigenvalues −0.1,0.7,0.7,0.7-0.1,0.7,0.7,0.7. It nevertheless sends every isolated qubit state to a valid state, because it shrinks its Bloch vector by −0.4-0.4. Only the reference-system test reveals the failure of complete positivity.

Check Hermiticity before calling eigvalsh: that routine assumes a Hermitian matrix and reads only one triangle. A mistaken non-Hermitian matrix must not be silently interpreted as Hermitian. For these small exact formulas, 10−1210^{-12} is a reasonable absolute tolerance; a genuine eigenvalue of −0.1-0.1 is not roundoff.

Try other channels

Section titled “Try other channels”

Replace Ks by [U][U] for a unitary channel, or by

[1+λ2I,1−λ2Z]\left[\sqrt{\frac{1+\lambda}{2}}I, \sqrt{\frac{1-\lambda}{2}}Z\right]

for dephasing with −1≤λ≤1-1\le\lambda\le1. Compare the outputs for ∣0⟩|0\rangle, ∣1⟩|1\rangle, ∣+⟩|+\rangle, ∣+y⟩|+y\rangle, and I/2I/2. The last state distinguishes unital noise from amplitude damping. The four matrices I/2I/2 and (I+σj)/2(I+\sigma_j)/2 span the qubit operator space, so agreement on this set determines a linear map, apart from numerical error.

For composition, matrix order follows operation order: if Ψ\Psi acts first and Φ\Phi second, the Kraus operators are KrLsK_rL_s. The final code lines check this and demonstrate that reversing the physical order generally changes the answer. Continue with Bloch-vector noise models to see the same maps geometrically.

Exercises

Section titled “Exercises”

Completeness. Verify the amplitude-damping completeness relation by hand.

Solution

K0†K0=diag⁡(1,1−p)K_0^\dagger K_0=\operatorname{diag}(1,1-p) and K1†K1=diag⁡(0,p)K_1^\dagger K_1=\operatorname{diag}(0,p). Their sum is II.

Unitality. Calculate the image of I/2I/2 at p=0.36p=0.36.

Solution

The result is diag⁡(0.68,0.32)\operatorname{diag}(0.68,0.32), so amplitude damping is not unital. Identity, unitary, dephasing, and depolarizing channels are unital.

Complete-positivity boundary. Derive the allowed depolarizing interval from its Choi eigenvalues.

Solution

The eigenvalue along the maximally entangled vector is (1+3λ)/2(1+3\lambda)/2; the three orthogonal eigenvalues are (1−λ)/2(1-\lambda)/2. Requiring both to be nonnegative gives −1/3≤λ≤1-1/3\le\lambda\le1.

A wrong tensor-factor trace. What does tracing the reference factor of the damping Choi matrix give?

Solution

It gives Φ(I)=diag⁡(1+p,1−p)\Phi(I)=\operatorname{diag}(1+p,1-p), which is not II for p>0p>0. This is a unitality test, not the trace-preservation test. The nonunital example makes a swapped partial trace detectable.

References

Section titled “References”