Bloch Sphere Mathematics: From State Vectors to Rotations
The mathematics behind the Bloch sphere: how the angles θ and φ map to quantum state vectors, why single-qubit gates are rotations, and how to read any single-qubit state geometrically.
The Bloch sphere is not just a diagram. It is an exact mathematical representation of every possible single-qubit state. This tutorial derives the parametrisation from first principles, explains why gates are rotations, connects the sphere to the group theory of SU(2) and SO(3), and builds up to noise models, tomography, and gate approximation theory.
Why a sphere?
A general single-qubit state is a unit vector in a two-dimensional complex Hilbert space:
|ψ⟩ = α|0⟩ + β|1⟩
where α and β are complex numbers satisfying |α|² + |β|² = 1. Complex numbers have two real parameters each, so naively |ψ⟩ seems to have four real degrees of freedom. Two constraints reduce this:
- Normalisation: |α|² + |β|² = 1 removes one degree of freedom.
- Global phase: |ψ⟩ and e^(iφ)|ψ⟩ are physically identical, because no measurement can distinguish them. This removes one more.
Four minus two equals two. A two-parameter family of unit vectors in ℝ³ is exactly a sphere. That is why a sphere works.
The parametrisation
Write α and β in polar form and absorb the global phase into α (making α real and non-negative):
α = cos(θ/2)
β = e^(iφ) sin(θ/2)
where θ ∈ [0, π] and φ ∈ [0, 2π). The full state is:
|ψ⟩ = cos(θ/2)|0⟩ + e^(iφ) sin(θ/2)|1⟩
The corresponding point on the Bloch sphere has Cartesian coordinates:
x = sin(θ) cos(φ)
y = sin(θ) sin(φ)
z = cos(θ)
This is just the standard spherical-to-Cartesian conversion with radius 1.
Why θ/2 and not θ?
This is the subtlety that trips people up. If the amplitude is cos(θ/2), then at θ=0 we get α=1, β=0, which is the state |0⟩ at the north pole. At θ=π we get α=0, β=1, which is the state |1⟩ at the south pole. The factor of 1/2 means a full rotation of the Bloch vector (360°) corresponds to a 720° rotation of the state vector. This is the defining property of spin-1/2 systems and is why fermions require two full turns to return to their original state.
Important points on the sphere
| State | θ | φ | Bloch vector |
|---|---|---|---|
| |0⟩ | 0 | any | (0, 0, 1) north pole |
| |1⟩ | π | any | (0, 0, −1) south pole |
| |+⟩ = (|0⟩+|1⟩)/√2 | π/2 | 0 | (1, 0, 0) +x axis |
| |−⟩ = (|0⟩−|1⟩)/√2 | π/2 | π | (−1, 0, 0) −x axis |
| |i⟩ = (|0⟩+i|1⟩)/√2 | π/2 | π/2 | (0, 1, 0) +y axis |
| |−i⟩ = (|0⟩−i|1⟩)/√2 | π/2 | 3π/2 | (0, −1, 0) −y axis |
You can verify each using the parametrisation above.
Reading measurements from the Bloch vector
The expectation values of the Pauli observables X, Y, Z equal the Bloch vector components directly:
⟨X⟩ = ⟨ψ|X|ψ⟩ = x = sin(θ) cos(φ)
⟨Y⟩ = ⟨ψ|Y|ψ⟩ = y = sin(θ) sin(φ)
⟨Z⟩ = ⟨ψ|Z|ψ⟩ = z = cos(θ)
The probability of measuring |0⟩ (outcome +1 for Z) is:
P(0) = |cos(θ/2)|² = (1 + cos(θ))/2 = (1 + z)/2
So if you know the z-coordinate of the Bloch vector, you know the measurement probability immediately.
Gates as rotations
Every single-qubit gate corresponds to a rotation of the Bloch vector. The general rotation by angle α around axis n = (nx, ny, nz) is:
R_n(α) = exp(−i α/2 (nx·X + ny·Y + nz·Z))
= cos(α/2)·I − i sin(α/2)·(nx·X + ny·Y + nz·Z)
This is the rotation-operator formula from the theory of angular momentum (SU(2)).
The standard single-qubit gates are:
Pauli gates: 180° rotations
X gate: rotation by π around x-axis
X|0⟩ = |1⟩ (north pole → south pole)
X|+⟩ = |+⟩ (fixed: +x axis)
Z gate: rotation by π around z-axis
Z|0⟩ = |0⟩ (fixed: north pole)
Z|+⟩ = |−⟩ (flips φ by π)
Y gate: rotation by π around y-axis
Y|0⟩ = i|1⟩ (north → south with phase)
Phase gates: rotations around z-axis
S gate: rotation by π/2 around z-axis (φ → φ + π/2)
T gate: rotation by π/4 around z-axis (φ → φ + π/4)
Rz(θ) gate: rotation by θ around z-axis
General rotation gates
from qiskit.circuit import QuantumCircuit
import numpy as np
qc = QuantumCircuit(1)
# Rx(π/2): rotate 90° around x-axis
qc.rx(np.pi / 2, 0)
# This takes |0⟩ (north pole) to (0, -1, 0) = |-i⟩
print(qc.draw())
The Hadamard gate as a specific rotation
The Hadamard gate deserves a detailed treatment because it appears in nearly every quantum algorithm and its rotation axis is less obvious than the Pauli gates.
The axis is n = (1, 0, 1)/√2, which points midway between the x and z directions. The rotation angle is π. Plugging into the rotation formula:
R_n(π) = cos(π/2)·I − i sin(π/2)·(n_x·X + n_z·Z)
= 0·I − i·1·(X/√2 + Z/√2)
= −i(X + Z)/√2
= −i·H
So R_n(π) = −i·H, which means H = i·R_{(x+z)/√2}(π). The factor of i is a global phase, so H and R_n(π) produce the same physical state. Let us verify this numerically:
import numpy as np
# Pauli matrices
I = np.eye(2, dtype=complex)
X = np.array([[0, 1], [1, 0]], dtype=complex)
Z = np.array([[1, 0], [0, -1]], dtype=complex)
# Rotation by π around (1, 0, 1)/√2
nx, nz = 1 / np.sqrt(2), 1 / np.sqrt(2)
alpha = np.pi
R = np.cos(alpha / 2) * I - 1j * np.sin(alpha / 2) * (nx * X + nz * Z)
# Standard Hadamard matrix
H = np.array([[1, 1], [1, -1]], dtype=complex) / np.sqrt(2)
print("R_n(π):")
print(np.round(R, 6))
print("\n-i * H:")
print(np.round(-1j * H, 6))
print("\nAre they equal?", np.allclose(R, -1j * H))
# Output: True
Geometrically, the rotation around (x+z)/√2 swaps the x and z axes. This is why:
- H|0⟩ = |+⟩ takes the north pole (z-axis) to the positive x-axis.
- H|+⟩ = |0⟩ takes the positive x-axis back to the north pole.
- H|1⟩ = |−⟩ takes the south pole to the negative x-axis.
The Hadamard is its own inverse (H² = I), which makes sense because two 180° rotations around the same axis return to the starting point.
The SU(2) to SO(3) double cover
The Bloch sphere connects two fundamental groups in physics: SU(2), the group of 2×2 unitary matrices with determinant 1, and SO(3), the group of 3D rotations. There is a 2-to-1 homomorphism from SU(2) to SO(3). This means that for every rotation in 3D, there are exactly two SU(2) matrices that produce it.
Why 2-to-1? Because |ψ⟩ and −|ψ⟩ differ only by a global phase of −1, so they map to the same Bloch vector. A 2π rotation of the Bloch vector (a full 360° turn) maps to a rotation by π in the state vector space (picking up a factor of −1). You need a 4π rotation (720°) to return the state vector to its original value.
import numpy as np
from qiskit import QuantumCircuit
from qiskit.quantum_info import Statevector
# Start with |0⟩
qc = QuantumCircuit(1)
sv_initial = Statevector(qc)
print("Initial state vector:", sv_initial.data)
# [1.+0.j, 0.+0.j]
# Apply Rz(2π): a full 360° rotation around z
qc.rz(2 * np.pi, 0)
sv_after_2pi = Statevector(qc)
print("After Rz(2π):", np.round(sv_after_2pi.data, 6))
# [-1.+0.j, -0.+0.j] (state vector is -|0⟩, NOT |0⟩)
# Compute the Bloch vectors for both
def bloch_from_sv(sv):
rho = np.outer(sv.data, sv.data.conj())
X = np.array([[0, 1], [1, 0]], dtype=complex)
Y = np.array([[0, -1j], [1j, 0]], dtype=complex)
Z = np.array([[1, 0], [0, -1]], dtype=complex)
x = float(np.real(np.trace(rho @ X)))
y = float(np.real(np.trace(rho @ Y)))
z = float(np.real(np.trace(rho @ Z)))
return (round(x, 6), round(y, 6), round(z, 6))
print("Bloch vector before:", bloch_from_sv(sv_initial))
# (0, 0, 1)
print("Bloch vector after Rz(2π):", bloch_from_sv(sv_after_2pi))
# (0, 0, 1) <- identical!
# Now apply another Rz(2π) for a total of 4π (720°)
qc.rz(2 * np.pi, 0)
sv_after_4pi = Statevector(qc)
print("After Rz(4π):", np.round(sv_after_4pi.data, 6))
# [1.+0.j, 0.+0.j] <- back to the original state vector
The state vector after Rz(2π) is −|0⟩, not |0⟩. But the Bloch vector is still (0, 0, 1). This is the double cover in action: two state vectors (|0⟩ and −|0⟩) map to the same physical state and the same point on the Bloch sphere.
This is not just mathematical curiosity. The −1 phase becomes detectable through interference. In neutron interferometry experiments, rotating a neutron by 360° introduces a measurable phase shift of π, confirming this SU(2) structure experimentally.
Euler angle decomposition
Any single-qubit unitary can be decomposed into three rotations: Rz(α) · Ry(β) · Rz(γ). This is the ZYZ Euler decomposition, and it is how hardware compilers break down abstract gates into the native gate set of a quantum processor.
For a general 2×2 unitary U with determinant 1:
U = [[a, -b*],
[b, a*]]
where |a|² + |b|² = 1, the Euler angles are:
α = arg(a) - arg(b) (first z-rotation)
β = 2 · arctan(|b|/|a|) (y-rotation)
γ = arg(a) + arg(b) (second z-rotation)
When b = 0, the unitary is already a z-rotation and β = 0. When a = 0, the y-rotation is π and the z-rotations handle the remaining phase.
Here is code that extracts the Euler angles from any 2×2 unitary and reconstructs it:
import numpy as np
def euler_angles_zyz(U):
"""Extract ZYZ Euler angles from a 2x2 unitary matrix.
Returns (alpha, beta, gamma, global_phase) such that:
U = exp(i * global_phase) * Rz(alpha) @ Ry(beta) @ Rz(gamma)
"""
# Remove global phase to get determinant 1
det = np.linalg.det(U)
phase = np.angle(det) / 2
V = U * np.exp(-1j * phase) # V now has det = 1
a = V[0, 0]
b = V[1, 0]
beta = 2 * np.arctan2(abs(b), abs(a))
if abs(a) < 1e-10:
# a ≈ 0: β = π, only (α - γ) is determined
alpha = -np.angle(b)
gamma = 0.0
elif abs(b) < 1e-10:
# b ≈ 0: β = 0, only (α + γ) is determined
alpha = np.angle(a)
gamma = 0.0
else:
alpha = np.angle(a) - np.angle(b)
gamma = np.angle(a) + np.angle(b)
return alpha, beta, gamma, phase
def rz(theta):
"""Rz(theta) rotation matrix."""
return np.array([
[np.exp(-1j * theta / 2), 0],
[0, np.exp(1j * theta / 2)]
])
def ry(theta):
"""Ry(theta) rotation matrix."""
return np.array([
[np.cos(theta / 2), -np.sin(theta / 2)],
[np.sin(theta / 2), np.cos(theta / 2)]
])
# Test with the Hadamard gate
H = np.array([[1, 1], [1, -1]]) / np.sqrt(2)
alpha, beta, gamma, phase = euler_angles_zyz(H)
print(f"Hadamard Euler angles:")
print(f" alpha = {alpha:.4f} ({alpha/np.pi:.4f}π)")
print(f" beta = {beta:.4f} ({beta/np.pi:.4f}π)")
print(f" gamma = {gamma:.4f} ({gamma/np.pi:.4f}π)")
print(f" global phase = {phase:.4f}")
# Reconstruct and verify
H_reconstructed = np.exp(1j * phase) * rz(alpha) @ ry(beta) @ rz(gamma)
print(f"\nOriginal H:\n{np.round(H, 6)}")
print(f"Reconstructed:\n{np.round(H_reconstructed, 6)}")
print(f"Match: {np.allclose(H, H_reconstructed)}")
# For H: α = π, β = π/2, γ = 0 (up to global phase)
# So H = e^(iφ) · Rz(π) · Ry(π/2) · Rz(0)
The Euler decomposition is why quantum hardware only needs two rotation axes to implement any single-qubit gate. Superconducting processors typically provide Rz (virtual, zero-cost) and a physical rotation like Rx or Ry. With two axes, any point on the Bloch sphere is reachable.
Berry phase and geometric phase
When a quantum state undergoes cyclic evolution (its Bloch vector traces a closed loop and returns to its starting point), the state acquires a phase that depends only on the geometry of the path. This is the Berry phase, and it equals half the solid angle enclosed by the path on the Bloch sphere.
For a path enclosing solid angle Ω:
Berry phase = Ω / 2
Consider a concrete example. Start with the state |+⟩, which sits at the equator at the point (1, 0, 0). Apply Rz(φ) for φ going from 0 to 2π. The Bloch vector traces a great circle on the equator, returning to its starting point.
The solid angle enclosed by the equator is 2π(1 − cos(θ₀)) where θ₀ = π/2 is the polar angle of the equator. So Ω = 2π(1 − cos(π/2)) = 2π(1 − 0) = 2π. The Berry phase is Ω/2 = π.
Let us verify this numerically:
import numpy as np
# |+⟩ state: (|0⟩ + |1⟩)/√2
psi_plus = np.array([1, 1], dtype=complex) / np.sqrt(2)
# Apply Rz(2π) to trace a full equatorial loop
def rz(theta):
return np.array([
[np.exp(-1j * theta / 2), 0],
[0, np.exp(1j * theta / 2)]
])
psi_final = rz(2 * np.pi) @ psi_plus
print("Initial |+⟩:", psi_plus)
print("After Rz(2π):", np.round(psi_final, 6))
# [-0.707107-0.j, -0.707107+0.j]
# This is -|+⟩, a phase of π (= e^(iπ) = -1)
# Extract the phase: psi_final = e^(i*phase) * psi_plus
# Compute overlap
overlap = np.vdot(psi_plus, psi_final)
phase = np.angle(overlap)
print(f"Phase acquired: {phase:.4f} rad = {phase/np.pi:.4f}π")
# Phase: -π (or equivalently π)
# This matches Berry phase = Ω/2 = 2π/2 = π
print(f"Expected Berry phase: π = {np.pi:.4f} rad")
print(f"Match: {np.isclose(abs(phase), np.pi)}")
The state |+⟩ acquires a phase of π after a full Rz(2π) rotation. Note that this phase has two contributions: a dynamical phase and a geometric (Berry) phase. For Rz rotations applied to equatorial states, the entire phase is geometric.
The Berry phase has practical applications. Holonomic quantum computing uses geometric phases as the basis for gate operations, providing natural robustness against certain types of noise because the phase depends only on the area enclosed, not on the speed or details of the path.
Computing the Bloch vector in Qiskit
from qiskit import QuantumCircuit
from qiskit.quantum_info import Statevector
import numpy as np
def bloch_vector(qc: QuantumCircuit) -> tuple:
"""Return (x, y, z) Bloch vector for a single-qubit circuit."""
sv = Statevector(qc)
rho = sv.to_operator().data # density matrix as 2x2 array
# Pauli matrices
X = np.array([[0, 1], [1, 0]], dtype=complex)
Y = np.array([[0, -1j], [1j, 0]], dtype=complex)
Z = np.array([[1, 0], [0, -1]], dtype=complex)
x = float(np.real(np.trace(rho @ X)))
y = float(np.real(np.trace(rho @ Y)))
z = float(np.real(np.trace(rho @ Z)))
return (x, y, z)
# |0⟩ state
qc0 = QuantumCircuit(1)
print(" |0⟩:", bloch_vector(qc0)) # (0, 0, 1)
# |+⟩ state
qcH = QuantumCircuit(1)
qcH.h(0)
print(" |+⟩:", bloch_vector(qcH)) # (1, 0, 0)
# After T gate from |+⟩: rotated π/4 around z
qcT = QuantumCircuit(1)
qcT.h(0)
qcT.t(0)
x, y, z = bloch_vector(qcT)
print("T|+⟩:", (round(x,4), round(y,4), round(z,4)))
# (cos(π/4), sin(π/4), 0) = (0.7071, 0.7071, 0)
Visualising gate sequences on the Bloch sphere
Tracing how the Bloch vector moves during a gate sequence builds strong geometric intuition. The following code applies a sequence of gates and plots the resulting trajectory on the sphere surface.
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
# Pauli matrices
I2 = np.eye(2, dtype=complex)
sx = np.array([[0, 1], [1, 0]], dtype=complex)
sy = np.array([[0, -1j], [1j, 0]], dtype=complex)
sz = np.array([[1, 0], [0, -1]], dtype=complex)
def state_to_bloch(psi):
"""Convert a state vector to a Bloch vector (x, y, z)."""
rho = np.outer(psi, psi.conj())
x = float(np.real(np.trace(rho @ sx)))
y = float(np.real(np.trace(rho @ sy)))
z = float(np.real(np.trace(rho @ sz)))
return np.array([x, y, z])
def rx(theta):
return np.cos(theta / 2) * I2 - 1j * np.sin(theta / 2) * sx
def ry(theta):
return np.cos(theta / 2) * I2 - 1j * np.sin(theta / 2) * sy
def rz(theta):
return np.cos(theta / 2) * I2 - 1j * np.sin(theta / 2) * sz
def trajectory(initial_state, gates):
"""Compute the Bloch vector after each gate in the sequence."""
psi = initial_state.copy()
points = [state_to_bloch(psi)]
for gate in gates:
psi = gate @ psi
points.append(state_to_bloch(psi))
return np.array(points)
# Initial state: |+⟩
plus = np.array([1, 1], dtype=complex) / np.sqrt(2)
# Three gate sequences
sequences = {
"Rz(π/4) × 8": [rz(np.pi / 4)] * 8,
"Rx(π/4) × 8": [rx(np.pi / 4)] * 8,
"Alternating Ry, Rz": [
ry(np.pi / 4) if i % 2 == 0 else rz(np.pi / 4)
for i in range(16)
],
}
fig = plt.figure(figsize=(18, 5))
for idx, (label, gates) in enumerate(sequences.items()):
ax = fig.add_subplot(1, 3, idx + 1, projection="3d")
# Draw wireframe sphere
u = np.linspace(0, 2 * np.pi, 30)
v = np.linspace(0, np.pi, 20)
xs = np.outer(np.cos(u), np.sin(v))
ys = np.outer(np.sin(u), np.sin(v))
zs = np.outer(np.ones_like(u), np.cos(v))
ax.plot_wireframe(xs, ys, zs, alpha=0.1, color="gray")
# Draw axes
for axis, color in [([1, 0, 0], "r"), ([0, 1, 0], "g"), ([0, 0, 1], "b")]:
ax.plot([0, axis[0]], [0, axis[1]], [0, axis[2]],
color=color, linewidth=1, alpha=0.5)
# Plot trajectory
pts = trajectory(plus, gates)
ax.plot(pts[:, 0], pts[:, 1], pts[:, 2], "o-", markersize=4, linewidth=2)
ax.scatter(*pts[0], color="green", s=80, zorder=5, label="start")
ax.scatter(*pts[-1], color="red", s=80, zorder=5, label="end")
ax.set_title(label)
ax.set_xlim([-1.2, 1.2])
ax.set_ylim([-1.2, 1.2])
ax.set_zlim([-1.2, 1.2])
ax.legend(fontsize=8)
plt.tight_layout()
plt.savefig("bloch_trajectories.png", dpi=150)
plt.show()
The first sequence (Rz rotations) traces a circle on the equator, since Rz only changes the azimuthal angle φ. The second sequence (Rx rotations) traces a great circle in the y-z plane, since Rx rotates around the x-axis. The third sequence (alternating Ry and Rz) produces a winding path that covers more of the sphere surface, illustrating how combining rotations around different axes can reach any point.
State tomography from the Bloch vector
If you can measure the three Pauli expectation values ⟨X⟩, ⟨Y⟩, and ⟨Z⟩, you can reconstruct the full density matrix of any single-qubit state. This is quantum state tomography for one qubit.
The procedure:
- Measure in the Z basis (computational basis) to estimate ⟨Z⟩.
- Apply H, then measure in the Z basis to estimate ⟨X⟩ (because HZH = X).
- Apply S†H, then measure in the Z basis to estimate ⟨Y⟩ (because HS · Z · S†H = Y).
From these three values, the density matrix is:
ρ = (I + ⟨X⟩·X + ⟨Y⟩·Y + ⟨Z⟩·Z) / 2
Here is code that simulates this process with realistic shot noise:
import numpy as np
# Pauli matrices
I2 = np.eye(2, dtype=complex)
X = np.array([[0, 1], [1, 0]], dtype=complex)
Y = np.array([[0, -1j], [1j, 0]], dtype=complex)
Z = np.array([[1, 0], [0, -1]], dtype=complex)
def simulate_measurement(psi, num_shots=1000):
"""Simulate Z-basis measurement, returning estimated ⟨Z⟩."""
p0 = abs(psi[0]) ** 2 # probability of |0⟩
outcomes = np.random.binomial(num_shots, p0)
# ⟨Z⟩ = P(0) - P(1) = 2*P(0) - 1
return 2 * outcomes / num_shots - 1
def tomography(psi_true, num_shots=10000):
"""Perform single-qubit state tomography with shot noise.
Returns the reconstructed density matrix.
"""
H = np.array([[1, 1], [1, -1]]) / np.sqrt(2)
Sdg = np.array([[1, 0], [0, -1j]]) # S†
# Measure ⟨Z⟩ directly
exp_z = simulate_measurement(psi_true, num_shots)
# Measure ⟨X⟩: apply H then measure Z
psi_x = H @ psi_true
exp_x = simulate_measurement(psi_x, num_shots)
# Measure ⟨Y⟩: apply S†H then measure Z
psi_y = Sdg @ H @ psi_true
exp_y = simulate_measurement(psi_y, num_shots)
# Reconstruct density matrix
rho_reconstructed = (I2 + exp_x * X + exp_y * Y + exp_z * Z) / 2
return rho_reconstructed, (exp_x, exp_y, exp_z)
def fidelity(rho, psi):
"""Compute fidelity F = ⟨ψ|ρ|ψ⟩ between a density matrix and a pure state."""
return float(np.real(psi.conj() @ rho @ psi))
# True state: Ry(π/3)|0⟩ (tilted toward +x)
theta = np.pi / 3
psi_true = np.array([np.cos(theta / 2), np.sin(theta / 2)], dtype=complex)
# True Bloch vector
rho_true = np.outer(psi_true, psi_true.conj())
true_bloch = (
float(np.real(np.trace(rho_true @ X))),
float(np.real(np.trace(rho_true @ Y))),
float(np.real(np.trace(rho_true @ Z))),
)
print(f"True Bloch vector: ({true_bloch[0]:.4f}, {true_bloch[1]:.4f}, {true_bloch[2]:.4f})")
# Perform tomography
np.random.seed(42)
rho_recon, measured_bloch = tomography(psi_true, num_shots=10000)
print(f"Measured Bloch vector: ({measured_bloch[0]:.4f}, {measured_bloch[1]:.4f}, {measured_bloch[2]:.4f})")
print(f"Fidelity: {fidelity(rho_recon, psi_true):.4f}")
# With 10,000 shots, typical fidelity is > 0.99
The accuracy improves with more shots. With N shots per basis, the standard error on each Pauli expectation value scales as 1/√N. For 10,000 shots, the uncertainty is about 0.01 per component, giving fidelities above 0.99 for most states.
Noise channels on the Bloch sphere
Noise transforms pure states on the sphere surface into mixed states in the interior. Each noise channel has a characteristic geometric effect on the Bloch vector.
Depolarising channel
The depolarising channel with parameter p replaces the state with the maximally mixed state with probability p:
ρ → (1 - p)ρ + p · I/2
On the Bloch vector, this is uniform shrinkage:
r_out = (1 - p) · r_in
Every component shrinks by the same factor. The Bloch sphere contracts uniformly toward the origin.
Amplitude damping
Amplitude damping models energy relaxation (T1 decay). A qubit in |1⟩ decays to |0⟩ with probability γ. The Bloch vector transforms as:
x_out = √(1 - γ) · x_in
y_out = √(1 - γ) · y_in
z_out = γ + (1 - γ) · z_in
This is not uniform shrinkage. The equatorial components shrink, and the z component is biased toward +1 (the |0⟩ state). The sphere becomes an ellipsoid that drifts toward the north pole.
Phase damping (dephasing)
Phase damping models T2 decay (loss of coherence without energy loss). Only the off-diagonal elements of the density matrix decay:
x_out = √(1 - λ) · x_in
y_out = √(1 - λ) · y_in
z_out = z_in
The z component is untouched. The equatorial plane shrinks while the poles remain fixed. The sphere becomes a pancake that collapses onto the z-axis.
Here is code implementing all three channels and verifying them against the Kraus operator formalism:
import numpy as np
# Pauli matrices
I2 = np.eye(2, dtype=complex)
X = np.array([[0, 1], [1, 0]], dtype=complex)
Y = np.array([[0, -1j], [1j, 0]], dtype=complex)
Z = np.array([[1, 0], [0, -1]], dtype=complex)
def bloch_from_rho(rho):
x = float(np.real(np.trace(rho @ X)))
y = float(np.real(np.trace(rho @ Y)))
z = float(np.real(np.trace(rho @ Z)))
return np.array([x, y, z])
def rho_from_bloch(r):
return (I2 + r[0] * X + r[1] * Y + r[2] * Z) / 2
# --- Depolarising channel ---
def depolarising_bloch(r, p):
return (1 - p) * r
def depolarising_kraus(rho, p):
K0 = np.sqrt(1 - 3 * p / 4) * I2
K1 = np.sqrt(p / 4) * X
K2 = np.sqrt(p / 4) * Y
K3 = np.sqrt(p / 4) * Z
return K0 @ rho @ K0.conj().T + K1 @ rho @ K1.conj().T + \
K2 @ rho @ K2.conj().T + K3 @ rho @ K3.conj().T
# --- Amplitude damping ---
def amplitude_damping_bloch(r, gamma):
return np.array([
np.sqrt(1 - gamma) * r[0],
np.sqrt(1 - gamma) * r[1],
gamma + (1 - gamma) * r[2]
])
def amplitude_damping_kraus(rho, gamma):
K0 = np.array([[1, 0], [0, np.sqrt(1 - gamma)]])
K1 = np.array([[0, np.sqrt(gamma)], [0, 0]])
return K0 @ rho @ K0.conj().T + K1 @ rho @ K1.conj().T
# --- Phase damping ---
def phase_damping_bloch(r, lam):
return np.array([
np.sqrt(1 - lam) * r[0],
np.sqrt(1 - lam) * r[1],
r[2]
])
def phase_damping_kraus(rho, lam):
K0 = np.array([[1, 0], [0, np.sqrt(1 - lam)]])
K1 = np.array([[0, 0], [0, np.sqrt(lam)]])
return K0 @ rho @ K0.conj().T + K1 @ rho @ K1.conj().T
# Test state: |+⟩ rotated by Ry(π/3), Bloch vector ≈ (0.866, 0, 0.5)
theta = np.pi / 3
psi = np.array([np.cos(theta / 2), np.sin(theta / 2)], dtype=complex)
rho = np.outer(psi, psi.conj())
r_in = bloch_from_rho(rho)
print(f"Input Bloch vector: ({r_in[0]:.4f}, {r_in[1]:.4f}, {r_in[2]:.4f})")
for name, bloch_fn, kraus_fn, param in [
("Depolarising (p=0.3)", depolarising_bloch, depolarising_kraus, 0.3),
("Amplitude damping (γ=0.4)", amplitude_damping_bloch, amplitude_damping_kraus, 0.4),
("Phase damping (λ=0.5)", phase_damping_bloch, phase_damping_kraus, 0.5),
]:
r_bloch = bloch_fn(r_in, param)
rho_kraus = kraus_fn(rho, param)
r_kraus = bloch_from_rho(rho_kraus)
print(f"\n{name}:")
print(f" Bloch formula: ({r_bloch[0]:.4f}, {r_bloch[1]:.4f}, {r_bloch[2]:.4f})")
print(f" Kraus result: ({r_kraus[0]:.4f}, {r_kraus[1]:.4f}, {r_kraus[2]:.4f})")
print(f" Match: {np.allclose(r_bloch, r_kraus)}")
Multiple qubits and the partial Bloch sphere
The Bloch sphere represents a single qubit. For a system of two or more qubits, there is no single 3D sphere that captures the full state. A 2-qubit system has 15 real parameters (a 4×4 density matrix minus normalisation), far more than the 2 parameters of a sphere.
However, you can always trace out one qubit and examine the “reduced Bloch vector” of the other. This gives the local state of that qubit, discarding all information about correlations.
For an entangled Bell state, tracing out one qubit gives the maximally mixed state:
import numpy as np
# Bell state |Φ+⟩ = (|00⟩ + |11⟩) / √2
bell = np.array([1, 0, 0, 1], dtype=complex) / np.sqrt(2)
rho_bell = np.outer(bell, bell.conj())
# Partial trace over qubit 1 to get reduced state of qubit 0
# For a 4x4 density matrix in the |00⟩, |01⟩, |10⟩, |11⟩ basis:
# rho_0 = Tr_1(rho) has elements rho_0[i,j] = sum_k rho[2*i+k, 2*j+k]
def partial_trace_1(rho_ab):
"""Trace out qubit 1 (second qubit) from a 2-qubit density matrix."""
rho_a = np.zeros((2, 2), dtype=complex)
for i in range(2):
for j in range(2):
for k in range(2):
rho_a[i, j] += rho_ab[2 * i + k, 2 * j + k]
return rho_a
rho_0 = partial_trace_1(rho_bell)
print("Reduced density matrix of qubit 0:")
print(np.round(rho_0, 4))
# [[0.5, 0. ],
# [0. , 0.5]] (this is I/2, the maximally mixed state)
# Bloch vector
X = np.array([[0, 1], [1, 0]], dtype=complex)
Y = np.array([[0, -1j], [1j, 0]], dtype=complex)
Z = np.array([[1, 0], [0, -1]], dtype=complex)
bx = float(np.real(np.trace(rho_0 @ X)))
by = float(np.real(np.trace(rho_0 @ Y)))
bz = float(np.real(np.trace(rho_0 @ Z)))
print(f"Bloch vector: ({bx:.4f}, {by:.4f}, {bz:.4f})")
# (0, 0, 0) (the center of the sphere)
# Contrast with a product state |0⟩ ⊗ |+⟩
zero = np.array([1, 0], dtype=complex)
plus = np.array([1, 1], dtype=complex) / np.sqrt(2)
product = np.kron(zero, plus)
rho_product = np.outer(product, product.conj())
rho_0_product = partial_trace_1(rho_product)
bx = float(np.real(np.trace(rho_0_product @ X)))
by = float(np.real(np.trace(rho_0_product @ Y)))
bz = float(np.real(np.trace(rho_0_product @ Z)))
print(f"\nProduct state |0⟩⊗|+⟩, qubit 0 Bloch vector: ({bx:.4f}, {by:.4f}, {bz:.4f})")
# (0, 0, 1) (pure state |0⟩, on the surface of the sphere)
This is a fundamental result: maximal entanglement produces maximal mixedness locally. If qubit 0 is maximally entangled with qubit 1, then qubit 0 alone looks completely random. Its Bloch vector is zero, sitting at the center of the sphere. Conversely, if a qubit’s reduced state is pure (Bloch vector on the surface), it cannot be entangled with anything else.
Gate error as Bloch sphere deviation
Real quantum hardware does not apply gates perfectly. A common systematic error is over-rotation: instead of Rz(θ), the hardware applies Rz(θ(1 + ε)) for some small fixed ε. This error accumulates linearly with the number of gates.
After n applications of the flawed gate, the intended rotation angle is nθ but the actual angle is nθ(1 + ε). The angular deviation is nθε. For ε = 0.01 and n = 100 applications of Rz(π/4), the accumulated over-rotation is 100 · (π/4) · 0.01 ≈ 0.785 radians, which is about 45°. That is an enormous error.
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
I2 = np.eye(2, dtype=complex)
sz = np.array([[1, 0], [0, -1]], dtype=complex)
def rz(theta):
return np.cos(theta / 2) * I2 - 1j * np.sin(theta / 2) * sz
def state_to_bloch(psi):
rho = np.outer(psi, psi.conj())
sx = np.array([[0, 1], [1, 0]], dtype=complex)
sy = np.array([[0, -1j], [1j, 0]], dtype=complex)
x = float(np.real(np.trace(rho @ sx)))
y = float(np.real(np.trace(rho @ sy)))
z = float(np.real(np.trace(rho @ sz)))
return np.array([x, y, z])
# Start in |+⟩
plus = np.array([1, 1], dtype=complex) / np.sqrt(2)
theta = np.pi / 4 # intended rotation per gate
epsilon = 0.01 # 1% over-rotation
n_gates = 100
# Track ideal and faulty trajectories
ideal_points = [state_to_bloch(plus)]
faulty_points = [state_to_bloch(plus)]
psi_ideal = plus.copy()
psi_faulty = plus.copy()
for i in range(n_gates):
psi_ideal = rz(theta) @ psi_ideal
psi_faulty = rz(theta * (1 + epsilon)) @ psi_faulty
ideal_points.append(state_to_bloch(psi_ideal))
faulty_points.append(state_to_bloch(psi_faulty))
ideal_points = np.array(ideal_points)
faulty_points = np.array(faulty_points)
# Compute angular deviation at each step
deviations = []
for i in range(n_gates + 1):
dot = np.clip(np.dot(ideal_points[i], faulty_points[i]), -1, 1)
deviations.append(np.arccos(dot))
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(14, 5))
# Left: accumulated error
ax1.plot(range(n_gates + 1), np.degrees(deviations))
ax1.set_xlabel("Number of gates")
ax1.set_ylabel("Angular deviation (degrees)")
ax1.set_title(f"Systematic over-rotation: ε = {epsilon}")
ax1.grid(True)
# Right: 3D trajectory comparison
ax2 = fig.add_subplot(122, projection="3d")
u = np.linspace(0, 2 * np.pi, 30)
v = np.linspace(0, np.pi, 20)
xs = np.outer(np.cos(u), np.sin(v))
ys = np.outer(np.sin(u), np.sin(v))
zs = np.outer(np.ones_like(u), np.cos(v))
ax2.plot_wireframe(xs, ys, zs, alpha=0.08, color="gray")
# Plot every 5th point for clarity
step = 5
ax2.plot(ideal_points[::step, 0], ideal_points[::step, 1],
ideal_points[::step, 2], "b.-", label="Ideal", alpha=0.7)
ax2.plot(faulty_points[::step, 0], faulty_points[::step, 1],
faulty_points[::step, 2], "r.-", label="Faulty", alpha=0.7)
ax2.set_title("Ideal vs faulty trajectories")
ax2.legend()
plt.tight_layout()
plt.savefig("gate_error_accumulation.png", dpi=150)
plt.show()
print(f"After {n_gates} gates:")
print(f" Accumulated over-rotation: {n_gates * theta * epsilon:.4f} rad "
f"= {np.degrees(n_gates * theta * epsilon):.1f}°")
print(f" Final angular deviation: {np.degrees(deviations[-1]):.1f}°")
This linear error accumulation is why techniques like dynamical decoupling (inserting corrective pulses to cancel systematic errors) and randomized compiling (randomly choosing equivalent gate decompositions to convert coherent errors into incoherent noise) are essential for running deep circuits on real hardware.
Solovay-Kitaev and gate approximation
The Solovay-Kitaev theorem states that any single-qubit unitary can be approximated to precision ε using O(log^c(1/ε)) gates from a universal gate set, where c is a small constant (approximately 3.97 for the original proof, improved to about 1.44 in recent versions). The standard universal gate set is Clifford+T: {H, S, T, and their inverses}.
This is a remarkable result. It means the Bloch sphere can be reached to arbitrary precision using only a finite set of discrete gates, and the cost grows only polylogarithmically with the desired precision.
Here is a simple demonstration that approximates Rz(π/5) using sequences of H and T gates. This is not the full Solovay-Kitaev algorithm (which is recursive and more sophisticated), but it illustrates the basic idea of gate set approximation:
import numpy as np
from itertools import product
# Define basic gates
H = np.array([[1, 1], [1, -1]]) / np.sqrt(2)
T = np.array([[1, 0], [0, np.exp(1j * np.pi / 4)]])
Tdg = T.conj().T # T†
S = T @ T # S = T²
# Target: Rz(π/5)
target_angle = np.pi / 5
Rz_target = np.array([
[np.exp(-1j * target_angle / 2), 0],
[0, np.exp(1j * target_angle / 2)]
])
def operator_distance(U, V):
"""Compute the operator norm distance ||U - V||.
Uses the fact that for unitaries, ||U - V|| = max singular value of (U - V).
"""
return np.linalg.norm(U - V, ord=2)
def matrix_up_to_phase(U, V):
"""Distance between U and V, minimised over global phase."""
# Try matching the (0,0) elements' phases
if abs(U[0, 0]) > 1e-10 and abs(V[0, 0]) > 1e-10:
phase = U[0, 0] / V[0, 0]
phase = phase / abs(phase)
return operator_distance(U, phase * V)
return operator_distance(U, V)
# Brute-force search over sequences of H, T, T† of increasing length
gate_dict = {"H": H, "T": T, "T†": Tdg}
gate_names = list(gate_dict.keys())
best_dist = float("inf")
best_seq = None
for length in range(1, 14):
for seq in product(gate_names, repeat=length):
# Build the unitary
U = np.eye(2, dtype=complex)
for name in seq:
U = gate_dict[name] @ U
dist = matrix_up_to_phase(Rz_target, U)
if dist < best_dist:
best_dist = dist
best_seq = seq
if dist < 0.01:
break
if best_dist < 0.01:
break
print(f"Target: Rz(π/5)")
print(f"Best sequence (length {len(best_seq)}): {' · '.join(best_seq)}")
print(f"Operator distance (up to global phase): {best_dist:.6f}")
# Verify
U_approx = np.eye(2, dtype=complex)
for name in best_seq:
U_approx = gate_dict[name] @ U_approx
print(f"\nTarget matrix:\n{np.round(Rz_target, 4)}")
print(f"Approximation:\n{np.round(U_approx, 4)}")
Note that this brute-force search is only feasible for short sequences. The actual Solovay-Kitaev algorithm uses a recursive approach that is far more efficient, finding near-optimal sequences in polynomial time (in log(1/ε)).
Geometric interpretation of universality
The Bloch sphere provides a clean way to understand why certain gate sets are universal and others are not.
H alone generates a finite group. The Hadamard gate has order 2 (H² = I). Starting from |0⟩ and applying H repeatedly, you can only reach two points: (0, 0, 1) and (1, 0, 0). Starting from other states, H still only generates a finite orbit. The group generated by H alone is too small to be dense in SU(2).
T alone generates a finite group on the z-axis. Since T = Rz(π/4), applying T repeatedly from a fixed starting point generates at most 8 distinct points (since 8 · π/4 = 2π). Alone, T cannot leave the z-rotation orbit of whatever state you start in.
H and T together generate a dense subgroup. The key insight is that H and T rotate around different axes. H rotates around the (x+z)/√2 axis, while T rotates around the z-axis. When two rotations have axes that are not aligned, and the rotation angles are irrational multiples of π (which is the case for the effective combined rotations), their compositions generate an infinite, dense subgroup of SU(2). This means any point on the Bloch sphere can be reached arbitrarily closely.
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
H = np.array([[1, 1], [1, -1]]) / np.sqrt(2)
T = np.array([[1, 0], [0, np.exp(1j * np.pi / 4)]])
sx = np.array([[0, 1], [1, 0]], dtype=complex)
sy = np.array([[0, -1j], [1j, 0]], dtype=complex)
sz = np.array([[1, 0], [0, -1]], dtype=complex)
def state_to_bloch(psi):
rho = np.outer(psi, psi.conj())
x = float(np.real(np.trace(rho @ sx)))
y = float(np.real(np.trace(rho @ sy)))
z = float(np.real(np.trace(rho @ sz)))
return (x, y, z)
# Apply random sequences of H and T to |0⟩
psi0 = np.array([1, 0], dtype=complex)
points = set()
rng = np.random.default_rng(42)
psi = psi0.copy()
for _ in range(5000):
gate = H if rng.random() < 0.5 else T
psi = gate @ psi
psi = psi / np.linalg.norm(psi) # numerical stability
b = state_to_bloch(psi)
points.add((round(b[0], 4), round(b[1], 4), round(b[2], 4)))
points = np.array(list(points))
print(f"Unique Bloch sphere points visited: {len(points)}")
fig = plt.figure(figsize=(8, 8))
ax = fig.add_subplot(111, projection="3d")
# Wireframe sphere
u = np.linspace(0, 2 * np.pi, 30)
v = np.linspace(0, np.pi, 20)
xs = np.outer(np.cos(u), np.sin(v))
ys = np.outer(np.sin(u), np.sin(v))
zs = np.outer(np.ones_like(u), np.cos(v))
ax.plot_wireframe(xs, ys, zs, alpha=0.05, color="gray")
ax.scatter(points[:, 0], points[:, 1], points[:, 2],
s=1, alpha=0.5, c="blue")
ax.set_title("Dense coverage of Bloch sphere by H + T sequences")
plt.tight_layout()
plt.savefig("ht_universality.png", dpi=150)
plt.show()
For full quantum universality (not just single-qubit operations), you also need a two-qubit entangling gate. Adding CNOT to {H, T} gives the standard universal gate set. H and T handle arbitrary single-qubit rotations (densely), while CNOT creates the entanglement that makes multi-qubit computation possible.
Mixed states: the interior of the sphere
So far we have considered pure states, which live on the surface of the sphere. A mixed state is described by a density matrix ρ and lives in the interior. The Bloch vector for a mixed state is:
r = (⟨X⟩, ⟨Y⟩, ⟨Z⟩) = (Tr(ρX), Tr(ρY), Tr(ρZ))
The length |r| = 1 for pure states, |r| < 1 for mixed states, |r| = 0 for the maximally mixed state I/2.
Physically, decoherence shrinks the Bloch vector toward the origin. T2 dephasing shrinks the x and y components (randomises the phase φ), while T1 relaxation drives the z component toward −1 (thermal ground state) or, more precisely, toward the thermal equilibrium point determined by the qubit temperature and frequency.
Common mistakes
These five mistakes appear frequently in quantum computing courses and textbooks. Being aware of them saves significant debugging time.
1. Confusing θ with θ/2 in the parametrisation.
The Bloch sphere parametrisation uses cos(θ/2) and sin(θ/2), not cos(θ) and sin(θ). At θ = π, the state is |1⟩ (the south pole), not at θ = π/2. The state at the equator (θ = π/2) has equal superposition amplitudes cos(π/4) = sin(π/4) = 1/√2. If you write cos(π) = −1 as a coefficient, you have made this error.
2. Treating the Bloch sphere as describing 2 qubits.
The Bloch sphere represents exactly one qubit. A 2-qubit system requires 15 real parameters to fully describe (a 4×4 Hermitian matrix with unit trace). You cannot draw entanglement on a single Bloch sphere. For two qubits, you need two reduced Bloch vectors plus the 9-component correlation tensor ⟨σ_i ⊗ σ_j⟩, giving 3 + 3 + 9 = 15 parameters total.
3. Representing mixed states as points on the sphere surface.
Mixed states live strictly inside the sphere (|r| < 1). If your qubit has undergone any decoherence, its state is mixed, and placing it on the sphere surface is incorrect. The purity of the state is Tr(ρ²) = (1 + |r|²)/2, so a Bloch vector of length 0.5 corresponds to a purity of 0.625, not a pure state.
4. Thinking φ has period 4π.
The angle θ uses θ/2 in the parametrisation, creating the 720° (4π) periodicity for polar rotations. But φ appears directly as e^(iφ), so it has ordinary 2π periodicity. The state |ψ(φ + 2π)⟩ = |ψ(φ)⟩ exactly. Only rotations that change θ exhibit the double-cover behaviour.
5. Assuming the Bloch sphere works for qutrits or higher dimensions.
A qutrit (3-level system) has 8 real parameters (a 3×3 Hermitian matrix with unit trace, minus one for normalisation). This is not a sphere in any meaningful sense. The space of qutrit states is a complicated 8-dimensional convex body. The Gell-Mann matrices replace the Pauli matrices, and the “generalised Bloch vector” lives in ℝ⁸, but its allowed region is not a ball. The Bloch sphere geometry is specific to the qubit.
Why this matters for gate design
Understanding the Bloch sphere geometry helps you:
-
Choose the right gate sequence: To move from |0⟩ to |i⟩ (the +y axis), you rotate from the north pole to the equator. One way: Ry(π/2). Another: H then S. Both rotate to the same point.
-
Minimise circuit depth: Any single-qubit operation can be decomposed as three Euler rotations: Rz(α) Ry(β) Rz(γ). Transpilers use this to decompose arbitrary unitaries into native gates.
-
Understand gate errors: On real hardware, control pulse imperfections cause the rotation axis to deviate slightly from the intended axis, or the rotation angle to be slightly off. These are systematic rotation errors on the Bloch sphere.
-
Interpret noise: T1 relaxation pushes the Bloch vector toward the south pole. T2 dephasing shrinks the equatorial plane without affecting the z-axis. You can see both effects in the shrinkage of the Bloch sphere after running circuits.
Related resources
- Interactive Bloch Sphere Simulator - visualise these rotations in real time
- Quantum Gates Explained - matrix representations of common gates
- Density Matrices and Mixed States - the full story on the sphere interior
- Glossary: Bloch Sphere
Was this tutorial helpful?