Quantum Cryptography Simulation in Python — Deep Dive

BB84 Implementation with Qiskit

Full Protocol Simulation

from qiskit import QuantumCircuit
from qiskit_aer import AerSimulator
import numpy as np

class BB84Protocol:
    def __init__(self, n_bits, eavesdrop=False):
        self.n_bits = n_bits
        self.eavesdrop = eavesdrop
        self.simulator = AerSimulator()

    def run(self):
        # Alice's random choices
        alice_bits = np.random.randint(0, 2, self.n_bits)
        alice_bases = np.random.randint(0, 2, self.n_bits)  # 0=Z, 1=X

        # Bob's random measurement bases
        bob_bases = np.random.randint(0, 2, self.n_bits)

        # Eve's random intercept bases (if eavesdropping)
        eve_bases = np.random.randint(0, 2, self.n_bits)

        bob_results = []

        for i in range(self.n_bits):
            qc = QuantumCircuit(1, 1)

            # Alice prepares qubit
            if alice_bits[i] == 1:
                qc.x(0)  # Set to |1⟩
            if alice_bases[i] == 1:
                qc.h(0)  # Rotate to X basis

            # Eve intercepts (if enabled)
            if self.eavesdrop:
                # Eve measures in her random basis
                if eve_bases[i] == 1:
                    qc.h(0)
                qc.measure(0, 0)
                # Eve must resend — simulated by re-preparing
                # based on her measurement result
                qc = self._eve_resend(qc, eve_bases[i])

            # Bob measures in his chosen basis
            if bob_bases[i] == 1:
                qc.h(0)
            qc.measure(0, 0)

            result = self.simulator.run(qc, shots=1).result()
            measured = int(list(result.get_counts().keys())[0])
            bob_results.append(measured)

        bob_results = np.array(bob_results)

        # Basis reconciliation: keep matching bases
        matching = alice_bases == bob_bases
        sifted_alice = alice_bits[matching]
        sifted_bob = bob_results[matching]

        return {
            'sifted_key_alice': sifted_alice,
            'sifted_key_bob': sifted_bob,
            'sifted_length': len(sifted_alice),
            'error_rate': self._compute_qber(sifted_alice, sifted_bob),
            'matching_bases': np.sum(matching),
        }

    def _eve_resend(self, original_qc, eve_basis):
        """Simulate Eve resending after measurement."""
        # In simulation: Eve's measurement collapses the state,
        # then she re-prepares based on her result
        qc = QuantumCircuit(1, 1)
        # The result of Eve's measurement determines her resend state
        # This is handled implicitly by the quantum circuit
        if eve_basis == 1:
            qc.h(0)
        return qc

    def _compute_qber(self, alice_key, bob_key):
        """Quantum Bit Error Rate."""
        if len(alice_key) == 0:
            return 0.0
        errors = np.sum(alice_key != bob_key)
        return errors / len(alice_key)

# Run without eavesdropper
bb84 = BB84Protocol(n_bits=1000, eavesdrop=False)
result = bb84.run()
print(f"No Eve - Sifted key length: {result['sifted_length']}")
print(f"No Eve - QBER: {result['error_rate']:.4f}")

# Run with eavesdropper
bb84_eve = BB84Protocol(n_bits=1000, eavesdrop=True)
result_eve = bb84_eve.run()
print(f"With Eve - Sifted key length: {result_eve['sifted_length']}")
print(f"With Eve - QBER: {result_eve['error_rate']:.4f}")
# Expected QBER with Eve: ~25%

Mathematical BB84 Simulation (Fast)

For large-scale analysis, skip the quantum circuit overhead:

def bb84_fast(n_bits, eavesdrop=False, channel_noise=0.0):
    """
    Fast mathematical simulation of BB84.
    Returns QBER and key material.
    """
    alice_bits = np.random.randint(0, 2, n_bits)
    alice_bases = np.random.randint(0, 2, n_bits)
    bob_bases = np.random.randint(0, 2, n_bits)

    # Start with perfect transmission
    received_bits = alice_bits.copy()

    if eavesdrop:
        eve_bases = np.random.randint(0, 2, n_bits)
        # Eve measures: wrong basis → 50% chance of flipping
        wrong_basis = eve_bases != alice_bases
        flip_mask = wrong_basis & (np.random.random(n_bits) < 0.5)
        received_bits = received_bits ^ flip_mask

    # Channel noise
    if channel_noise > 0:
        noise_mask = np.random.random(n_bits) < channel_noise
        received_bits = received_bits ^ noise_mask

    # Bob's measurement: wrong basis → random result
    wrong_basis_bob = bob_bases != alice_bases
    if eavesdrop:
        # After Eve, Bob's wrong-basis measurement on Eve's
        # re-prepared state is more complex
        random_mask = wrong_basis_bob & (np.random.random(n_bits) < 0.5)
    else:
        random_mask = wrong_basis_bob & (np.random.random(n_bits) < 0.5)
    bob_bits = received_bits ^ random_mask

    # Sifting
    matching = alice_bases == bob_bases
    sifted_alice = alice_bits[matching]
    sifted_bob = bob_bits[matching]

    qber = np.mean(sifted_alice != sifted_bob)

    return sifted_alice, sifted_bob, qber

# Statistical analysis over many runs
qbers_no_eve = []
qbers_with_eve = []
for _ in range(100):
    _, _, q = bb84_fast(10000, eavesdrop=False)
    qbers_no_eve.append(q)
    _, _, q = bb84_fast(10000, eavesdrop=True)
    qbers_with_eve.append(q)

print(f"QBER without Eve: {np.mean(qbers_no_eve):.4f} ± {np.std(qbers_no_eve):.4f}")
print(f"QBER with Eve: {np.mean(qbers_with_eve):.4f} ± {np.std(qbers_with_eve):.4f}")

E91 Protocol Implementation

from qiskit import QuantumCircuit
from qiskit_aer import AerSimulator
import numpy as np

class E91Protocol:
    """
    E91 QKD using entangled pairs and Bell inequality.
    Alice measures in bases: 0°, 45°, 90°
    Bob measures in bases: 45°, 90°, 135°
    """
    def __init__(self, n_pairs):
        self.n_pairs = n_pairs
        self.simulator = AerSimulator()

        # Measurement angles (in units of π/4)
        self.alice_angles = [0, np.pi/4, np.pi/2]
        self.bob_angles = [np.pi/4, np.pi/2, 3*np.pi/4]

    def create_bell_pair(self):
        """Create a Bell state |Φ⁺⟩ = (|00⟩ + |11⟩)/√2."""
        qc = QuantumCircuit(2, 2)
        qc.h(0)
        qc.cx(0, 1)
        return qc

    def measure_in_basis(self, qc, qubit, angle):
        """Measure qubit along a direction specified by angle."""
        qc.ry(-2 * angle, qubit)
        return qc

    def run(self):
        alice_basis_choices = np.random.randint(0, 3, self.n_pairs)
        bob_basis_choices = np.random.randint(0, 3, self.n_pairs)

        alice_results = []
        bob_results = []

        for i in range(self.n_pairs):
            qc = self.create_bell_pair()

            # Apply measurement rotations
            a_angle = self.alice_angles[alice_basis_choices[i]]
            b_angle = self.bob_angles[bob_basis_choices[i]]

            self.measure_in_basis(qc, 0, a_angle)
            self.measure_in_basis(qc, 1, b_angle)

            qc.measure([0, 1], [0, 1])

            result = self.simulator.run(qc, shots=1).result()
            bits = list(result.get_counts().keys())[0]
            alice_results.append(int(bits[1]))
            bob_results.append(int(bits[0]))

        alice_results = np.array(alice_results)
        bob_results = np.array(bob_results)

        # Key generation: use pairs where both chose matching angles
        # Alice basis 1 (45°) matches Bob basis 0 (45°)
        key_mask = (alice_basis_choices == 1) & (bob_basis_choices == 0)
        key_alice = alice_results[key_mask]
        key_bob = bob_results[key_mask]

        # Bell test: use other basis combinations
        bell_value = self._compute_bell(
            alice_results, bob_results,
            alice_basis_choices, bob_basis_choices
        )

        return {
            'key_alice': key_alice,
            'key_bob': key_bob,
            'key_length': len(key_alice),
            'bell_value': bell_value,
            'qber': np.mean(key_alice != key_bob) if len(key_alice) > 0 else 0,
        }

    def _compute_bell(self, a_res, b_res, a_bases, b_bases):
        """Compute CHSH Bell inequality value."""
        def correlation(a_base, b_base):
            mask = (a_bases == a_base) & (b_bases == b_base)
            if np.sum(mask) == 0:
                return 0
            a = 2 * a_res[mask] - 1  # Convert 0,1 to -1,+1
            b = 2 * b_res[mask] - 1
            return np.mean(a * b)

        # S = E(a1,b1) - E(a1,b3) + E(a3,b1) + E(a3,b3)
        S = (correlation(0, 0) - correlation(0, 2)
             + correlation(2, 0) + correlation(2, 2))
        return S

# Run E91
e91 = E91Protocol(n_pairs=5000)
result = e91.run()
print(f"Key length: {result['key_length']}")
print(f"QBER: {result['qber']:.4f}")
print(f"Bell value: {result['bell_value']:.3f}")
print(f"Bell violation (|S| > 2): {abs(result['bell_value']) > 2}")
# Classical limit: |S| ≤ 2
# Quantum maximum: |S| = 2√2 ≈ 2.83

Eavesdropping Strategies

Intercept-Resend Attack

The simplest attack — Eve measures and resends:

def intercept_resend_analysis(n_bits=10000, n_trials=50):
    """Analyze the intercept-resend attack detection probability."""
    qbers = []
    for _ in range(n_trials):
        _, _, qber = bb84_fast(n_bits, eavesdrop=True)
        qbers.append(qber)

    mean_qber = np.mean(qbers)
    print(f"Expected QBER: 25% (theoretical)")
    print(f"Measured QBER: {mean_qber*100:.1f}%")

    # Detection probability for checking n_check bits
    for n_check in [10, 20, 50, 100]:
        # Probability of detecting Eve with n_check test bits
        p_detect = 1 - (1 - mean_qber) ** n_check
        print(f"  Detection prob with {n_check} check bits: {p_detect:.6f}")

Beam-Splitting Attack

Eve splits each photon and measures her copy:

def beam_splitting_attack(n_bits, split_ratio=0.5):
    """
    Eve splits the photon beam, keeping a fraction.
    She can only measure her portion; the rest reaches Bob undisturbed.
    """
    alice_bits = np.random.randint(0, 2, n_bits)
    alice_bases = np.random.randint(0, 2, n_bits)
    bob_bases = np.random.randint(0, 2, n_bits)

    # Eve captures some photons entirely (split_ratio)
    # and forwards the rest untouched
    eve_intercepts = np.random.random(n_bits) < split_ratio

    # For intercepted photons: Eve measures, Bob gets nothing (or noise)
    eve_info = np.zeros(n_bits, dtype=int)
    eve_correct = np.zeros(n_bits, dtype=bool)

    eve_bases = np.random.randint(0, 2, n_bits)
    eve_correct_basis = eve_bases == alice_bases
    eve_info[eve_intercepts & eve_correct_basis] = alice_bits[
        eve_intercepts & eve_correct_basis
    ]

    # Bob receives photons Eve didn't intercept
    bob_bits = np.full(n_bits, -1)
    forwarded = ~eve_intercepts
    correct_basis = bob_bases == alice_bases
    bob_bits[forwarded & correct_basis] = alice_bits[forwarded & correct_basis]
    # Wrong basis → random
    wrong_basis = forwarded & ~correct_basis
    bob_bits[wrong_basis] = np.random.randint(0, 2, np.sum(wrong_basis))

    # Bob needs to report which photons he received
    received = forwarded
    sifted = received & correct_basis

    print(f"Photons forwarded: {np.sum(forwarded)}/{n_bits}")
    print(f"Sifted key bits: {np.sum(sifted)}")
    print(f"Eve's information on sifted bits: "
          f"{np.sum(eve_intercepts & correct_basis)} intercepted in key")

    return np.sum(sifted), np.sum(eve_intercepts & correct_basis)

Privacy Amplification

After error correction, Alice and Bob perform privacy amplification to eliminate any information Eve might have:

import hashlib

def privacy_amplification(key, eve_info_bits, security_parameter=10):
    """
    Reduce key length to remove Eve's partial information.

    key: raw sifted key (numpy array of bits)
    eve_info_bits: estimated bits of information Eve has
    security_parameter: additional bits to sacrifice for security
    """
    raw_length = len(key)
    final_length = raw_length - eve_info_bits - security_parameter

    if final_length <= 0:
        raise ValueError("Not enough key material after privacy amplification")

    # Use universal hash function (Toeplitz matrix method)
    # Simplified: use SHA-256 and truncate
    key_bytes = np.packbits(key).tobytes()
    hash_input = key_bytes

    # Generate final_length bits using hash
    final_key_bits = []
    counter = 0
    while len(final_key_bits) < final_length:
        h = hashlib.sha256(hash_input + counter.to_bytes(4, 'big')).digest()
        bits = np.unpackbits(np.frombuffer(h, dtype=np.uint8))
        final_key_bits.extend(bits)
        counter += 1

    return np.array(final_key_bits[:final_length])

# Example
raw_key = np.random.randint(0, 2, 500)
final_key = privacy_amplification(raw_key, eve_info_bits=50)
print(f"Raw key: {len(raw_key)} bits → Final key: {len(final_key)} bits")

Information-Theoretic Security Proof

The security of BB84 can be quantified:

def secure_key_rate(qber, basis_efficiency=0.5):
    """
    Calculate the asymptotic secure key rate for BB84.

    Rate = basis_efficiency * [1 - 2*h(QBER)]

    where h is the binary entropy function.
    """
    if qber >= 0.5:
        return 0.0

    def binary_entropy(p):
        if p == 0 or p == 1:
            return 0
        return -p * np.log2(p) - (1-p) * np.log2(1-p)

    rate = basis_efficiency * (1 - 2 * binary_entropy(qber))
    return max(0, rate)

# Key rate vs QBER
qbers = np.linspace(0, 0.15, 100)
rates = [secure_key_rate(q) for q in qbers]

# Security threshold: rate drops to 0 at QBER ≈ 11%
threshold = qbers[np.argmin(np.array(rates) > 0) - 1]
print(f"Security threshold QBER: {threshold:.3f}")
# Above this QBER, no secure key can be extracted

Channel Noise vs. Eavesdropping

A key challenge: how do you distinguish noise from eavesdropping?

def noise_vs_eve_analysis():
    """Compare QBER from channel noise vs eavesdropping."""
    noise_levels = [0.01, 0.03, 0.05, 0.08, 0.10]

    print("Scenario                  | QBER    | Secure?")
    print("-" * 52)

    for noise in noise_levels:
        # Pure channel noise
        _, _, q = bb84_fast(50000, eavesdrop=False, channel_noise=noise)
        rate = secure_key_rate(q)
        status = "✓" if rate > 0 else "✗"
        print(f"Noise={noise:.2f}, no Eve       | {q:.4f} | {status} (rate={rate:.3f})")

    print()
    # Eve + channel noise
    for noise in [0.00, 0.02, 0.05]:
        _, _, q = bb84_fast(50000, eavesdrop=True, channel_noise=noise)
        rate = secure_key_rate(q)
        status = "✓" if rate > 0 else "✗"
        print(f"Noise={noise:.2f}, with Eve     | {q:.4f} | {status} (rate={rate:.3f})")

noise_vs_eve_analysis()

In practice, QKD systems must distinguish between physical noise and eavesdropping. Error correction handles noise; privacy amplification handles Eve. The total QBER from both combined must stay below the security threshold.

Real-World QKD Simulation

def realistic_qkd_simulation(
    distance_km=50,
    fiber_loss_db_per_km=0.2,
    detector_efficiency=0.1,
    dark_count_rate=1e-6,
    source_rate_hz=1e9,
    protocol='bb84'
):
    """
    Simulate QKD performance with realistic physical parameters.
    """
    # Channel transmission
    total_loss_db = fiber_loss_db_per_km * distance_km
    channel_transmission = 10 ** (-total_loss_db / 10)

    # Detection probability
    p_detect = channel_transmission * detector_efficiency

    # QBER from dark counts
    p_dark = dark_count_rate / source_rate_hz
    qber_dark = p_dark / (2 * (p_detect + p_dark))

    # Raw key rate
    raw_rate = source_rate_hz * p_detect * 0.5  # 0.5 for basis matching

    # Secure key rate
    total_qber = qber_dark  # Add other noise sources as needed
    skr_fraction = secure_key_rate(total_qber)
    secure_rate = raw_rate * skr_fraction

    return {
        'distance_km': distance_km,
        'channel_loss_db': total_loss_db,
        'detection_prob': p_detect,
        'qber': total_qber,
        'raw_key_rate_bps': raw_rate,
        'secure_key_rate_bps': secure_rate,
    }

# Performance vs distance
for d in [10, 50, 100, 200, 300]:
    result = realistic_qkd_simulation(distance_km=d)
    print(f"{d:>3} km: loss={result['channel_loss_db']:.0f}dB, "
          f"QBER={result['qber']:.6f}, "
          f"secure rate={result['secure_key_rate_bps']:.0f} bps")

One thing to remember: Simulating quantum cryptography in Python isn’t just an academic exercise — it’s how researchers design real QKD systems, model attack strategies, optimize protocols for specific hardware, and prove security bounds before building expensive physical infrastructure. The simulation code and the physics it models are two sides of the same coin.