Level 3 · Quantum Computing: Entanglement

Level 2 had one qubit: two amplitudes, 2 by 2 gates, and the Born rule. This level puts qubits side by side in a register. Three qubits have eight bit strings, so a register of three qubits is a list of eight amplitudes, and n qubits need 2^n. Your simulator grows up here: it joins single qubits into one state with the Kronecker product, applies a gate to one qubit of many without ever building the giant matrix, runs the two-qubit gate CNOT, and measures one qubit while the others carry on.

Then comes the new idea. Most states of two qubits cannot be split into a state for each qubit: the pair has a state, but neither qubit has one of its own. That is entanglement. You will learn to recognise it (a 2 by 2 table whose determinant is not 0), to measure it (singular values), and to use it in the first protocols: teleportation, superdense coding, and the BB84 key exchange, which needs no entanglement but lives on the fact that measuring disturbs.

The paths meet here. The bit tricks of the Algorithms path at Level 5 (<<, &, |, ^) are how a simulator finds one qubit inside an index. The Statistics path's Level 3 has independence and conditional probability from tables of counts: a product state is independence, and measuring one qubit is conditioning on what it showed. The Signals, Systems & Control path meets the singular values and the rank of a matrix at Level 6; here they say how entangled two parts of a register are.

A register, its bit order and bit masks

A register of n qubits has one basis state for each string of n bits, written |b_0 b_1 \dots b_{n-1}\rangle, and its state is a list of 2^n complex amplitudes. Measuring every qubit gives basis state number i with probability |a_i|^2, the Born rule of Level 2.

The bit order is Level 2's outcome labels: |b_0 b_1 \dots b_{n-1}\rangle has index int("b0b1...", 2), so qubit 0 is the leftmost bit, the most significant one, and qubit n-1 is the rightmost. |011\rangle is index 3, and in it qubit 0 reads 0 while qubits 1 and 2 read 1.

Read along the second table: qubit 0 reads 1 in the right half (indices 4 to 7), qubit 2 reads 1 at every odd index, and qubit 1 in between, in blocks of two.

A simulator rarely needs the whole string format(i, f"0{n}b"); it needs one qubit's bit of an index, fast. Qubit q of n sits n - 1 - q places from the right, so its mask, the number with only that bit set, is 1 << (n - 1 - q). For 3 qubits the masks are 4, 2 and 1, written 100, 010 and 001. Four operations do everything:

• read qubit q: (i >> (n - 1 - q)) & 1, or test i & mask (nonzero means 1);
• flip it: i ^ mask;
• set it to 1: i | mask;
• clear it to 0: i & ~mask.

n = 4
i = int("1101", 2)
for q in range(n):
    s = n - 1 - q
    mask = 1 << s
    print("qubit", q, "mask", format(mask, "04b"), "reads", (i >> s) & 1,
          " flipped", format(i ^ mask, "04b"), " set", format(i | mask, "04b"), " cleared", format(i & ~mask, "04b"))
print("right:", [(i >> (n - 1 - q)) & 1 for q in range(n)], " wrong:", [(i >> q) & 1 for q in range(n)])
print("qubit 1 of 3 reads 1 at", [k for k in range(8) if k & 0b010])

qubit 0 mask 1000 reads 1  flipped 0101  set 1101  cleared 0101
qubit 1 mask 0100 reads 1  flipped 1001  set 1101  cleared 1001
qubit 2 mask 0010 reads 0  flipped 1111  set 1111  cleared 1101
qubit 3 mask 0001 reads 1  flipped 1100  set 1101  cleared 1100
right: [1, 1, 0, 1]  wrong: [1, 0, 1, 1]
qubit 1 of 3 reads 1 at [2, 3, 6, 7]

The wrong list reads the bits from the right, so it numbers the qubits backwards: here it claims qubit 1 reads 0. (Some toolkits do number the qubits from the right, as Level 2's aside said; never mix the two.) The last line lists the indices where one qubit reads 1, which is the first step of every measurement this level.

Common mistakes: 1 << q instead of 1 << (n - 1 - q). Forgetting that n qubits have 2^n amplitudes, not n: recover n from a state with len(state).bit_length() - 1. Writing i & mask == 1: when the bit is set, i & mask is the mask itself (2, 4, ...), not 1, so the test fails for every qubit but the last. Use != 0, or shift first. Using ~ to flip: ~i flips every bit and makes a negative number.

Practice: Which Lamps Are Lit on the Register?, Asking One Qubit Out of Many

Product states: kron, qubit 0 first

Prepare qubit 0 in [a, b] and qubit 1 in [c, d], separately. The register is in the product state whose amplitude for |b_0 b_1\rangle is qubit 0's amplitude for b_0 times qubit 1's amplitude for b_1:

That is the Kronecker product of Level 1, with qubit 0 as the left factor. numpy has it built in as np.kron, in the same order:

import numpy as np
def show(v):
    v = np.round(np.asarray(v, dtype=complex), 4) + 0       # rounded; + 0 turns -0.0 into 0.0
    return [x.real if x.imag == 0 else x for x in v.tolist()]

q0 = [0, 1]                  # |1>
q1 = [0.8, 0.6j]
print(show(np.kron(q0, q1)))
print(show(np.kron(q1, q0)))
one, zero = [0, 1], [1, 0]
print(show(np.kron(np.kron(one, zero), one)), "index", int("101", 2))
print(round(float(np.sum(np.abs(np.kron(q0, q1)) ** 2)), 6))

[0.0, 0.0, 0.8, 0.6j]
[0.0, 0.8, 0.0, 0.6j]
[0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0] index 5
1.0

Qubit 0 reads 1, so the amplitudes sit at indices 2 and 3 (10 and 11). The swapped order puts them at 1 and 3, as if qubit 1 were the one reading 1. For three qubits, multiply on: q_0 \otimes q_1 \otimes q_2 is np.kron(np.kron(q0, q1), q2), and the grouping does not matter, only the order. |1\rangle \otimes |0\rangle \otimes |1\rangle is |101\rangle, index 5. Weights multiply, so a product of unit vectors is a unit vector.

Common mistakes: the factors in reverse order, which makes qubit 0 the rightmost bit. np.outer(q0, q1), which gives the same four numbers as a 2 by 2 table instead of one list (useful later, but not a state). Adding the qubit states instead of multiplying.

Practice: Qubits Prepared One by One

A gate on one qubit, as a full matrix

A gate G on qubit k that leaves every other qubit alone is the Kronecker product of n matrices, in the order of the qubits, with G in place k and the 2 by 2 identity I everywhere else:

The reason is a rule of the Kronecker product, (A \otimes B)(u \otimes v) = (Au) \otimes (Bv): each factor acts on its own qubit. Here is X on qubit 1 of 3, applied to |001\rangle:

import numpy as np
I = np.eye(2)
X = np.array([[0, 1], [1, 0]])
state = np.zeros(8); state[1] = 1                 # |001>
M = np.kron(np.kron(I, X), I)                     # X on qubit 1
print(M.shape, int(np.argmax(M @ state)), format(int(np.argmax(M @ state)), "03b"))
wrong = np.kron(np.kron(I, I), X)                 # X on qubit 2 instead
print(format(int(np.argmax(wrong @ state)), "03b"))
print(M.astype(int).tolist()[:4])

(8, 8) 3 011
000
[[0, 0, 1, 0, 0, 0, 0, 0], [0, 0, 0, 1, 0, 0, 0, 0], [1, 0, 0, 0, 0, 0, 0, 0], [0, 1, 0, 0, 0, 0, 0, 0]]

Flipping qubit 1 of 001 gives 011, index 3. The matrix with X last flips qubit 2 instead. The first rows of M show its shape: it is mostly zeros, with one 1 per row.

This route is correct, and hopeless for big registers. The state has 2^n amplitudes; the matrix has 2^n \times 2^n = 4^n entries, at 16 bytes each for complex numbers:

def size(nbytes):
    for unit in ("B", "kB", "MB", "GB", "TB", "PB"):
        if nbytes < 1000:
            return f"{nbytes:.3g} {unit}"
        nbytes /= 1000
    return f"{nbytes:.3g} EB"

for n in (3, 10, 18, 30):
    print(f"{n:2d} qubits: state {size(16 * 2 ** n):>8}, full matrix {size(16 * 4 ** n):>8}")

 3 qubits: state    128 B, full matrix  1.02 kB
10 qubits: state  16.4 kB, full matrix  16.8 MB
18 qubits: state  4.19 MB, full matrix   1.1 TB
30 qubits: state  17.2 GB, full matrix  18.4 EB

At 18 qubits the state is a few megabytes and the matrix a terabyte, almost all of it zeros. The next two sections apply a gate the way the matrix does, without the matrix.

Common mistakes: putting G in place n - 1 - k. Building np.kron(G, np.eye(2 ** (n - 1))) for every k, which always acts on qubit 0. Building a full matrix in a loop over many gates of a large register.

Practice: Padding a Gate with Identities

Pairs of amplitudes

Look at what G on qubit k does to one index. It changes only qubit k's bit, so the new amplitude at i depends only on the old amplitudes at i and at i with that bit flipped. The indices fall into pairs that differ only in qubit k: i0, with qubit k reading 0, and i1 = i0 | mask. On each pair G acts exactly as it acts on one qubit:

Here are the pairs for qubit 1 of 3, and the full matrix of H on that qubit agreeing with them: each row has its nonzero entries exactly in the columns of its pair.

import numpy as np
n, k = 3, 1
mask = 1 << (n - 1 - k)
pairs = [(i0, i0 | mask) for i0 in range(2 ** n) if not i0 & mask]
print(pairs)
print([(format(a, "03b"), format(b, "03b")) for a, b in pairs])
H = np.array([[1, 1], [1, -1]]) / np.sqrt(2)
M = np.kron(np.kron(np.eye(2), H), np.eye(2))
print([np.flatnonzero(row).tolist() for row in M])

[(0, 2), (1, 3), (4, 6), (5, 7)]
[('000', '010'), ('001', '011'), ('100', '110'), ('101', '111')]
[[0, 2], [1, 3], [0, 2], [1, 3], [4, 6], [5, 7], [4, 6], [5, 7]]

Each pair is visited once, from its i0; 2 ** (n - 1) pairs cover all 2 ** n indices. Work per gate drops from 4^n to 2^n.

To watch amplitudes move, push the list 0, 1, ..., 7 through X on qubit 1 (a gate is linear, so it can be applied to any list, normalised or not). X swaps the two members of every pair. Read from the old list and write into a new one:

n, k = 3, 1
mask = 1 << (n - 1 - k)
old = list(range(8))
new = list(old)
for i0 in range(2 ** n):
    if not i0 & mask:
        i1 = i0 | mask
        new[i0], new[i1] = old[i1], old[i0]
print(new)

bad = list(old)                 # every index, in place: each pair is swapped twice
for i in range(2 ** n):
    j = i ^ mask
    bad[i], bad[j] = bad[j], bad[i]
print(bad)

[2, 3, 0, 1, 6, 7, 4, 5]
[0, 1, 2, 3, 4, 5, 6, 7]

The second loop visits each pair from both ends, so it swaps it back and changes nothing. For a general gate there is a second trap: both new values need both old values. Here H acts on one pair holding (1, 0):

import math
s = 1 / math.sqrt(2)
x, y = 1.0, 0.0
x = s * x + s * y               # overwrites x ...
y = s * x - s * y               # ... and then uses the new x
print(round(x, 4), round(y, 4))

x, y = 1.0, 0.0
x, y = s * x + s * y, s * x - s * y      # the right side is computed first, from the old values
print(round(x, 4), round(y, 4))

0.7071 0.5
0.7071 0.7071

The right answer is (0.7071, 0.7071). Python's tuple assignment computes the whole right side before it writes anything, which is exactly what a pair needs.

Common mistakes: the mask 1 << k. Visiting both i0 and i1 as starting points. Overwriting the first amplitude of a pair before computing the second.

Practice: One Gate, Eighteen Qubits, No Giant Matrix

numpy: one axis per qubit

numpy can do the pairing for you. An index splits into three parts: the bits of qubits 0 \dots k-1 (a number x below 2^k), qubit k's bit b, and the bits of the qubits after it (a number y):

reshape(2 ** k, 2, -1) cuts the list exactly this way, because numpy fills an array row by row with the last index running fastest, the same order as the bits. The middle axis is qubit k: slice [:, 0, :] holds the amplitudes where it reads 0, slice [:, 1, :] their partners.

import numpy as np
n, k = 3, 1
psi = np.arange(8)
t = psi.reshape(2 ** k, 2, -1)
print(t[:, 0, :].tolist(), t[:, 1, :].tolist())
print(t[:, ::-1, :].reshape(-1).tolist())                    # X on qubit 1: swap the two slices
z = t.copy(); z[:, 1, :] *= -1
print(z.reshape(-1).tolist())                                # Z on qubit 1: negate slice 1
print(psi.reshape(2, 2, 2).transpose(2, 1, 0).reshape(-1).tolist())   # qubits 0 and 2 exchanged

[[0, 1], [4, 5]] [[2, 3], [6, 7]]
[2, 3, 0, 1, 6, 7, 4, 5]
[0, 1, -2, -3, 4, 5, -6, -7]
[0, 4, 2, 6, 1, 5, 3, 7]

The slices are the pairs of the previous section: [0, 1] with [2, 3], [4, 5] with [6, 7]. Swapping the slices is X, negating slice 1 is Z, and a general gate mixes them: the new slice 0 is g_{00} times slice 0 plus g_{01} times slice 1, and the new slice 1 is g_{10} times slice 0 plus g_{11} times slice 1. -1 in a reshape means "whatever is left", and reshape(-1) flattens back to a list. k = 0 works too: the first axis then has length 1.

With reshape([2] * n) every qubit gets its own axis, qubit 0 first, and transpose renumbers the qubits: the last line exchanges qubits 0 and 2, so 001 (index 1) and 100 (index 4) trade places. That is the SWAP gate.

Common mistakes: reshape(2, 2 ** k, -1), which cuts at the wrong bit. Forgetting .copy() before changing a reshaped array in place: a reshape shares memory with the original. Printing arrays straight away instead of through .tolist().

Practice: Rewiring the Register

CNOT: flip one qubit if another reads 1

One-qubit gates never make qubits depend on each other. The controlled-NOT, CNOT(control c, target t), does: in every basis state where qubit c reads 1 it flips qubit t, and where qubit c reads 0 it does nothing. On two qubits, CNOT(0, 1) maps |00\rangle \to |00\rangle, |01\rangle \to |01\rangle, |10\rangle \to |11\rangle, |11\rangle \to |10\rangle; its matrix swaps the last two entries:

In terms of indices: an index whose control bit is set moves to i ^ target_mask; all others stay. The two qubits need not be neighbours, and the control may be the larger number. Here is where CNOT(1, 2) sends each basis state of 3 qubits:

n, c, t = 3, 1, 2
cmask, tmask = 1 << (n - 1 - c), 1 << (n - 1 - t)
for i in range(2 ** n):
    j = i ^ tmask if i & cmask else i
    print(format(i, "03b"), "->", format(j, "03b"), "(moved)" if i != j else "(stays)")

000 -> 000 (stays)
001 -> 001 (stays)
010 -> 011 (moved)
011 -> 010 (moved)
100 -> 100 (stays)
101 -> 101 (stays)
110 -> 111 (moved)
111 -> 110 (moved)

Only the four states with qubit 1 reading 1 move, and they swap in pairs (010 with 011, 110 with 111). CNOT never mixes amplitudes; it moves them, so it is its own inverse: applying it twice moves everything back. All the care of the pairs section applies: read from the old state, and do not swap a pair twice.

CNOT(1, 0) on two qubits is a different gate from CNOT(0, 1): it swaps |01\rangle and |11\rangle. And CNOT is not a Kronecker product A \otimes B of two one-qubit gates; that is exactly why it can do what they cannot.

Common mistakes: masks 1 << c, which number the qubits from the right. Flipping when the control reads 0. Swapping control and target.

Practice: Flip It Only If the Other Reads 1, The Simulator Runs a Whole Schedule

Bell states, and equal up to a phase

Start from |00\rangle, apply H to qubit 0 and then CNOT(0, 1). H gives (|00\rangle + |10\rangle)/\sqrt2, and CNOT turns the |10\rangle into |11\rangle:

As matrices, the circuit is CNOT @ np.kron(H, I) (the gate applied first stands on the right, Level 2). Its columns are what it makes from |00\rangle, |01\rangle, |10\rangle, |11\rangle, the four Bell states:

import numpy as np
H = np.array([[1, 1], [1, -1]]) / np.sqrt(2)
CNOT = np.array([[1, 0, 0, 0], [0, 1, 0, 0], [0, 0, 0, 1], [0, 0, 1, 0]])
U = CNOT @ np.kron(H, np.eye(2))
for j, name in enumerate(["00", "01", "10", "11"]):
    print("|" + name + "> ->", (np.round(U[:, j], 4) + 0.0).tolist())

|00> -> [0.7071, 0.0, 0.0, 0.7071]
|01> -> [0.0, 0.7071, 0.7071, 0.0]
|10> -> [0.7071, 0.0, 0.0, -0.7071]
|11> -> [0.0, 0.7071, -0.7071, 0.0]

so the inputs |00\rangle, |01\rangle, |10\rangle, |11\rangle give \Phi^+, \Psi^+, \Phi^-, \Psi^-. Measure \Phi^+ and you get 00 or 11, each with probability 1/2: each qubit alone is a fair coin, yet the two always agree. \Psi^\pm always disagree.

When are two register states the same? Level 2's rule carries over: two states that differ by one factor e^{i\phi} on all amplitudes are the same state. For unit vectors there is a one-line test,

since |\langle a|b\rangle| \le \|a\|\,\|b\| = 1 with equality only when one is a multiple of the other. np.vdot(a, b) is \langle a|b\rangle, conjugate on the first vector:

import numpy as np
r = 1 / np.sqrt(2)
phi_minus = np.array([r, 0, 0, -r])
tests = {"e^{0.3i} Phi-": np.exp(0.3j) * phi_minus,
         "Phi+": np.array([r, 0, 0, r]),
         "relative phase": np.array([r, 0, 0, 1j * r])}
for name, b in tests.items():
    f = abs(np.vdot(phi_minus, b)) ** 2
    print(f"{name:15s} {round(float(f), 4):6}  same state: {bool(f > 1 - 1e-9)}")

e^{0.3i} Phi-      1.0  same state: True
Phi+               0.0  same state: False
relative phase     0.5  same state: False

A factor on all four amplitudes passes; a factor on one amplitude (a relative phase) gives a different state, here halfway between.

Aside (optional): what is entangled, physically? In a lab the two qubits of \Phi^+ can be two photons from one crystal, or two ions in one trap, later moved apart. Measured far apart, they agree every time, and each side on its own sees a fair coin. No signal passes between them: nothing one side does changes the odds the other side sees, which is why entanglement alone cannot send a message. What it can do, as the protocols below show, is make a few classical bits go much further.

Common mistakes: comparing amplitudes directly, which fails for a global phase. Testing == 1 instead of a tolerance. Forgetting the conjugate (np.dot instead of np.vdot).

Practice: Which Two Bits Made This Pair?, One Coin Flip, Passed Down the Line

Measuring one qubit of many

Level 2 already found the probability that one qubit of a register reads 1: the total probability of all basis states in which it reads 1, since those outcomes exclude each other:

The new part is what happens next. After the reading the register collapses: it is still a register of n qubits, but every amplitude that disagrees with the outcome becomes 0, and the survivors are divided by \sqrt{p}, the square root of the outcome's probability, so that the weights add up to 1 again. Here is a state of 3 qubits whose qubit 2 reads 1:

import numpy as np
def show(v):
    v = np.round(np.asarray(v, dtype=complex), 4) + 0
    return [x.real if x.imag == 0 else x for x in v.tolist()]

psi = np.array([0.5, 0, 0.5j, 0.5, 0, -0.5, 0, 0])
n, q = 3, 2
mask = 1 << (n - 1 - q)
reads1 = np.array([(i & mask) != 0 for i in range(2 ** n)])
p1 = float(np.sum(np.abs(psi[reads1]) ** 2))
print("qubit 2 reads 1 at", np.flatnonzero(reads1).tolist(), " p0 =", round(1 - p1, 4), " p1 =", round(p1, 4))
after = np.where(reads1, psi, 0) / np.sqrt(p1)
print(show(after), round(float(np.sum(np.abs(after) ** 2)), 6))
print("new odds:", (np.round(np.abs(after) ** 2, 4) + 0.0).tolist())

qubit 2 reads 1 at [1, 3, 5, 7]  p0 = 0.5  p1 = 0.5
[0.0, 0.0, 0.0, 0.7071, 0.0, -0.7071, 0.0, 0.0] 1.0
new odds: [0.0, 0.0, 0.0, 0.5, 0.0, 0.5, 0.0, 0.0]

Before the reading, 011 and 101 each had probability 1/4; after it they have 1/2 each. That is conditional probability from the Statistics path: P(\text{basis state} \mid \text{qubit 2 read } 1) = |a_i|^2 / p_1. Amplitudes are divided by \sqrt{p_1} so that their squares are divided by p_1.

With entanglement this has a striking effect. Measure qubit 0 of \Phi^+ = (|00\rangle + |11\rangle)/\sqrt2 and see 1: only 11 survives, and qubit 1, never touched, now reads 1 for certain. Before the reading it was a fair coin.

Common mistakes: dividing by p instead of \sqrt p. Dropping the measured qubit and returning half the list: the register keeps all n qubits. Dividing by zero when the outcome has probability 0: then there is no state to collapse to.

Practice: Reading One Qubit, Keeping the Rest

Measurements in a row, with a seed

To simulate a reading you must pick the outcome at random, and to make the run reproducible you use Level 2's procedure: rng = random.Random(seed), created once; for each measurement draw r = rng.random() and pick the first outcome whose running total of probabilities is greater than r. For one qubit: outcome 0 if p0 > r, else outcome 1 (and if rounding leaves no such outcome, the last one with nonzero probability). Then collapse, and continue with the collapsed state.

Every measurement draws exactly one number, even when its outcome is certain. Here \Phi^+ is measured, qubit 0 first; after its collapse qubit 1 is certain; a third qubit in |+\rangle is read last. The version that skips the certain draw gives other answers:

import random
def readouts(seed, skip_certain_draw):
    rng = random.Random(seed)
    r = rng.random()
    m0 = 0 if 0.5 > r else 1                 # qubit 0 of Phi+: p0 = 0.5
    p0 = 1.0 if m0 == 0 else 0.0             # collapsed to |00> or |11>: qubit 1 is certain
    if skip_certain_draw:
        m1 = m0
    else:
        r = rng.random()
        m1 = 0 if p0 > r else 1
    r = rng.random()
    m2 = 0 if 0.5 > r else 1                 # a third qubit in |+>
    return f"{m0}{m1}{m2}"

for seed in range(6):
    print(seed, readouts(seed, False), readouts(seed, True))

0 110 111
1 001 001
2 110 111
3 000 001
4 000 000
5 111 111

The first two bits always agree, as \Phi^+ promises. Skipping the certain draw hands the third measurement a different random number, and for seeds 0, 2 and 3 that changes its outcome, so a seeded test tells the two programs apart; the right one draws.

Common mistakes: creating the generator inside the loop, so every measurement gets the same first number. Skipping the draw for a certain outcome. Picking outcome 1 when r < p1: that is a different procedure, with different results.

Practice: Readouts in the Middle of the Schedule

Product or entangled?

Write the four amplitudes of two qubits as a 2 by 2 table, row = qubit 0's bit, column = qubit 1's bit. state.reshape(2, 2) does exactly that, because qubit 0 is the most significant bit. For a product [a, b] \otimes [c, d] the table is

Both rows are multiples of the same row [c, d]: the table has rank 1. For a 2 by 2 table that happens exactly when its determinant a_{00}a_{11} - a_{01}a_{10} is 0, and conversely a determinant of 0 lets you read off [c, d] from a nonzero row. So a two-qubit state is a product state when the determinant is 0, and entangled when it is not. The amplitudes are floats, so "0" means "below a tolerance":

import numpy as np
r = 1 / np.sqrt(2)
states = {
    "kron([0.8, 0.6j], [r, -r])": np.kron([0.8, 0.6j], [r, -r]),
    "Phi+": np.array([r, 0, 0, r]),
    "[0, 0, 0.6, 0.8]": np.array([0, 0, 0.6, 0.8]),
    "[0.5, 0.5, 0.5j, -0.5j]": np.array([0.5, 0.5, 0.5j, -0.5j]),
}
for name, s in states.items():
    T = s.reshape(2, 2)
    det = T[0, 0] * T[1, 1] - T[0, 1] * T[1, 0]
    print(f"{name:27s} |det| = {abs(det):.4f}   below 1e-9: {bool(abs(det) < 1e-9)}")

kron([0.8, 0.6j], [r, -r])  |det| = 0.0000   below 1e-9: True
Phi+                        |det| = 0.5000   below 1e-9: False
[0, 0, 0.6, 0.8]            |det| = 0.0000   below 1e-9: True
[0.5, 0.5, 0.5j, -0.5j]     |det| = 0.5000   below 1e-9: False

Zeros do not decide it: [0, 0, 0.6, 0.8] is |1\rangle \otimes [0.6, 0.8], a product with two zeros, and the last state has no zero at all and is entangled. |\det| runs from 0 for a product to 1/2 for a Bell state, and it ignores global phases, since a factor e^{i\phi} on all four amplitudes multiplies the determinant by e^{2i\phi}.

The Statistics path's independence is the same picture with probabilities: two events are independent exactly when their joint table is the product of the two separate tables, one row times one column. A product state has that kind of table; the probability table of \Phi^+, \begin{pmatrix}1/2 & 0\\ 0 & 1/2\end{pmatrix}, does not.

Common mistakes: testing det == 0. Comparing a complex determinant with a number: for Python's complex numbers det > 1e-9 raises an error, and numpy's silently compare the real parts first. Take abs first. Deciding by "some amplitude is 0".

Practice: Can the Pair Be Pulled Apart?

Schmidt coefficients: how entangled?

The determinant answers yes or no for two qubits. For more qubits, and for "how much", cut the register into part A, qubits 0 \dots m-1, and part B, the rest. Every index is A's bits followed by B's bits, so the table with A's basis states as rows and B's as columns is just state.reshape(2 ** m, 2 ** (n - m)). Its singular values above a small tolerance, largest first, are the Schmidt coefficients s_k of the cut, and how many there are is the Schmidt rank: rank 1 means A and B are a product, a larger rank means they are entangled across the cut. numpy computes them with np.linalg.svd(M, compute_uv=False).

For a normalised state the squares \lambda_k = s_k^2 add up to 1, like probabilities, and the entanglement entropy

is 0 for a product and 1 for a Bell pair. Here are two states of 3 qubits, each cut in its two places:

import math
import numpy as np
r = 1 / math.sqrt(2)
examples = {"|0> and Phi+ on qubits 1, 2": np.kron([1, 0], [r, 0, 0, r]),
            "(|000> + |111>) / sqrt 2": np.array([r, 0, 0, 0, 0, 0, 0, r])}
n = 3
for name, state in examples.items():
    for m in (1, 2):
        s = np.linalg.svd(state.reshape(2 ** m, 2 ** (n - m)), compute_uv=False)
        lam = [float(x) ** 2 for x in s if x > 1e-9]
        E = -sum(l * math.log2(l) for l in lam)
        print(f"{name:28s} cut after qubit {m - 1}: s = {[round(math.sqrt(l), 4) for l in lam]}, E = {round(E, 4) + 0.0}")

|0> and Phi+ on qubits 1, 2  cut after qubit 0: s = [1.0], E = 0.0
|0> and Phi+ on qubits 1, 2  cut after qubit 1: s = [0.7071, 0.7071], E = 1.0
(|000> + |111>) / sqrt 2     cut after qubit 0: s = [0.7071, 0.7071], E = 1.0
(|000> + |111>) / sqrt 2     cut after qubit 1: s = [0.7071, 0.7071], E = 1.0

In the first state qubit 0 stands alone, so the cut after qubit 0 finds a product, and the cut after qubit 1 splits the Bell pair. In the three-qubit state (|000\rangle + |111\rangle)/\sqrt2 every cut finds entanglement.

Between product and Bell pair lies a whole range. For two qubits \cos t\,|00\rangle + \sin t\,|11\rangle:

import math
import numpy as np
ts = [k * math.pi / 40 for k in range(21)]
ents = []
for t in ts:
    s = np.linalg.svd(np.array([[math.cos(t), 0], [0, math.sin(t)]]), compute_uv=False)
    lam = [float(x) ** 2 for x in s if x > 1e-9]
    ents.append(-sum(l * math.log2(l) for l in lam) + 0.0)
print("t = 0:", round(ents[0], 4), " t = pi/5:", round(ents[8], 4), " t = pi/4:", round(ents[10], 4))
s = np.linalg.svd(np.array([[0.8, 0], [0, 0.6]]), compute_uv=False)
print(s.tolist(), round(float(s[0] * s[1]), 4))
plot(ts, ents, title="entanglement entropy of cos t |00> + sin t |11>", xlabel="t", ylabel="E")

t = 0: 0.0  t = pi/5: 0.93  t = pi/4: 1.0
[0.8, 0.6] 0.48

Press Run to see the entropy climb from 0 at t = 0 to 1 at t = \pi/4, the Bell state, and fall back. For two qubits the product of the two coefficients is the |\det| of the previous section: 0.8 \cdot 0.6 = 0.48 for [0.8, 0, 0, 0.6].

Common mistakes: reshape(2 ** (n - m), 2 ** m), which cuts in another place. Keeping singular values like 1e-17, which adds fake coefficients and makes math.log2(0) fail. Forgetting to square the coefficients before the entropy.

Practice: Cutting the Register in Two

Expectation values and correlations

Level 2 averaged one qubit's reading. With a matrix M that describes a measurement (an observable), the average result is the expectation value

which numpy writes as np.vdot(psi, M @ psi), the conjugate on the left. Z \otimes Z is diagonal with +1 where the two bits agree and -1 where they differ, so \langle Z \otimes Z\rangle is p_{\text{same}} - p_{\text{diff}}: +1 if the two qubits always agree, -1 if they never do, 0 if agreeing is a coin toss.

import numpy as np
r = 1 / np.sqrt(2)
Z = np.diag([1, -1])
X = np.array([[0, 1], [1, 0]])
ZZ = np.kron(Z, Z)
print(np.diag(ZZ).tolist())
for name, psi in (("Phi+", [r, 0, 0, r]), ("Psi+", [0, r, r, 0]), ("|0>|+>", [r, r, 0, 0])):
    psi = np.array(psi)
    print(f"{name:7s} <ZZ> = {round(float(np.vdot(psi, ZZ @ psi).real), 4) + 0.0}")
psi = np.array([r, 0, 0, 1j * r])
XX = np.kron(X, X)
print("with the conjugate:", round(float(np.vdot(psi, XX @ psi).real), 4) + 0.0, "  without:", complex(np.round(psi @ XX @ psi, 4)))

[1, -1, -1, 1]
Phi+    <ZZ> = 1.0
Psi+    <ZZ> = -1.0
|0>|+>  <ZZ> = 0.0
with the conjugate: 0.0   without: 1j

Without the conjugate, (|00\rangle + i|11\rangle)/\sqrt2 gives the "average" i, which no list of \pm1 readings can have.

Measuring along an angle \theta instead of 0/1 uses the observable A(\theta) = \cos\theta\, Z + \sin\theta\, X = \begin{pmatrix}\cos\theta & \sin\theta\\ \sin\theta & -\cos\theta\end{pmatrix}, with A(0) = Z and A(\pi/2) = X. If one lab measures qubit 0 along a and another qubit 1 along b, their correlation is E(a, b) = \langle A(a) \otimes A(b)\rangle. For \Phi^+ it works out to \cos(a - b):

import math
import numpy as np
r = 1 / math.sqrt(2)
phi = np.array([r, 0, 0, r])
def A(t):
    return np.array([[math.cos(t), math.sin(t)], [math.sin(t), -math.cos(t)]])
def E(a, b):
    return float(np.vdot(phi, np.kron(A(a), A(b)) @ phi).real)
print(round(E(0.3, 1.0), 4), round(math.cos(0.3 - 1.0), 4))
bs = [k * math.pi / 20 for k in range(-20, 21)]
plot(bs, [E(0, b) for b in bs], title="E(0, b) for Phi+", xlabel="b", ylabel="E")

0.7648 0.7648

Press Run to see the correlation fall from +1 at b = 0 to -1 at b = \pm\pi. The CHSH test combines four such correlations, two settings per lab:

Any product state, and any pair of labs following answers agreed in advance, keeps |S| \le 2. Entangled pairs reach 2\sqrt2 \approx 2.83 at well-chosen angles. Real experiments measure S above 2, so no plan agreed in advance can explain what entangled pairs do.

Common mistakes: leaving out the conjugate. Putting qubit 1's matrix on the left of the Kronecker product. Using degrees.

Practice: Do the Two Counters Click Together?, Four Dial Settings and a Score Above Two

Teleportation

A sender wants to give a receiver the state |\psi\rangle = [\alpha, \beta] of a qubit, without knowing \alpha and \beta and with only a phone line for classical bits. They share a Bell pair \Phi^+. The three qubits:

• qubit 0 holds |\psi\rangle (the sender's);
• qubit 1 is the sender's half of the pair, qubit 2 the receiver's half.

The start state is kron(psi, pair). The protocol:

1. The sender applies CNOT(0, 1), then H to qubit 0.
2. The sender measures qubits 0 and 1, getting bits m_0 m_1, and phones them over.
3. The receiver corrects qubit 2: X if m_1 = 1, then Z if m_0 = 1.

After step 2, outcome m_0 m_1 leaves qubit 2 holding the two amplitudes at |m_0 m_1 0\rangle and |m_0 m_1 1\rangle (indices 4m_0 + 2m_1 and one more), divided by their length. With full matrices on 8 amplitudes:

import numpy as np
def show(v):
    v = np.round(np.asarray(v, dtype=complex), 4) + 0
    return [x.real if x.imag == 0 else x for x in v.tolist()]

r = 1 / np.sqrt(2)
X = np.array([[0, 1], [1, 0]])
Z = np.diag([1, -1])
H = np.array([[1, 1], [1, -1]]) * r
CNOT = np.array([[1, 0, 0, 0], [0, 1, 0, 0], [0, 0, 0, 1], [0, 0, 1, 0]])

psi = np.array([0.8, 0.6j])
state = np.kron(psi, [r, 0, 0, r])              # qubit 0: psi; qubits 1 and 2: Phi+
state = np.kron(CNOT, np.eye(2)) @ state        # CNOT(0, 1)
state = np.kron(H, np.eye(4)) @ state           # H on qubit 0
for m0 in (0, 1):
    for m1 in (0, 1):
        i = 4 * m0 + 2 * m1
        u = state[i:i + 2]
        p = float(np.sum(np.abs(u) ** 2))
        q = u / np.sqrt(p)
        fixed = q
        if m1:
            fixed = X @ fixed
        if m0:
            fixed = Z @ fixed
        print(f"{m0}{m1}: p = {p:.2f}  qubit 2 {show(q)}  corrected {show(fixed)}")

00: p = 0.25  qubit 2 [0.8, 0.6j]  corrected [0.8, 0.6j]
01: p = 0.25  qubit 2 [0.6j, 0.8]  corrected [0.8, 0.6j]
10: p = 0.25  qubit 2 [0.8, -0.6j]  corrected [0.8, 0.6j]
11: p = 0.25  qubit 2 [-0.6j, 0.8]  corrected [0.8, 0.6j]

Each outcome has probability 1/4 whatever |\psi\rangle is, so the two bits on the phone line say nothing about \alpha and \beta. Before the correction the receiver holds |\psi\rangle with an X, a Z, or both applied; after it, |\psi\rangle itself.

Corrections only matter up to a global phase, so the same card can be written in more than one way. "X, then Z" is the matrix ZX, and ZX = iY, while XZ = -iY: both orders are Y up to a factor of modulus 1, and the receiver may apply any of the three.

import numpy as np
X = np.array([[0, 1], [1, 0]])
Y = np.array([[0, -1j], [1j, 0]])
Z = np.diag([1, -1])
print((Z @ X).tolist(), bool(np.allclose(Z @ X, 1j * Y)), bool(np.allclose(X @ Z, -1j * Y)))

[[0, 1], [-1, 0]] True True

A correction must work for every |\psi\rangle. Checking it on |\psi\rangle = |0\rangle alone proves little: Z|0\rangle = |0\rangle, so there Z and I look the same. Test with a state that has two nonzero amplitudes, a complex one included, as above.

Aside (optional): no copying. Teleportation moves a state, it does not copy it: after the measurement qubit 0 holds |m_0\rangle and |\psi\rangle is gone from the sender's side. That is not a weakness of this protocol. No circuit can copy an unknown qubit (the no-cloning theorem), and the maths says why: a gate that sent |0\rangle|0\rangle \to |0\rangle|0\rangle and |1\rangle|0\rangle \to |1\rangle|1\rangle, like CNOT, is linear, so it sends |+\rangle|0\rangle to \Phi^+, not to |+\rangle|+\rangle. The same fact makes eavesdropping detectable in the key exchange below.

Common mistakes: reading the receiver's qubit as the leftmost bit. Forgetting to divide by the length. Applying the corrections to the wrong outcome bits (X belongs to m_1, Z to m_0 here). Testing a correction on |0\rangle only.

Practice: The Card Taped to the Receiver, Sending a Qubit Down the Corridor

Superdense coding

Teleportation sends one qubit with two classical bits and a shared pair. Superdense coding runs the other way: two classical bits with one qubit and a shared pair. The sender holds qubit 0 of \Phi^+ and encodes two bits by one gate on it: I for 00, X for 01, Z for 10, X then Z for 11. Each turns \Phi^+ into a different Bell state:

The sender posts qubit 0. The receiver, now holding both, runs the Bell circuit backwards: CNOT(0, 1), then H on qubit 0. That undoes the circuit that made the Bell states from |00\rangle, |01\rangle, |10\rangle, |11\rangle, so each Bell state becomes its basis state and a measurement reads the two bits for certain:

import numpy as np
r = 1 / np.sqrt(2)
I = np.eye(2)
X = np.array([[0, 1], [1, 0]])
Z = np.diag([1, -1])
H = np.array([[1, 1], [1, -1]]) * r
CNOT = np.array([[1, 0, 0, 0], [0, 1, 0, 0], [0, 0, 0, 1], [0, 0, 1, 0]])
make = CNOT @ np.kron(H, I)                    # H on qubit 0, then CNOT(0, 1)
undo = np.kron(H, I) @ CNOT                    # CNOT(0, 1), then H on qubit 0
print(bool(np.allclose(undo @ make, np.eye(4))))
bell = {"Phi+": make[:, 0], "Psi+": make[:, 1], "Phi-": make[:, 2], "Psi-": make[:, 3]}
for name, G in (("I", I), ("X", X), ("Z", Z), ("ZX", Z @ X)):
    sent = np.kron(G, I) @ bell["Phi+"]
    match = [b for b, v in bell.items() if abs(np.vdot(v, sent)) ** 2 > 1 - 1e-9]
    print(name, "turns Phi+ into", match)

True
I turns Phi+ into ['Phi+']
X turns Phi+ into ['Psi+']
Z turns Phi+ into ['Phi-']
ZX turns Phi+ into ['Psi-']

Both circuits are made of gates that are their own inverses, so the undoing circuit is the same gates in reverse order. A gate outside the code book, H or T for example, turns \Phi^+ into a state that is not one Bell state, and the receiver's reading becomes uncertain.

Common mistakes: decoding with H before the CNOT. Applying the sender's gate to qubit 1. Expecting XZ and ZX to send different messages: they differ only by the factor -1.

Practice: Two Bits Folded into One Qubit

BB84: a key from single qubits

Two stations, A and B, want a shared secret key of random bits. BB84 uses no entanglement, only single qubits and two bases:

• A picks a random bit b and a random basis. In basis "+" A sends |b\rangle; in basis "x" A sends H|b\rangle, which is |+\rangle or |-\rangle.
• B picks a random basis too and measures in it: for "+" an ordinary 0/1 reading, for "x" first H, then the reading.

When the bases agree, B reads A's bit for certain (H undoes H). When they differ, B gets a fair coin.

import numpy as np
H = np.array([[1, 1], [1, -1]]) / np.sqrt(2)
def prepare(bit, basis):
    q = np.eye(2)[bit]
    return H @ q if basis == "x" else q
def odds(q, basis):                        # Born rule for reading 0 and 1 in this basis
    if basis == "x":
        q = H @ q
    return np.abs(q) ** 2

for a in "+x":
    print("A sends 0 in", a, "-> B reads 0 with probability",
          {b: round(float(odds(prepare(0, a), b)[0]), 4) for b in "+x"})

err = 0.0
for e in "+x":                             # a spy measures in her own basis, each with probability 1/2,
    for e_bit, pe in enumerate(odds(prepare(0, "+"), e)):
        resent = prepare(e_bit, e)         # and sends a fresh qubit prepared from her result
        err += 0.5 * pe * odds(resent, "+")[1]   # B measures in A's basis "+" and reads 1
print("errors in sifted bits with a spy on every qubit:", round(float(err), 4))

A sends 0 in + -> B reads 0 with probability {'+': 1.0, 'x': 0.5}
A sends 0 in x -> B reads 0 with probability {'+': 0.5, 'x': 1.0}
errors in sifted bits with a spy on every qubit: 0.25

Afterwards A and B announce their bases in public, never their bits, and keep only the positions where the bases agreed. That is sifting: about half the positions survive, and there the two keys should be equal.

A spy who intercepts a qubit, measures it and resends what she saw must guess the basis. With the right guess she learns the bit and leaves no trace; with the wrong one (half the time) her resent qubit makes B's sifted bit a coin toss. That is the 0.25 above: an error in a quarter of the intercepted sifted bits, so a spy on a fraction f of the qubits causes an error rate (the QBER) of about f/4. To catch her, A and B publish a random sample of their sifted bits, count mismatches, and abort if the error rate is too high; the published bits are then dropped from the key. Without a spy (and without noise) a sifted bit never disagrees.

Common mistakes: sifting on equal bits instead of equal bases. Measuring in "x" without applying H first. Keeping the sample bits that were announced in public.

Practice: Keeping Only the Columns That Match, Burn the Key or Keep It?, Someone Is Tapping the Fibre

Practice

These problems are required on this path at this level; the path adds one problem from each of the level's topics, picked for you and kept.

• Qubits Prepared One by One
• One Gate, Eighteen Qubits, No Giant Matrix
• Flip It Only If the Other Reads 1
• Reading One Qubit, Keeping the Rest
• Can the Pair Be Pulled Apart?
• The Card Taped to the Receiver

Checklist

You can now:

• Convert between a basis index and its bit string with qubit 0 as the leftmost, most significant bit, and read, flip, set and clear one qubit's bit with the mask 1 << (n - 1 - q).
• Build a product state with np.kron, qubit 0 first, and the full matrix I \otimes \dots \otimes G \otimes \dots \otimes I of a gate on one qubit, and say why that matrix is out of reach at 18 qubits.
• Apply a one-qubit gate on pairs of indices (each pair once, both old values read first), or with reshape(2 ** k, 2, -1) in numpy.
• Apply CNOT between any two qubits by moving amplitudes with bit masks, and make the four Bell states.
• Test whether two states are equal up to a global phase with |\langle a|b\rangle|^2 = 1 and a tolerance.
• Measure one qubit of a register: its probability, the collapsed state of all n qubits, and seeded outcomes with one draw per measurement.
• Tell a two-qubit product state from an entangled one by the determinant of its 2 by 2 table, and measure entanglement across any cut with singular values and the entropy.
• Compute expectation values \langle\psi|M|\psi\rangle with the conjugate, correlations of two qubits, and say what the CHSH value shows.
• Explain teleportation, superdense coding and BB84 step by step, including why corrections only matter up to a global phase and why a spy raises the error rate.