Quantum-Computation-course-Basel/extra_resources/Entanglement_extra/12_Purification.ipynb

In [1]:

"""
### Overview
Simulate a simple entanglement purification step (BBPSSW-like) on isotropic states and track fidelity gains.

Exercises:
- Chain multiple rounds and observe convergence threshold behavior.
"""

Out [1]:

'\n### Overview\nSimulate a simple entanglement purification step (BBPSSW-like) on isotropic states and track fidelity gains.\n\nExercises:\n- Chain multiple rounds and observe convergence threshold behavior.\n'

In [1]:

import numpy as np
def bbpssw_round(F):
    # Toy mapping for isotropic states under BBPSSW (not exact; illustrative)
    # Successful output fidelity (heuristic):
    F_out = (F**2 + ((1-F)/3)**2) / (F**2 + 2*F*(1-F)/3 + 5*((1-F)/3)**2)
    p_succ = F**2 + 2*F*(1-F)/3 + 5*((1-F)/3)**2
    return F_out, p_succ
for F in [0.5,0.6,0.7,0.8,0.9]:
    Fo, ps = bbpssw_round(F)
    print(f"F_in={F:.2f}  F_out≈{Fo:.3f}  p_succ≈{ps:.3f}")

F_in=0.50  F_out≈0.500  p_succ≈0.556
F_in=0.60  F_out≈0.620  p_succ≈0.609
F_in=0.70  F_out≈0.735  p_succ≈0.680
F_in=0.80  F_out≈0.838  p_succ≈0.769
F_in=0.90  F_out≈0.926  p_succ≈0.876

In [4]:

import numpy as np
from qiskit import QuantumCircuit
from qiskit.quantum_info import Statevector, DensityMatrix, partial_trace

qc_phip = QuantumCircuit(2)
qc_phip.h(0)
qc_phip.cx(0, 1)
rho_phip = DensityMatrix.from_instruction(qc_phip)

qc_phim = QuantumCircuit(2)
qc_phim.x(0)
qc_phim.h(0)
qc_phim.cx(0, 1)
rho_phim = DensityMatrix.from_instruction(qc_phim)


qc_psip = QuantumCircuit(2)
qc_psip.x(1)
qc_psip.h(0)
qc_psip.cx(0, 1)
rho_psip = DensityMatrix.from_instruction(qc_psip)


qc_psim = QuantumCircuit(2)
qc_psim.x(0)
qc_psim.x(1)
qc_psim.h(0)
qc_psim.cx(0, 1)
rho_psim = DensityMatrix.from_instruction(qc_psim)

print("Bell States Density Matrices:")
print("Phi+:\n", rho_phip)
print("Phi-:\n", rho_phim)
print("Psi+:\n", rho_psip)
print("Psi-:\n", rho_psim)

Bell States Density Matrices:
Phi+:
 DensityMatrix([[0.5+0.j, 0. +0.j, 0. +0.j, 0.5+0.j],
               [0. +0.j, 0. +0.j, 0. +0.j, 0. +0.j],
               [0. +0.j, 0. +0.j, 0. +0.j, 0. +0.j],
               [0.5+0.j, 0. +0.j, 0. +0.j, 0.5+0.j]],
              dims=(2, 2))
Phi-:
 DensityMatrix([[ 0.5+0.j,  0. +0.j,  0. +0.j, -0.5+0.j],
               [ 0. +0.j,  0. +0.j,  0. +0.j,  0. +0.j],
               [ 0. +0.j,  0. +0.j,  0. +0.j,  0. +0.j],
               [-0.5+0.j,  0. +0.j,  0. +0.j,  0.5+0.j]],
              dims=(2, 2))
Psi+:
 DensityMatrix([[0. +0.j, 0. +0.j, 0. +0.j, 0. +0.j],
               [0. +0.j, 0.5+0.j, 0.5+0.j, 0. +0.j],
               [0. +0.j, 0.5+0.j, 0.5+0.j, 0. +0.j],
               [0. +0.j, 0. +0.j, 0. +0.j, 0. +0.j]],
              dims=(2, 2))
Psi-:
 DensityMatrix([[ 0. +0.j,  0. +0.j,  0. +0.j,  0. +0.j],
               [ 0. +0.j,  0.5+0.j, -0.5+0.j,  0. +0.j],
               [ 0. +0.j, -0.5+0.j,  0.5+0.j,  0. +0.j],
               [ 0. +0.j,  0. +0.j,  0. +0.j,  0. +0.j]],
              dims=(2, 2))

In [5]:

rho_id   = DensityMatrix(np.eye(4) / 4)

def isotropic_state(F):
	"""Create an isotropic state with fidelity F to |Φ+>."""
	return F * rho_phip + (1 - F) / 3 * (rho_phim + rho_psip + rho_psim)

def bbpssw_circuit():
	"""Construct the BBPSSW purification circuit."""
	circ = QuantumCircuit(4, 2)
	# CNOTs
	circ.cx(0, 2)
	circ.cx(1, 3)
	return circ

def simulate_bbpssw(F):
	"""Simulate one round of BBPSSW purification on isotropic states."""
	# Prepare initial state
	rho1 = isotropic_state(F)
	rho2 = isotropic_state(F)
	rho_in = rho1.tensor(rho2)
	proj = rho_id.tensor(rho_phip+rho_phim)

	# Create BBPSSW circuit
	circ = bbpssw_circuit()

	# Apply circuit to the state
	rho_out = DensityMatrix(circ).evolve(rho_in)

	# Trace out measured qubits (2 and 3) over only successful outcomes and renormalize
	rho_reduced = partial_trace(DensityMatrix(proj.data @ rho_out.data), [2, 3])
	norm = np.trace(rho_reduced.data)
	rho_reduced = DensityMatrix(rho_reduced.data / norm)

	# Calculate fidelity with |Φ+>
	fid = np.real(np.trace(rho_reduced.data @ rho_phip.data))

	return fid, rho_reduced

for F in [0.5, 0.6, 0.7, 0.8, 0.9]:
	Fo, rho_final = simulate_bbpssw(F)
	print(f"F_in={F:.2f}  F_out={Fo:.3f}")

F_in=0.50  F_out=0.900
F_in=0.60  F_out=0.953
F_in=0.70  F_out=0.980
F_in=0.80  F_out=0.993
F_in=0.90  F_out=0.999

In [ ]: