What is Quantum Entanglement?

Quantum entanglement is one of the most discussed and most misunderstood concepts in all of physics. This tutorial gives you a precise, practical understanding of what entanglement actually is, how to recognize it mathematically, and why it matters for quantum computing. Along the way, you will build working code that creates, measures, and quantifies entanglement across a range of scenarios.

What Entanglement Is Not

Before defining entanglement, it helps to address five specific misconceptions head-on.

Misconception 1: “Entanglement allows faster-than-light signaling.” When Alice measures her half of an entangled pair, she instantly learns something about what Bob would get if he measured his half. But Alice cannot choose her measurement outcome. Her results look completely random to her, and Bob’s results look completely random to him. Only when they later compare notes (using a classical channel, limited by the speed of light) do they discover that their results are correlated. Because Alice cannot control her outcome, she cannot encode a message for Bob. The no-signaling theorem makes this precise: the reduced density matrix on Bob’s side is completely independent of what Alice chooses to measure.

Misconception 2: “Entanglement is spooky action at a distance.” Einstein used this phrase in frustration, arguing that quantum mechanics must be incomplete. Decades of Bell test experiments have confirmed that quantum correlations are real and that no local hidden variable theory can reproduce them. Yet quantum mechanics is “local” in the operational sense: it satisfies no-signaling. The correlations look non-local when you compare measurement records after the fact, but they cannot transmit information. The correct framing is that entanglement reveals correlations that have no classical explanation, not that it involves any physical action propagating faster than light.

Misconception 3: “Measuring one qubit tells you the other’s exact state.” This is only true for maximally entangled states measured in the right basis. If Alice and Bob share the Bell state (|00> + |11>)/sqrt(2) and Alice measures in the computational basis, she learns Bob’s state with certainty: if she gets |0>, Bob has |0>. But for a partially entangled state like cos(theta)|00> + sin(theta)|11>, Alice’s measurement gives her a conditional probability distribution over Bob’s outcomes, not a definite value (unless theta = pi/4). And if she measures in a different basis entirely, the correlation structure changes.

Misconception 4: “Entanglement is fragile and instantly destroyed by any interaction.” Entanglement is sensitive to noise, but it does not vanish at the first hint of decoherence. Weak interactions degrade entanglement gradually. In fact, as we will see in the section on entanglement sudden death, concurrence can remain nonzero through many rounds of noise before dropping to zero at a finite threshold. Quantum error correction exploits this: by encoding logical qubits in entangled states of many physical qubits, you can protect entanglement against realistic noise rates.

Misconception 5: “More entanglement always means better quantum computation.” Entanglement is necessary for quantum speedup (any quantum computation that stays in product states throughout can be efficiently simulated classically). But entanglement alone is not sufficient. Some highly entangled states, such as stabilizer states, can be simulated efficiently on a classical computer via the Gottesman-Knill theorem. The power of quantum computing comes from the combination of entanglement, superposition, and carefully structured interference.

What Entanglement Is

Two qubits are entangled when their joint quantum state cannot be written as a product of two individual qubit states. A product state is one where each qubit has its own independent description:

import numpy as np

ket_0 = np.array([1, 0], dtype=complex)
ket_1 = np.array([0, 1], dtype=complex)
plus  = np.array([1/np.sqrt(2), 1/np.sqrt(2)], dtype=complex)

# Product state: qubit A in |+>, qubit B in |0>
# This CAN be written as tensor product of individual states
product_state = np.kron(plus, ket_0)
print("Product state:", np.round(product_state, 3))
# [0.707, 0., 0.707, 0.]

An entangled state cannot be factored this way:

# Bell state: (|00> + |11>) / sqrt(2)
# This CANNOT be written as tensor product of two single-qubit states
bell = np.array([1/np.sqrt(2), 0, 0, 1/np.sqrt(2)], dtype=complex)

# Check: try to factor it as [a, b] x [c, d]
# You need a*c = 1/sqrt(2), a*d = 0, b*c = 0, b*d = 1/sqrt(2)
# From a*d = 0 and a*c = 1/sqrt(2): a != 0, so d = 0
# But d = 0 contradicts b*d = 1/sqrt(2)
# Contradiction: the Bell state has no product form

Schmidt Decomposition

The Schmidt decomposition provides a systematic way to determine whether a pure two-qubit state is entangled and to quantify how much entanglement it contains.

For any two-qubit pure state |psi>, you can reshape its coefficient vector into a 2x2 matrix M and apply the singular value decomposition M = U S V†. The singular values (the diagonal entries of S) are the Schmidt coefficients. Their squares sum to 1 (by normalization). If only one Schmidt coefficient is nonzero, the state is a product state. If two or more are nonzero, the state is entangled.

The entanglement entropy quantifies the entanglement:

E = -sum_k (lambda_k^2) * log2(lambda_k^2)

where lambda_k are the Schmidt coefficients. For a maximally entangled Bell state, E = 1 ebit (one “bit” of entanglement). For a product state, E = 0.

import numpy as np

def schmidt_decomposition(psi_2q):
    """
    Compute the Schmidt decomposition of a 2-qubit pure state.
    
    Returns:
        schmidt_coeffs: array of Schmidt coefficients (singular values)
        U: left unitary (Schmidt basis for qubit A)
        Vh: right unitary conjugate transpose (Schmidt basis for qubit B)
    """
    # Reshape the 4-component state vector into a 2x2 matrix
    # Row index = qubit A, column index = qubit B
    M = psi_2q.reshape(2, 2)
    
    # Singular value decomposition
    U, singular_values, Vh = np.linalg.svd(M)
    
    return singular_values, U, Vh

def entanglement_entropy(schmidt_coeffs):
    """
    Compute the entanglement entropy from Schmidt coefficients.
    E = -sum_k lambda_k^2 * log2(lambda_k^2)
    """
    probs = schmidt_coeffs**2
    # Filter out zero probabilities to avoid log(0)
    probs = probs[probs > 1e-15]
    return -np.sum(probs * np.log2(probs))

# --- Test on several states ---

# 1. Maximally entangled Bell state: (|00> + |11>) / sqrt(2)
bell = np.array([1/np.sqrt(2), 0, 0, 1/np.sqrt(2)], dtype=complex)
coeffs, U, Vh = schmidt_decomposition(bell)
print(f"Bell state Schmidt coefficients: {np.round(coeffs, 4)}")
print(f"Bell state entanglement entropy: {entanglement_entropy(coeffs):.4f} ebits")
# Coefficients: [0.7071, 0.7071], Entropy: 1.0000 ebits

# 2. Partially entangled state: cos(pi/8)|00> + sin(pi/8)|11>
theta = np.pi / 8
partial = np.array([np.cos(theta), 0, 0, np.sin(theta)], dtype=complex)
coeffs, U, Vh = schmidt_decomposition(partial)
print(f"\nPartially entangled Schmidt coefficients: {np.round(coeffs, 4)}")
print(f"Partially entangled entropy: {entanglement_entropy(coeffs):.4f} ebits")
# Coefficients: [0.9239, 0.3827], Entropy: ~0.6009 ebits

# 3. Product state: |+> tensor |0>
product = np.kron(
    np.array([1/np.sqrt(2), 1/np.sqrt(2)], dtype=complex),
    np.array([1, 0], dtype=complex)
)
coeffs, U, Vh = schmidt_decomposition(product)
print(f"\nProduct state Schmidt coefficients: {np.round(coeffs, 4)}")
print(f"Product state entanglement entropy: {entanglement_entropy(coeffs):.4f} ebits")
# Coefficients: [1.0, 0.0], Entropy: 0.0000 ebits

The Four Bell States

The four Bell states form a complete orthonormal basis for the two-qubit Hilbert space. Each is maximally entangled (entanglement entropy = 1 ebit):

Name	State
Phi+	(
Phi-	(
Psi+	(
Psi-	(

Creating All Four Bell States in Qiskit

Each Bell state can be prepared with a short circuit starting from a computational basis state:

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

def make_bell_circuit(name):
    """
    Create a circuit that prepares the specified Bell state.
    
    Phi+ : H on q0, CNOT(q0, q1) applied to |00>
    Phi- : X on q0, H on q0, CNOT(q0, q1) applied to |00>
           (equivalently: H-CNOT on |00>, then Z on q0)
    Psi+ : X on q1, H on q0, CNOT(q0, q1) applied to |00>
           (equivalently: H-CNOT on |01>)
    Psi- : X on q1, X on q0, H on q0, CNOT(q0, q1) applied to |00>
           (equivalently: H-CNOT on |11>)
    """
    qc = QuantumCircuit(2)
    
    if name == "Phi+":
        qc.h(0)
        qc.cx(0, 1)
    elif name == "Phi-":
        qc.h(0)
        qc.cx(0, 1)
        qc.z(0)
    elif name == "Psi+":
        qc.x(1)
        qc.h(0)
        qc.cx(0, 1)
    elif name == "Psi-":
        qc.x(1)
        qc.h(0)
        qc.cx(0, 1)
        qc.z(0)
    
    qc.save_statevector()
    return qc

# Expected Bell state vectors (Qiskit uses little-endian qubit ordering)
# In Qiskit, the state vector index is |q1 q0>, so:
# |00> -> index 0, |01> -> index 1, |10> -> index 2, |11> -> index 3
expected = {
    "Phi+": np.array([1, 0, 0, 1]) / np.sqrt(2),
    "Phi-": np.array([1, 0, 0, -1]) / np.sqrt(2),
    "Psi+": np.array([0, 1, 1, 0]) / np.sqrt(2),
    "Psi-": np.array([0, 1, -1, 0]) / np.sqrt(2),
}

sim = AerSimulator(method="statevector")

for name in ["Phi+", "Phi-", "Psi+", "Psi-"]:
    qc = make_bell_circuit(name)
    result = sim.run(qc).result()
    sv = np.array(result.get_statevector(qc))
    
    # Check that the simulated state matches the expected Bell state
    # (up to a global phase)
    overlap = abs(np.dot(sv.conj(), expected[name]))
    print(f"{name}: statevector = {np.round(sv, 4)},  overlap = {overlap:.4f}")
    print(qc.draw())
    print()

Distinguishing Bell States by Measurement

A single measurement in the computational basis cannot distinguish all four Bell states. To perform a complete Bell state measurement, apply the inverse of the preparation circuit (CNOT followed by H) and then measure:

from qiskit import QuantumCircuit
from qiskit_aer import AerSimulator

def bell_measurement_circuit(bell_name):
    """Prepare a Bell state, then undo it to identify which one it is."""
    qc = QuantumCircuit(2, 2)
    
    # --- Prepare the Bell state ---
    if bell_name == "Phi+":
        qc.h(0); qc.cx(0, 1)
    elif bell_name == "Phi-":
        qc.h(0); qc.cx(0, 1); qc.z(0)
    elif bell_name == "Psi+":
        qc.x(1); qc.h(0); qc.cx(0, 1)
    elif bell_name == "Psi-":
        qc.x(1); qc.h(0); qc.cx(0, 1); qc.z(0)
    
    qc.barrier()
    
    # --- Bell measurement: reverse the Bell creation ---
    qc.cx(0, 1)
    qc.h(0)
    qc.measure([0, 1], [0, 1])
    return qc

sim = AerSimulator()

# Each Bell state maps to a unique 2-bit outcome
# Phi+ -> 00, Phi- -> 01 (phase flip on q0), Psi+ -> 10, Psi- -> 11
for name in ["Phi+", "Phi-", "Psi+", "Psi-"]:
    qc = bell_measurement_circuit(name)
    result = sim.run(qc, shots=1024).result()
    counts = result.get_counts()
    print(f"{name} -> measured: {counts}")

Classical Correlations vs. Quantum Correlations: the CHSH Inequality

Classical correlations can be explained by shared prior information. If you put one red ball and one blue ball in separate boxes and ship them to opposite sides of the planet, opening one box tells you the color of the other. But the colors were determined when you packed the boxes.

Quantum entanglement is different. Bell’s theorem, proved in 1964, shows that the correlations in entangled quantum states are stronger than any classical system could produce, regardless of what hidden information the particles might carry.

The CHSH inequality is a testable form of Bell’s theorem. For any classical hidden-variable theory, the CHSH parameter S satisfies:

|S| = |E(a,b) + E(a,b') + E(a',b) - E(a',b')| <= 2

where E(x,y) is the correlation when Alice measures along direction x and Bob measures along direction y. Quantum mechanics predicts that entangled states can achieve S = 2*sqrt(2), known as the Tsirelson bound.

Full CHSH Computation

For the Bell state |Phi+> = (|00> + |11>)/sqrt(2), the correlation function E(a,b) when Alice measures her qubit rotated by angle a from the Z axis and Bob measures his qubit rotated by angle b is:

E(a, b) = <Phi+| (R(a) tensor R(b))† (Z tensor Z) (R(a) tensor R(b)) |Phi+>

where R(theta) = exp(-i * theta/2 * Y) is a rotation about the Y axis.

import numpy as np

# Pauli matrices
I2 = np.eye(2, dtype=complex)
Z = np.array([[1, 0], [0, -1]], dtype=complex)
Y = np.array([[0, -1j], [1j, 0]], dtype=complex)

def Ry(theta):
    """Rotation about Y axis by angle theta."""
    return np.cos(theta/2) * I2 - 1j * np.sin(theta/2) * Y

def correlator(psi, angle_a, angle_b):
    """
    Compute E(a, b) = <psi| (Za tensor Zb) |psi>
    where Za = R(-a) Z R(a) is Z rotated by angle a,
    and Zb = R(-b) Z R(b) is Z rotated by angle b.
    """
    # Rotated measurement operators
    Za = Ry(-angle_a) @ Z @ Ry(angle_a)
    Zb = Ry(-angle_b) @ Z @ Ry(angle_b)
    
    # Tensor product of the two operators
    ZaZb = np.kron(Za, Zb)
    
    # Expectation value
    return np.real(psi.conj() @ ZaZb @ psi)

# Bell state |Phi+>
phi_plus = np.array([1/np.sqrt(2), 0, 0, 1/np.sqrt(2)], dtype=complex)

# Optimal CHSH angles
# Alice: a = 0, a' = pi/2
# Bob:   b = pi/4, b' = -pi/4
a, a_prime = 0, np.pi/2
b, b_prime = np.pi/4, -np.pi/4

E_ab   = correlator(phi_plus, a, b)
E_ab2  = correlator(phi_plus, a, b_prime)
E_a2b  = correlator(phi_plus, a_prime, b)
E_a2b2 = correlator(phi_plus, a_prime, b_prime)

S_quantum = E_ab + E_ab2 + E_a2b - E_a2b2

print("Quantum CHSH:")
print(f"  E(a,  b)  = {E_ab:.4f}")
print(f"  E(a,  b') = {E_ab2:.4f}")
print(f"  E(a', b)  = {E_a2b:.4f}")
print(f"  E(a', b') = {E_a2b2:.4f}")
print(f"  S = {S_quantum:.4f}")
print(f"  2*sqrt(2) = {2*np.sqrt(2):.4f}")
# S = 2.8284 = 2*sqrt(2), violating the classical bound of 2

# --- Classical simulation: correlated coins ---
# Best classical strategy: Alice and Bob use a shared random variable lambda
# and deterministic response functions A(a, lambda), B(b, lambda) in {+1, -1}
rng = np.random.default_rng(42)
n_trials = 1_000_000

# Shared random variable: uniform on [0, 2*pi)
lambdas = rng.uniform(0, 2*np.pi, n_trials)

# Deterministic strategy: A(a, lam) = sign(cos(a - lam)), B(b, lam) = sign(cos(b - lam))
def classical_response(angle, lam):
    return np.sign(np.cos(angle - lam))

A_a    = classical_response(a, lambdas)
A_a2   = classical_response(a_prime, lambdas)
B_b    = classical_response(b, lambdas)
B_b2   = classical_response(b_prime, lambdas)

S_classical = (
    np.mean(A_a * B_b) + np.mean(A_a * B_b2)
    + np.mean(A_a2 * B_b) - np.mean(A_a2 * B_b2)
)

print(f"\nClassical CHSH (Monte Carlo, {n_trials} trials):")
print(f"  S = {S_classical:.4f}")
print(f"  Classical bound: |S| <= 2")
# The classical S will be approximately <= 2

Concurrence and Entanglement Measures

For pure states, the Schmidt decomposition gives a complete picture of entanglement. For mixed states (statistical mixtures described by density matrices), we need different tools.

The concurrence C is an entanglement measure for two-qubit states, ranging from 0 (separable) to 1 (maximally entangled). For a density matrix rho:

1. Compute rho_tilde = (Y tensor Y) rho* (Y tensor Y), where rho* is the complex conjugate
2. Compute R = rho * rho_tilde
3. Find the eigenvalues of R in decreasing order: lambda_1, lambda_2, lambda_3, lambda_4
4. C = max(0, sqrt(lambda_1) - sqrt(lambda_2) - sqrt(lambda_3) - sqrt(lambda_4))

import numpy as np

Y = np.array([[0, -1j], [1j, 0]], dtype=complex)
YY = np.kron(Y, Y)

def concurrence(rho):
    """Compute the concurrence of a 2-qubit density matrix."""
    # rho_tilde = (Y x Y) rho* (Y x Y)
    rho_tilde = YY @ rho.conj() @ YY
    
    # R = rho @ rho_tilde
    R = rho @ rho_tilde
    
    # Eigenvalues of R (take real parts; imaginary parts are numerical noise)
    eigenvalues = np.real(np.linalg.eigvals(R))
    
    # Sort in decreasing order and take square roots
    eigenvalues = np.sort(np.maximum(eigenvalues, 0))[::-1]
    sqrt_eigs = np.sqrt(eigenvalues)
    
    return max(0.0, sqrt_eigs[0] - sqrt_eigs[1] - sqrt_eigs[2] - sqrt_eigs[3])

# --- Test 1: Bell state |Phi+> (C = 1) ---
phi_plus = np.array([1/np.sqrt(2), 0, 0, 1/np.sqrt(2)], dtype=complex)
rho_bell = np.outer(phi_plus, phi_plus.conj())
print(f"Bell state concurrence: {concurrence(rho_bell):.4f}")
# 1.0000

# --- Test 2: Product state |00> (C = 0) ---
ket_00 = np.array([1, 0, 0, 0], dtype=complex)
rho_product = np.outer(ket_00, ket_00.conj())
print(f"Product state concurrence: {concurrence(rho_product):.4f}")
# 0.0000

# --- Test 3: Werner state rho_W = p |Phi+><Phi+| + (1-p) I/4 ---
# Analytical result: C = max(0, (3p - 1) / 2)
# Entangled for p > 1/3
I4 = np.eye(4, dtype=complex)

for p in [0.2, 1/3, 0.5, 0.8, 1.0]:
    rho_werner = p * rho_bell + (1 - p) * I4 / 4
    C_numerical = concurrence(rho_werner)
    C_analytical = max(0, (3*p - 1) / 2)
    print(f"Werner p={p:.2f}: C_numerical={C_numerical:.4f}, C_analytical={C_analytical:.4f}")

Detecting Entanglement in Mixed States

For pure states, the Schmidt rank test is definitive. For mixed states, the partial transpose test (also called the Peres-Horodecki criterion) is widely used: if the partial transpose of a density matrix has any negative eigenvalue, the state is entangled.

import numpy as np

def partial_transpose_qubit1(rho_2q):
    """Partial transpose on qubit 1 (the second qubit)."""
    rho = rho_2q.reshape(2, 2, 2, 2)
    # Transpose indices of qubit 1: swap axis 1 and axis 3
    rho_pt = rho.transpose(0, 3, 2, 1).reshape(4, 4)
    return rho_pt

# Bell state density matrix
phi_plus = np.array([1/np.sqrt(2), 0, 0, 1/np.sqrt(2)], dtype=complex)
rho_bell = np.outer(phi_plus, phi_plus.conj())

rho_pt = partial_transpose_qubit1(rho_bell)
eigenvalues = np.linalg.eigvalsh(rho_pt)
print("Bell state, eigenvalues of partial transpose:", np.round(eigenvalues, 4))
# [-0.5,  0.5,  0.5,  0.5]
# Negative eigenvalue confirms entanglement

# Product state: no negative eigenvalues
ket_00 = np.array([1, 0, 0, 0], dtype=complex)
rho_product = np.outer(ket_00, ket_00.conj())
rho_pt_prod = partial_transpose_qubit1(rho_product)
eigenvalues_prod = np.linalg.eigvalsh(rho_pt_prod)
print("Product state, eigenvalues of partial transpose:", np.round(eigenvalues_prod, 4))
# [0., 0., 0., 1.] -- all non-negative, consistent with separability

For two-qubit systems, the partial transpose test is both necessary and sufficient: a two-qubit state is entangled if and only if its partial transpose has a negative eigenvalue. For higher-dimensional systems, the test is necessary but not sufficient (there exist entangled states with positive partial transpose, called “bound entangled” states).

Entanglement in Multipartite Systems

Moving beyond two qubits reveals a richer landscape. For three or more qubits, there are fundamentally different types of entanglement that cannot be converted into one another, even with local operations and classical communication.

GHZ and W States

The two most important classes of three-qubit entanglement are:

GHZ state: (|000> + |111>) / sqrt(2)

The GHZ (Greenberger-Horne-Zeilinger) state is the natural generalization of the Bell state to three qubits. It has a striking property: if you trace out (lose access to) any one qubit, the remaining two qubits are in a completely unentangled classical mixture. The entanglement is “all or nothing.”

W state: (|001> + |010> + |100>) / sqrt(3)

The W state distributes entanglement more democratically. If you trace out any one qubit, the remaining two qubits are still entangled. This makes W-state entanglement more robust against particle loss.

import numpy as np

ket_0 = np.array([1, 0], dtype=complex)
ket_1 = np.array([0, 1], dtype=complex)

def ket(bits):
    """Create a computational basis state from a bit string like '010'."""
    state = np.array([1], dtype=complex)
    for b in bits:
        state = np.kron(state, ket_0 if b == '0' else ket_1)
    return state

# GHZ state: (|000> + |111>) / sqrt(2)
ghz = (ket('000') + ket('111')) / np.sqrt(2)

# W state: (|001> + |010> + |100>) / sqrt(3)
w = (ket('001') + ket('010') + ket('100')) / np.sqrt(3)

def trace_out_qubit(psi_3q, qubit_to_remove):
    """
    Trace out one qubit from a 3-qubit pure state.
    Returns the 4x4 reduced density matrix of the remaining 2 qubits.
    
    qubit_to_remove: 0 (first), 1 (middle), or 2 (last)
    """
    # Full density matrix
    rho = np.outer(psi_3q, psi_3q.conj())
    
    # Reshape to (2,2,2, 2,2,2) for qubit indices
    rho_tensor = rho.reshape(2, 2, 2, 2, 2, 2)
    
    # Trace over the specified qubit
    # For qubit k, sum over axes k and k+3
    if qubit_to_remove == 0:
        # Trace over axes 0 and 3
        rho_reduced = np.trace(rho_tensor, axis1=0, axis2=3)  # shape (2,2,2,2)
    elif qubit_to_remove == 1:
        # Trace over axes 1 and 4
        rho_reduced = np.trace(rho_tensor, axis1=1, axis2=4)  # shape (2,2,2,2)
    elif qubit_to_remove == 2:
        # Trace over axes 2 and 5
        rho_reduced = np.trace(rho_tensor, axis1=2, axis2=5)  # shape (2,2,2,2)
    
    return rho_reduced.reshape(4, 4)

def purity(rho):
    """Purity = Tr(rho^2). Equal to 1 for pure states, < 1 for mixed states."""
    return np.real(np.trace(rho @ rho))

print("=== GHZ state: trace out each qubit ===")
for q in range(3):
    rho_reduced = trace_out_qubit(ghz, q)
    p = purity(rho_reduced)
    # Check if the reduced state is entangled via partial transpose
    rho_pt = rho_reduced.reshape(2,2,2,2).transpose(0,3,2,1).reshape(4,4)
    min_eig = np.min(np.linalg.eigvalsh(rho_pt))
    print(f"  Trace out qubit {q}: purity = {p:.4f}, min PT eigenvalue = {min_eig:.4f}")
# Purity = 0.5 for all (maximally mixed 2-qubit state)
# min PT eigenvalue >= 0 (no entanglement in the remaining pair)

print("\n=== W state: trace out each qubit ===")
for q in range(3):
    rho_reduced = trace_out_qubit(w, q)
    p = purity(rho_reduced)
    rho_pt = rho_reduced.reshape(2,2,2,2).transpose(0,3,2,1).reshape(4,4)
    min_eig = np.min(np.linalg.eigvalsh(rho_pt))
    print(f"  Trace out qubit {q}: purity = {p:.4f}, min PT eigenvalue = {min_eig:.4f}")
# Purity > 0.5 (not maximally mixed)
# min PT eigenvalue < 0 (remaining pair IS still entangled)

The GHZ and W states represent genuinely distinct entanglement classes. No amount of local operations and classical communication can convert a GHZ state into a W state or vice versa. This classification becomes even more complex for four or more qubits, where infinitely many distinct entanglement classes exist.

Monogamy of Entanglement

Entanglement obeys a fundamental constraint called monogamy: the more entangled qubit A is with qubit B, the less entangled it can be with qubit C. Quantitatively, for three qubits:

tau(A:BC) >= tau(A:B) + tau(A:C)

where tau is the tangle, defined as the square of the concurrence (tau = C^2). If A is maximally entangled with B (tangle = 1), it cannot share any entanglement with C.

This is crucial for quantum cryptography. In quantum key distribution (QKD), if Alice and Bob verify that their qubits are maximally entangled, the monogamy constraint guarantees that no eavesdropper Eve can share any entanglement with their qubits. Any eavesdropping necessarily disturbs the entanglement, which Alice and Bob can detect.

import numpy as np

# Reusing the concurrence function from above
Y = np.array([[0, -1j], [1j, 0]], dtype=complex)
YY = np.kron(Y, Y)

def concurrence(rho):
    """Compute the concurrence of a 2-qubit density matrix."""
    rho_tilde = YY @ rho.conj() @ YY
    R = rho @ rho_tilde
    eigenvalues = np.real(np.linalg.eigvals(R))
    eigenvalues = np.sort(np.maximum(eigenvalues, 0))[::-1]
    sqrt_eigs = np.sqrt(eigenvalues)
    return max(0.0, sqrt_eigs[0] - sqrt_eigs[1] - sqrt_eigs[2] - sqrt_eigs[3])

def ket(bits):
    """Computational basis state from bit string."""
    ket_0 = np.array([1, 0], dtype=complex)
    ket_1 = np.array([0, 1], dtype=complex)
    state = np.array([1], dtype=complex)
    for b in bits:
        state = np.kron(state, ket_0 if b == '0' else ket_1)
    return state

def partial_trace_to_pair(psi_3q, keep):
    """
    Trace out one qubit from a 3-qubit state, keeping the specified pair.
    keep: tuple of two qubit indices to keep, e.g. (0,1) keeps qubits 0 and 1.
    """
    rho = np.outer(psi_3q, psi_3q.conj()).reshape(2,2,2,2,2,2)
    remove = list(set([0,1,2]) - set(keep))[0]
    rho_red = np.trace(rho, axis1=remove, axis2=remove+3)
    return rho_red.reshape(4,4)

# GHZ state: (|000> + |111>) / sqrt(2)
ghz = (ket('000') + ket('111')) / np.sqrt(2)

# Compute pairwise tangles
rho_AB = partial_trace_to_pair(ghz, (0, 1))
rho_AC = partial_trace_to_pair(ghz, (0, 2))

C_AB = concurrence(rho_AB)
C_AC = concurrence(rho_AC)
tau_AB = C_AB**2
tau_AC = C_AC**2

# Tangle of A with BC (for a pure 3-qubit state, this equals
# the linear entropy of the reduced state of A)
rho_full = np.outer(ghz, ghz.conj()).reshape(2,2,2,2,2,2)
rho_A = np.trace(np.trace(rho_full, axis1=2, axis2=5), axis1=1, axis2=2)
# rho_A is a 2x2 matrix
tau_A_BC = 4 * np.real(np.linalg.det(rho_A))  # = 4 * det(rho_A) for a qubit

print("=== Monogamy for GHZ state ===")
print(f"tau(A:B)  = {tau_AB:.4f}")
print(f"tau(A:C)  = {tau_AC:.4f}")
print(f"tau(A:BC) = {tau_A_BC:.4f}")
print(f"tau(A:B) + tau(A:C) = {tau_AB + tau_AC:.4f}")
print(f"Monogamy: tau(A:BC) >= tau(A:B) + tau(A:C)?  {tau_A_BC >= tau_AB + tau_AC - 1e-10}")
# For GHZ: tau(A:B) = tau(A:C) = 0 (pairwise entanglement is zero!)
# But tau(A:BC) = 1 (A is maximally entangled with BC jointly)
# The "residual tangle" tau(A:BC) - tau(A:B) - tau(A:C) = 1
# is called the 3-tangle, a measure of genuine 3-party entanglement

Quantum Teleportation

Quantum teleportation transfers an unknown quantum state from Alice to Bob using a shared entangled pair and two classical bits of communication. No physical qubit travels from Alice to Bob, and no faster-than-light signaling occurs.

The protocol works as follows:

1. Alice and Bob share a Bell pair (qubits 1 and 2 in state |Phi+>).
2. Alice has an unknown state |psi> on qubit 0.
3. Alice performs a Bell measurement on qubits 0 and 1 (CNOT then H, then measure).
4. Alice sends her two measurement bits to Bob over a classical channel.
5. Bob applies corrections to qubit 2: X if bit 0 is 1, Z if bit 1 is 1.
6. Bob’s qubit 2 is now in state |psi>.

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

def teleportation_circuit(theta, phi):
    """
    Teleport the state cos(theta/2)|0> + e^(i*phi)*sin(theta/2)|1>
    from qubit 0 to qubit 2.
    """
    qc = QuantumCircuit(3, 2)
    
    # Step 0: Prepare the state to teleport on qubit 0
    qc.ry(theta, 0)
    qc.rz(phi, 0)
    
    # Step 1: Create Bell pair between qubit 1 (Alice) and qubit 2 (Bob)
    qc.h(1)
    qc.cx(1, 2)
    
    qc.barrier()
    
    # Step 2: Alice performs Bell measurement on qubits 0 and 1
    qc.cx(0, 1)
    qc.h(0)
    qc.measure(0, 0)
    qc.measure(1, 1)
    
    qc.barrier()
    
    # Step 3: Bob applies conditional corrections
    # X gate on qubit 2 if classical bit 1 is 1
    qc.x(2).c_if(1, 1)
    # Z gate on qubit 2 if classical bit 0 is 1
    qc.z(2).c_if(0, 1)
    
    return qc

# Choose a state to teleport: |psi> = cos(pi/6)|0> + sin(pi/6)|1>
theta = np.pi / 3  # Ry rotation angle
phi = np.pi / 5    # Rz phase angle

qc = teleportation_circuit(theta, phi)

# To verify teleportation, save the statevector after corrections
# We use a separate circuit without measurement for verification
qc_verify = QuantumCircuit(3)
qc_verify.ry(theta, 0)
qc_verify.rz(phi, 0)
qc_verify.h(1)
qc_verify.cx(1, 2)
qc_verify.cx(0, 1)
qc_verify.h(0)
# For verification, assume measurement gave 00 (no corrections needed)
qc_verify.save_statevector()

sim = AerSimulator(method="statevector")
result = sim.run(qc_verify).result()
sv = np.array(result.get_statevector())

# The original state on qubit 0
original = np.array([
    np.cos(theta/2),
    np.exp(1j * phi) * np.sin(theta/2)
], dtype=complex)

print("Teleportation circuit:")
print(qc.draw())
print(f"\nOriginal state |psi>: {np.round(original, 4)}")
print(f"\nFull 3-qubit statevector (before measurement):")
print(np.round(sv, 4))
print("\nAfter Alice measures 00, Bob's qubit 2 will be in state |psi>.")
print("After Alice measures 01, Bob applies X to recover |psi>.")
print("After Alice measures 10, Bob applies Z to recover |psi>.")
print("After Alice measures 11, Bob applies X then Z to recover |psi>.")

# Run the full circuit with measurement to show it works statistically
qc_test = teleportation_circuit(0, 0)  # Teleport |0>
qc_test.measure_all()
result = AerSimulator().run(qc_test, shots=1024).result()
print(f"\nTeleporting |0>: measurement results = {result.get_counts()}")

Superdense Coding

Superdense coding is the complement of teleportation. Instead of using entanglement and classical bits to send a quantum state, superdense coding uses entanglement and one qubit to send two classical bits.

The protocol:

1. Alice and Bob share a Bell pair |Phi+> (qubit 0 is Alice’s, qubit 1 is Bob’s).
2. Alice encodes two classical bits by applying a gate to her qubit:
   • 00: Identity (do nothing)
   • 01: X gate (bit flip)
   • 10: Z gate (phase flip)
   • 11: X then Z (both flips)
3. Alice sends her qubit to Bob.
4. Bob decodes by applying CNOT(0,1) then H(0), and measuring both qubits.

from qiskit import QuantumCircuit
from qiskit_aer import AerSimulator

def superdense_coding(bits):
    """
    Superdense coding: encode 2 classical bits using 1 qubit transmission.
    bits: string "00", "01", "10", or "11"
    """
    qc = QuantumCircuit(2, 2)
    
    # Step 1: Create shared Bell pair
    qc.h(0)
    qc.cx(0, 1)
    
    qc.barrier()
    
    # Step 2: Alice encodes her 2 bits on qubit 0
    if bits == "00":
        pass  # Identity
    elif bits == "01":
        qc.x(0)
    elif bits == "10":
        qc.z(0)
    elif bits == "11":
        qc.x(0)
        qc.z(0)
    
    qc.barrier()
    
    # Step 3: Bob decodes (Alice sends qubit 0 to Bob)
    qc.cx(0, 1)
    qc.h(0)
    qc.measure([0, 1], [0, 1])
    
    return qc

sim = AerSimulator()

for bits in ["00", "01", "10", "11"]:
    qc = superdense_coding(bits)
    result = sim.run(qc, shots=1024).result()
    counts = result.get_counts()
    print(f"Encoded: {bits}  ->  Decoded: {counts}")
    # Each case should decode to exactly the encoded bits

print("\nCircuit for encoding '11':")
print(superdense_coding("11").draw())

Entanglement Swapping

Entanglement swapping creates entanglement between two particles that have never directly interacted. This is the foundation of quantum repeaters, which will enable long-distance quantum communication.

Setup: There are four qubits. Qubits 0 and 1 form Bell pair 1. Qubits 2 and 3 form Bell pair 2. Initially, qubits 0 and 3 share no entanglement.

Protocol: Perform a Bell state measurement on qubits 1 and 2 (one from each pair). After this measurement, qubits 0 and 3, which never interacted, become entangled.

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

# Entanglement swapping circuit
qc = QuantumCircuit(4, 2)

# Create Bell pair 1: qubits 0 and 1
qc.h(0)
qc.cx(0, 1)

# Create Bell pair 2: qubits 2 and 3
qc.h(2)
qc.cx(2, 3)

qc.barrier()

# Bell state measurement on qubits 1 and 2
qc.cx(1, 2)
qc.h(1)
qc.measure(1, 0)
qc.measure(2, 1)

qc.barrier()

# Conditional corrections on qubit 3 (based on measurement of qubits 1 and 2)
qc.x(3).c_if(1, 1)  # If qubit 2 measured 1, apply X to qubit 3
qc.z(3).c_if(0, 1)  # If qubit 1 measured 1, apply Z to qubit 3

print("Entanglement swapping circuit:")
print(qc.draw())

# Verify using statevector simulation (without mid-circuit measurement)
qc_verify = QuantumCircuit(4)
qc_verify.h(0)
qc_verify.cx(0, 1)
qc_verify.h(2)
qc_verify.cx(2, 3)
qc_verify.cx(1, 2)
qc_verify.h(1)
qc_verify.save_statevector()

sim = AerSimulator(method="statevector")
result = sim.run(qc_verify).result()
sv = np.array(result.get_statevector())

# Compute reduced density matrix of qubits 0 and 3
rho_full = np.outer(sv, sv.conj()).reshape(2,2,2,2, 2,2,2,2)
# Trace out qubits 1 and 2 (axes 1,2 and 5,6)
rho_03 = np.trace(np.trace(rho_full, axis1=2, axis2=6), axis1=1, axis2=4)
rho_03 = rho_03.reshape(4, 4)

# Check entanglement via partial transpose
rho_pt = rho_03.reshape(2,2,2,2).transpose(0,3,2,1).reshape(4,4)
min_eig = np.min(np.linalg.eigvalsh(rho_pt))
purity_03 = np.real(np.trace(rho_03 @ rho_03))

print(f"\nReduced state of qubits 0 and 3:")
print(f"  Purity: {purity_03:.4f}")
print(f"  Min eigenvalue of partial transpose: {min_eig:.4f}")
print(f"  Entangled? {min_eig < -1e-10}")
# After tracing out qubits 1 and 2 without conditioning on their measurement,
# qubits 0 and 3 are in a mixed but entangled state.
# Conditioned on a specific measurement outcome, they are in a pure Bell state.

Entanglement in Quantum Algorithms

Entanglement is not just a theoretical curiosity; it is an essential ingredient in quantum algorithms that achieve speedups over classical computation.

Grover’s Algorithm

In Grover’s search algorithm, the oracle marks the target state with a phase flip, and the diffusion operator amplifies the amplitude of that state. The diffusion operator creates entanglement across all qubits through a multi-controlled phase gate. Without entanglement (without the CNOT and multi-qubit gates), each qubit evolves independently, and the interference that concentrates amplitude on the target state cannot occur.

Here is a demonstration: a working 2-qubit Grover circuit that finds |11>, compared with a broken version where the entangling gates are removed.

from qiskit import QuantumCircuit
from qiskit_aer import AerSimulator

def grover_2qubit_working():
    """Grover's algorithm to find |11> among 4 items (1 iteration)."""
    qc = QuantumCircuit(2, 2)
    
    # Initialize uniform superposition
    qc.h([0, 1])
    
    # Oracle: mark |11> with a phase flip (CZ gate)
    qc.cz(0, 1)  # This is an entangling gate
    
    # Diffusion operator
    qc.h([0, 1])
    qc.z([0, 1])
    qc.cz(0, 1)  # Another entangling gate
    qc.h([0, 1])
    
    qc.measure([0, 1], [0, 1])
    return qc

def grover_2qubit_broken():
    """Same circuit but with entangling (CZ) gates removed."""
    qc = QuantumCircuit(2, 2)
    
    # Initialize uniform superposition
    qc.h([0, 1])
    
    # Oracle: CZ removed (no entanglement)
    # qc.cz(0, 1)  # REMOVED
    
    # Diffusion operator: CZ removed
    qc.h([0, 1])
    qc.z([0, 1])
    # qc.cz(0, 1)  # REMOVED
    qc.h([0, 1])
    
    qc.measure([0, 1], [0, 1])
    return qc

sim = AerSimulator()

print("Working Grover (with entanglement):")
result = sim.run(grover_2qubit_working(), shots=1024).result()
print(f"  Counts: {result.get_counts()}")
# Should find |11> with 100% probability

print("\nBroken Grover (entangling gates removed):")
result = sim.run(grover_2qubit_broken(), shots=1024).result()
print(f"  Counts: {result.get_counts()}")
# Uniform distribution, no speedup at all

Shor’s Algorithm

In Shor’s factoring algorithm, the quantum Fourier transform (QFT) creates entanglement across all n qubits of the register simultaneously. This entanglement allows the QFT to extract the period of a modular exponentiation function, which is the key step in factoring. Without the controlled rotation gates in the QFT (which are entangling gates), each qubit would carry only its own Fourier coefficient, and the period-finding interference pattern would not form.

Quantum Error Correction

Quantum error correction encodes a single logical qubit into an entangled state of multiple physical qubits. The smallest code that corrects arbitrary single-qubit errors, the 5-qubit code, uses a highly entangled state of 5 physical qubits. Syndrome measurements project the physical qubits onto a subspace that reveals which error occurred, all without measuring (and thus disturbing) the logical information. This is only possible because the syndrome operators act on the entangled state as a whole, extracting error information while leaving the encoded state intact.

Decoherence and Entanglement Sudden Death

In real quantum hardware, noise gradually degrades entanglement. A remarkable phenomenon called entanglement sudden death shows that entanglement can vanish completely at a finite noise level, rather than decaying smoothly to zero.

To see this, start with a Bell state and apply a depolarizing channel to both qubits at each step. The depolarizing channel with parameter p replaces the state with the maximally mixed state with probability p:

rho -> (1 - p) * rho + p * I/4

After each step, compute the concurrence. You will find that the concurrence drops to exactly zero at a finite step count.

import numpy as np
import warnings
warnings.filterwarnings("ignore")

Y = np.array([[0, -1j], [1j, 0]], dtype=complex)
YY = np.kron(Y, Y)
I4 = np.eye(4, dtype=complex)

def concurrence(rho):
    """Compute the concurrence of a 2-qubit density matrix."""
    rho_tilde = YY @ rho.conj() @ YY
    R = rho @ rho_tilde
    eigenvalues = np.real(np.linalg.eigvals(R))
    eigenvalues = np.sort(np.maximum(eigenvalues, 0))[::-1]
    sqrt_eigs = np.sqrt(eigenvalues)
    return max(0.0, sqrt_eigs[0] - sqrt_eigs[1] - sqrt_eigs[2] - sqrt_eigs[3])

def depolarizing_channel_2q(rho, p):
    """Apply 2-qubit depolarizing channel with error rate p."""
    return (1 - p) * rho + p * I4 / 4

# Start with Bell state |Phi+>
phi_plus = np.array([1/np.sqrt(2), 0, 0, 1/np.sqrt(2)], dtype=complex)
rho_init = np.outer(phi_plus, phi_plus.conj())

n_steps = 80

for p in [0.05, 0.10, 0.20]:
    rho = rho_init.copy()
    concurrences = [concurrence(rho)]
    death_step = None
    
    for step in range(1, n_steps + 1):
        rho = depolarizing_channel_2q(rho, p)
        C = concurrence(rho)
        concurrences.append(C)
        if C == 0 and death_step is None:
            death_step = step
    
    print(f"p = {p:.2f}: entanglement dies at step {death_step}")
    # Print concurrence at a few checkpoints
    checkpoints = [0, 5, 10, 20, 40, 60, 80]
    for s in checkpoints:
        if s <= n_steps:
            print(f"  Step {s:3d}: C = {concurrences[s]:.4f}")
    print()

# To create a proper plot, uncomment below (requires matplotlib)
# import matplotlib.pyplot as plt
# fig, ax = plt.subplots()
# for p in [0.05, 0.10, 0.20]:
#     rho = rho_init.copy()
#     cs = [concurrence(rho)]
#     for _ in range(n_steps):
#         rho = depolarizing_channel_2q(rho, p)
#         cs.append(concurrence(rho))
#     ax.plot(range(n_steps + 1), cs, label=f"p = {p:.2f}")
# ax.set_xlabel("Noise steps")
# ax.set_ylabel("Concurrence")
# ax.set_title("Entanglement Sudden Death")
# ax.legend()
# ax.set_ylim(0, 1.05)
# plt.savefig("entanglement_sudden_death.png", dpi=150)
# plt.show()

The key observation: concurrence reaches exactly zero at a finite number of steps. It does not asymptotically approach zero. This is entanglement sudden death, first predicted by Yu and Eberly in 2004 and experimentally confirmed in 2007. The practical implication is that quantum error correction must act before entanglement dies, not merely slow its decay.

Entanglement as a Resource

Entanglement is treated as a fungible resource in quantum information theory, much like energy in thermodynamics. The standard unit is the ebit: one ebit is the entanglement contained in a single Bell pair.

Resource accounting makes the costs of quantum protocols explicit:

Protocol	Entanglement cost	Classical communication	Quantum communication	Result
Teleportation	1 ebit consumed	2 classical bits sent	0 qubits	1 qubit state transferred
Superdense coding	1 ebit consumed	0 classical bits	1 qubit sent	2 classical bits transferred
Entanglement swapping	2 ebits consumed	2 classical bits sent	0 qubits	1 ebit created (remote pair)

Entanglement distillation takes many copies of a noisy (weakly entangled) state and produces fewer copies of a nearly perfect Bell pair. The reverse process, entanglement dilution, converts Bell pairs into copies of a target entangled state. The rate of interconversion is governed by the entanglement entropy, connecting the abstract measure to operational protocols.

Key Takeaways

• Entanglement means a multi-qubit state cannot be written as a product of individual qubit states
• The Schmidt decomposition quantifies entanglement for pure states via entanglement entropy
• Entangled correlations are stronger than any classical correlation, as proven by Bell’s theorem and the CHSH inequality (quantum bound 2*sqrt(2) vs classical bound 2)
• Entanglement cannot be used for faster-than-light communication
• Concurrence quantifies entanglement for mixed states, with the Werner state providing a clean parametric example
• Entanglement enables quantum teleportation, superdense coding, entanglement swapping, and quantum algorithms
• Multipartite entanglement (GHZ vs W states) reveals fundamentally different entanglement classes
• Monogamy constrains how entanglement can be shared, underpinning quantum cryptography
• Entanglement sudden death means noise destroys entanglement at a finite threshold, not asymptotically
• Entanglement is necessary but not sufficient for quantum computational advantage

Next Steps

• Quantum Teleportation Explained, the classic application of Bell pairs
• Quantum Measurement Explained, explaining how measurement affects entangled systems
• How Quantum Algorithms Work, covering the role of entanglement in quantum speedup