Quantum-Computation-course-Basel/practice/PracticeExercise3.ipynb

Quantum teleportation and dense coding with GHZ and W states

Each task requires you to do use the code proposed below, or, more directly, use Qiskit

1. Can standard teleportation be efficiently done with GHZ state?

In standard quantum teleportation, Alice and Bob share a Bell pair, and Alice can teleport an unknown qubit state to Bob using two classical bits.

In this exercise, you will explore how teleportation can (or cannot) be achieved using a GHZ state shared among three parties — Alice, Bob, and Charlie.

A 3-qubit GHZ state is defined as:

|\text{GHZ}\rangle = \frac{1}{\sqrt{2}} (|000\rangle + |111\rangle).

Tasks:

1. Prepare a GHZ state.
2. Build the state with Qubit layout: |\psi\rangle |\text{ghz}_a\,\, \text{ghz}_c\, \,\text{ghz}_b\rangle. |\psi\rangle=\alpha|0\rangle + \beta|1\rangle (qubit 0) is the state Alice wants to send. We place GHZ on qubits 1,2,3 where Alice holds qubit 1, Charlie qubit 2, Bob qubit 3. We give Alice access to |\psi\rangle and her GHZ qubit (|\psi\rangle and |\text{ghz}_a\rangle), while Bob holds |\text{ghz}_b\rangle. Charlie holds |\text{ghz}_c\rangle but we will consider the case where Charlie does not cooperate.
3. Let Alice perform a Bell-basis measurement on (|\psi\rangle, |\text{ghz}_a\rangle) — same as in standard teleportation.
4. Condition on Alice's measurement outcome, allow Bob to apply the usual Pauli corrections based only on Alice's classical bits, but without any action from Charlie.
5. Examine Bob's reduced state and compute fidelity with the original |\psi\rangle.

In [ ]:

import numpy as np
from numpy import sqrt, exp, pi
from qiskit.quantum_info import Statevector, partial_trace, DensityMatrix
from itertools import product

In [ ]:

# single qubit state |psi> parameterized by angles (theta, phi)
def single_qubit_state(theta=0.7, phi=0.3):
    return np.array([np.cos(theta/2), np.exp(1j*phi)*np.sin(theta/2)], dtype=complex)

# prepare initial global state: |psi> (q0) ⊗ |GHZ(1,2,3)>
def prepare_global_state(psi_vec):
    # GHZ on qubits 1,2,3: (|000> + |111>)/sqrt(2)
    ghz = (1/sqrt(2)) * np.array([1,0,0,0,0,0,0,1], dtype=complex)  # 3 qubits
    #global_state = kronecker product the two parts to have 1 ⊗ 3 = 4 qubits total
    return global_state

In [ ]:

# Pauli matrices and Bell basis projectors
X = np.array([[0,1],[1,0]], dtype=complex)
# define Y, Z, I similarly

# Bell states
phi_plus = (1/sqrt(2)) * np.array([1,0,0,1], dtype=complex)
#define phi_minus, psi_plus and psi_minus similarly
bells = [phi_plus, phi_minus, psi_plus, psi_minus]

In [ ]:

#define an arbitrary single qubit state and the resultant 4 qubit resource
# psi_vec =   single_qubit_state(?)
# global_state = prepare_global_state(?)

In [ ]:

# projector onto Bell state for qubits (0,1) in a 4-qubit system: build as kron(P_bell, I_2, I_2)
def bell_projector_on_01(bvec):
    P = np.outer(bvec.conj(), bvec)
    # embed into 4-qubit space: qubits 0 and 1 are the Bell qubits
    return np.kron(P, np.kron(I, I))

In [ ]:

# fidelity between pure psi and density matrix rho on a single qubit
def fidelity_single_qubit(psi_vec, rho):
    # psi_vec is vector, rho is 2x2 density matrix
    return np.real(np.vdot(psi_vec, rho @ psi_vec))

In [ ]:

results = []
for idx, bvec in enumerate(bells):
    P = bell_projector_on_01(bvec)
    # Project global state onto Bell outcome (unnormalized)
    proj_state = P @ global_state
    prob = np.vdot(proj_state, proj_state).real
    if prob < 1e-12:
        continue
    # Normalize post-measurement global state
    proj_state = proj_state / sqrt(prob)
    # After Alice's measurement outcome idx, usually Bob would apply a correction depending on idx:
    # Mapping: phi_plus -> I, phi_minus -> Z, psi_plus -> X, psi_minus -> XZ (up to phase)
    if idx == 0:  # hint: phi_plus -> I (no correction)
        corr = np.kron(np.eye(1), np.kron(np.eye(1), np.eye(2)))  # placeholder (unused)
        # For 4-qubit vector, we will apply correction on qubit 3 (Bob's qubit)
        corr_on_bob = ?
    elif idx == 1: 
        corr_on_bob = ?
    elif idx == 2: 
        corr_on_bob = ?
    else:  
        corr_on_bob = ?

    # Apply correction on qubit 3: form operator I(0) ⊗ I(1) ⊗ I(2) ⊗ corr_on_bob
    Ucorr = np.kron(np.eye(8), corr_on_bob)  # 4x4 ⊗ 2x2 = 8x8
    state_after_corr = Ucorr @ proj_state

    # Now compute Bob's reduced state by tracing out qubits 0,1,2 -> keep qubit 3
    # reshape to 4-qubit density and partial trace
    rho = np.outer(state_after_corr, state_after_corr.conj())
    # trace out qubits 0,1,2 by reshaping
    rho_bob = np.zeros((2,2), dtype=complex)
    # indices: total 4 qubits -> size 16; iterate basis states and sum over traced subsystems
    for i in range(16):
        for j in range(16):
            # bitstrings for i and j
            bi = format(i, '04b')
            bj = format(j, '04b')
            # keep only elements where first three bits equal (we are tracing them)
            if bi[0:3] == bj[0:3]:
                # accumulate contributions to bob density matrix element (last bit)
                ib = int(bi[3])
                jb = int(bj[3])
                rho_bob[ib, jb] += rho[i, j]
    F = fidelity_single_qubit(psi_vec, rho_bob)
    results.append((idx, prob, F))

In [ ]:

# Print results
print("Bell outcome index, probability, fidelity of Bob's reduced state with original |psi>")
for r in results:
    print(r)
# compute average fidelity over outcomes (weighted)
avgF = sum(prob * F for (_, prob, F) in results) / sum(prob for (_, prob, F) in results)
print("\nAverage fidelity (with no Charlie cooperation) = {:.6f}".format(avgF))

6. why is the fidelity <1 without Charlie's intervention?
7. If instead of Charlie, the third qubit of the GHZ state actually modelizes the entanglement of the qubits a and b with the environment (and therefore, this qubit is not accessible to any operator), conclude on the qualitative effect of this entanglement (i.e. a kind of noise) on the teleportation protocol.

2. Can superdense coding (particularly the case of two classical bits sent by sharing one qubit) be efficiently done with GHZ and W states?

Examining the original superdense coding architecture, the basic principle is that Alice, who wants to transmit two classical bits of information, performs local unitaries on her qubit and four orthogoonal states (Bell states) are produced. We'll explicitly test whether Alice's local unitaries on her single qubit can produce 4 mutually orthogonal global states when three qubit GHZ and W states are used as resources.

The W state is defined as: |W\rangle = \frac{1}{\sqrt{3}} (|100\rangle + |010\rangle + |001\rangle)

Tasks:

1. Build the GHZ and W states on 3 qubits (labelled 0, 1, 2).
2. On Alice's qubit (qubit 0) apply the four usual encoding unitaries: I, X, Z, Y=(iXZ) on both the GHZ and W states.
3. Compute the eight resulting global statevectors and form the Gram matrix of overlaps to test orthogonality (for W state, repeat and observe that the states are not mutually orthogonal).
4. Conclude.

In [ ]:

# Build GHZ state with ordering |q0 q1 q2> where index = q0*4 + q1*2 + q2
ghz = np.zeros(8, dtype=complex)
ghz[0] = 1.0
ghz[7] = 1.0
ghz = ghz / np.linalg.norm(ghz)   # (|000> + |111>)/sqrt(2)

In [ ]:

# Build W state similar to above (|100> + |010> + |001>)/sqrt(3)
#w = np.zeros(...)
#...
#...
#...
w = w / np.linalg.norm(w)

In [ ]:

# Alice has qubit 0; Bob have qubits 1 and 2
# Define encoding unitaries on Alice's qubit: I, X, Z, XZ
encs = [I, X, Z, X @ Z]
enc_names = ['I','X','Z','XZ']

In [ ]:

def apply_on_alice(U, state):
    # apply U on qubit 0 of 3-qubit state: U ⊗ I ⊗ I
    return np.kron(U, np.kron(I, I)) @ state

def gram_matrix(states):
    n = len(states)
    G = np.zeros((n,n), dtype=complex)
    for i in range(n):
        for j in range(n):
            G[i,j] = np.vdot(states[i], states[j])
    return G

In [ ]:

# GHZ case
#ghz_states =?... apply single qubit unitaries U on qubit 0 and identities on qubit 1 and 2
#G_ghz = ?... calculate the Gram matrix
print("GHZ Gram matrix (absolute values):\n", np.round(np.abs(G_ghz), 6))

# Check if GHZ states are mutually orthogonal
orth_ghz = np.allclose(G_ghz, np.diag(np.diag(G_ghz)), atol=1e-8)
print("GHZ: are encoded states orthogonal?", orth_ghz)

In [ ]:

# W case
# go through the same steps as above, are the encoded states orthogonal?

3. Can superdense coding of three classical bits be done by sharing two qubits with the GHZ state?

If Alice now has access to two qubits of the GHZ state (the 0th and 1st qubit) and Bob holds the last qubit (2nd), then Alice can send three classical bits by encoding them in an orthogonal set of 8 states.​ Conceptually, this means Alice can performs any two qubit unitaries from \{I,X,Y,Z\}_0\otimes \{I,X,Y,Z\}_1 (or any 1 and 2 qubits gates, but we restrict ourselves to Pauli in this exercise) on her qubits and the result should contain 8 orthogonal states.

1. Apply the all the Pauli unitaries above to the GHZ state.

In [ ]:

# kronecker for two-qubit operator acting on qubits (0,1) with qubit order q0 q1 q2.
def apply_on_alice(U_ab, state):
    # U_ab is 4x4 acting on qubits 0 and 1; full operator is U_ab ⊗ I (Bob qubit is last)
    #return [?]

In [ ]:

# Build two-qubit operator list: all tensor products of single-qubit Paulis
single_paulis = [('I',I), ('X',X), ('Y',Y), ('Z',Z)]
op_list = []
op_names = []
for (n1, M1), (n2, M2) in product(single_paulis, repeat=2):
    name = f"{n1}⊗{n2}"
    op = np.kron(M1, M2)
    op_list.append(op)
    op_names.append(name)

In [ ]:

# Compute encoded states
encoded_states = [apply_on_alice(U, ghz) for U in op_list]
# Ensure normalization (should be unitary applied to normalized state)
for v in encoded_states:
    assert np.allclose(np.linalg.norm(v), 1.0, atol=1e-9)

# Compute Gram matrix
G = np.array([[np.vdot(sj, si) for sj in encoded_states] for si in encoded_states], dtype=complex)
absG = np.round(np.abs(G), 8)
print("Absolute Gram matrix (|<psi_i|psi_j>|):")
print(absG)

In [ ]:

# Rank (how many orthogonal states)
# Use SVD tolerance
u,s,vt = np.linalg.svd(G)
rank = np.sum(s > 1e-8)
print("\n Numeric rank of Gram matrix (max # orthogonal states):", rank)

Understandably, not all 16 combinations are needed, and 8 are enough.

2. Find such a set.

In [ ]:

# pick an orthogonal subset
orth_indices = []
tol = 1e-8
for i in range(len(encoded_states)):
    v = encoded_states[i]
    ok = True
    for j in orth_indices:
        ov = abs(np.vdot(encoded_states[j], v))
        if ov > tol:
            ok = False
            break
    if ok:
        orth_indices.append(i)
    if len(orth_indices) >= rank:
        break

print("\n Found orthogonal subset of size", len(orth_indices))
for k in orth_indices:
    print(" -", op_names[k])

3. Conclude