MRI is, at its core, a quantum measurement device. One of the few that operates at room temperature on bulk matter and still yields genuinely quantum-mechanical signal.
Atomic nuclei with an odd number of protons or neutrons carry an intrinsic quantum property called spin angular momentum, not because they are literally spinning like a top but because they possess a quantised internal angular momentum with no classical equivalent.
Think of it as a built-in immutable “twistiness” that each nucleus carries. This spin makes the nucleus behave like a tiny magnetic dipole; it has a north and a south pole, and it responds to external magnetic fields. Hydrogen-1 (a single proton) is the nucleus MRI exploits because it is the most abundant nucleus in the human body and has the strongest magnetic response of any stable nucleus.
Spin is quantised meaning it doesn’t take continuous values. For ¹H, the spin quantum number s = ½, which gives exactly two allowed orientations in an applied magnetic field:
|↑⟩ → mₛ = +½ (spin-up, "parallel")
|↓⟩ → mₛ = −½ (spin-down, "anti-parallel")In isolation, no external field in these two spin states have identical energy. They are degenerate. Apply a magnetic field and that degeneracy breaks. The two orientations now have different energies, with the parallel state sitting slightly lower. This field-induced energy splitting is the Zeeman effect, named after the Dutch physicist who first observed analogous splitting in atomic spectral lines in 1896. The stronger the applied field, the wider the gap.
The MRI scanner’s static field B₀ (typically 1.5–7 T) does exactly this:
ΔE = ℏ γ B₀The constant γ here is the gyromagnetic ratio, a nucleus-specific number that captures how strongly that nucleus’s spin couples to an external magnetic field. It sets the proportionality between field strength and precession frequency. A higher γ means the nucleus responds more vigorously to the same field, like a tighter gear ratio. For ¹H, γ/2π ≈ 42.577 MHz/T. This is why you can’t directly substitute a different nucleus without retuning the entire RF system - carbon-13, for instance, has γ/2π ≈ 10.7 MHz/T, nearly four times lower.
At 3T, the energy gap ΔE corresponds to a Larmor precession frequency of ~127.7 MHz, squarely in the radiofrequency band. This precession is the quantum mechanical analogue of a gyroscope wobbling under gravity. The spin axis rotates around B₀ at this characteristic frequency. It is the resonant frequency you tune the scanner to, and the frequency at which you must broadcast energy to actually perturb the spin states.
You’re not manipulating individual qubits, you’re working with an ensemble of ~10²³ nuclear spins, and thermal fluctuations at 310 K mean the population difference between |↑⟩ and |↓⟩ is tiny.
The Boltzmann distribution gives the fractional excess of spin-up nuclei:
ΔN/N ≈ ℏγB₀ / (2kT)At 3T and 310K, this is roughly 1 in 10⁶ - a vanishingly small polarisation. But with ~10²³ spins in a voxel, the coherent classical observable, the bulk magnetisation vector M is large enough to induce a measurable EMF in a receiver coil.
M aligns with B₀ at thermal equilibrium:
M₀ = ρ · (γ²ℏ²B₀) / (4kT)where ρ is the spin density. This is the quantity you manipulate and detect.
Once you accept the ensemble picture, the quantum dynamics of M are governed by the Bloch equations, a semi-classical description that emerges from the density matrix formalism.
Each equation has two structurally distinct parts doing different jobs.
The code below works in the rotating frame, a coordinate system that itself rotates at the Larmor frequency ω₀. This is the classical equivalent of moving to the interaction picture: in this frame, B₀ vanishes from the equations entirely, and you only track the slow envelope of the magnetisation dynamics relative to the carrier frequency. The delta_omega term captures any off-resonance deviation from ω₀; when on-resonance it is zero, and the equations simplify to pure T1/T2 decay after any applied B1 pulse.
dMx/dt = γ(My·Bz − Mz·By) − Mx/T2
dMy/dt = γ(Mz·Bx − Mx·Bz) − My/T2
dMz/dt = γ(Mx·By − My·Bx) − (Mz − M₀)/T1T1 (spin-lattice relaxation) governs how quickly Mz recovers to equilibrium. Think of it as the timescale on which the spin ensemble dumps its excess energy into the surrounding molecular environment - the “lattice”, and returns to its low-energy ground state. T2 (spin-spin relaxation) governs transverse dephasing. After the RF pulse tips spins into the transverse plane, they all start precessing together but local magnetic field variations (from neighbouring nuclei, susceptibility gradients, magnet imperfections) cause individual spins to precess at slightly different rates. They fan out in phase, the coherent signal cancels, and Mxy decays - even before any energy has been exchanged with the lattice. T2 ≤ T1, always, because dephasing is faster than thermalisation.
The signal you detect is the oscillating transverse magnetisation Mxy = Mx + iMy, picked up by Faraday induction.
Below is a clean numerical integration of the Bloch equations, simulating a 90° RF pulse followed by free induction decay (FID):
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp
# ── Tissue and scanner parameters
gamma = 2 * np.pi * 42.577e6 # rad/s/T (¹H gyromagnetic ratio)
B0 = 3.0 # Tesla
omega0 = gamma * B0 # Larmor frequency (rad/s)
T1 = 1.0 # seconds (e.g. white matter at 3T)
T2 = 0.08 # seconds
M0 = 1.0 # normalised equilibrium magnetisation
# In the rotating frame at ω₀, B0 drops out; only off-resonance Δω and B1 remain.
delta_omega = 0.0 # on-resonance
def bloch_rotating(t, M, B1x, B1y, delta_omega, T1, T2, M0):
Mx, My, Mz = M
dMx = delta_omega * My - Mx / T2 + gamma * My * 0 # simplified
dMy = -delta_omega * Mx - My / T2 - gamma * Mz * B1x
dMz = gamma * My * B1x - (Mz - M0) / T1
# Full rotating-frame Bloch:
dMx = delta_omega*My + gamma*Mz*B1y - Mx/T2
dMy = -delta_omega*Mx - gamma*Mz*B1x - My/T2
dMz = gamma*(Mx*B1y - My*B1x) - (Mz - M0)/T1
return [dMx, dMy, dMz]
# ── Phase 1: 90° hard RF pulse
# Flip angle α = γ·B1·τ_pulse → B1 = α / (γ·τ_pulse)
alpha = np.pi / 2 # 90° flip
tau_rf = 1e-3 # 1 ms pulse duration
B1x_on = alpha / (gamma * tau_rf)
B1y_on = 0.0
t_rf = np.linspace(0, tau_rf, 500)
sol_rf = solve_ivp(
bloch_rotating,
[0, tau_rf],
[0.0, 0.0, M0], # start at thermal equilibrium
t_eval=t_rf,
args=(B1x_on, B1y_on, delta_omega, T1, T2, M0),
max_step=tau_rf/200
)
M_after_pulse = sol_rf.y[:, -1] # state immediately post-pulse
# ── Phase 2: Free induction decay (FID)
t_fid = np.linspace(0, 0.5, 5000)
sol_fid = solve_ivp(
bloch_rotating,
[0, 0.5],
M_after_pulse,
t_eval=t_fid,
args=(0.0, 0.0, delta_omega, T1, T2, M0),
max_step=1e-4
)
Mx_fid, My_fid, Mz_fid = sol_fid.y
# ── Plot
fig, axes = plt.subplots(2, 1, figsize=(10, 7), sharex=False)
# RF pulse dynamics
ax = axes[0]
ax.plot(sol_rf.t * 1e3, sol_rf.y[0], label='Mx')
ax.plot(sol_rf.t * 1e3, sol_rf.y[1], label='My')
ax.plot(sol_rf.t * 1e3, sol_rf.y[2], label='Mz')
ax.set_xlabel('Time (ms)')
ax.set_ylabel('Magnetisation (normalised)')
ax.set_title('90° RF Pulse in Rotating Frame')
ax.legend(); ax.grid(True, alpha=0.3)
# FID
ax = axes[1]
signal = np.sqrt(Mx_fid**2 + My_fid**2)
ax.plot(sol_fid.t * 1e3, signal, color='steelblue', label='|Mxy| (FID envelope)')
ax.plot(sol_fid.t * 1e3, M0 * np.exp(-sol_fid.t / T2),
'r--', alpha=0.7, label=f'M₀·exp(−t/T2), T2={T2*1e3:.0f}ms')
ax.plot(sol_fid.t * 1e3, M0 * (1 - np.exp(-sol_fid.t / T1)),
'g--', alpha=0.7, label=f'Mz recovery (T1={T1*1e3:.0f}ms)')
ax.plot(sol_fid.t * 1e3, Mz_fid, 'g', alpha=0.4)
ax.set_xlabel('Time (ms)')
ax.set_ylabel('Magnetisation (normalised)')
ax.set_title('Free Induction Decay')
ax.legend(); ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig('bloch_simulation.png', dpi=150)
plt.show()Running this will show the magnetisation tipping from +Z into the transverse plane during the pulse, then the characteristic exponential T2 decay of the FID signal - the raw quantum signal before any spatial encoding.
The scanner detects signal from the entire excited volume simultaneously. To reconstruct an image, it uses linear gradient fields, small spatially varying perturbations superimposed on B₀:
B(r, t) = B₀ + G(t) · rThis makes the local Larmor frequency position-dependent: ω(r, t) = γ(B₀ + G·r). A spin at position r accumulates phase:
φ(r, t) = γ ∫₀ᵗ G(τ)·r dτ = k(t) · rwhere k(t) = γ∫G(τ)dτ is the k-space trajectory - a spatial frequency coordinate. The total received signal is:
S(k) = ∫ ρ(r) · e^{i k·r} drThis is a Fourier transform. The entire field of view is sampled in k-space, and an inverse FFT yields the image. Gradient waveform design is fundamentally k-space trajectory design.
import numpy as np
def simulate_kspace_readout(rho, fov=0.24, N=256, gamma_bar=42.577e6):
"""
Simulate a Cartesian EPI k-space acquisition.
rho : 2D array, spin density (the 'true' image)
fov : field of view in metres
N : matrix size
Returns: k-space data S(kx, ky)
"""
# k-space is just the 2D FFT of spin density
# (ignoring relaxation, off-resonance, and noise for clarity)
kspace = np.fft.fftshift(np.fft.fft2(rho))
return kspace
def reconstruct_image(kspace):
"""Standard inverse FFT reconstruction."""
img = np.fft.ifft2(np.fft.ifftshift(kspace))
return np.abs(img)
# Example: simulate a simple phantom and round-trip through k-space
N = 256
phantom = np.zeros((N, N))
cy, cx = N // 2, N // 2
Y, X = np.ogrid[:N, :N]
# Elliptical brain-like region
phantom[(X-cx)**2/90**2 + (Y-cy)**2/110**2 < 1] = 1.0
# Inner structure (ventricles, CSF - long T2)
phantom[(X-cx)**2/20**2 + (Y-cy)**2/25**2 < 1] = 0.3
kspace = simulate_kspace_readout(phantom)
# Demonstrate k-space undersampling (acceleration factor R=4)
kspace_undersampled = kspace.copy()
kspace_undersampled[::4, :] = 0 # zero every 4th ky line
img_full = reconstruct_image(kspace)
img_undersampled = reconstruct_image(kspace_undersampled)
print(f"k-space shape: {kspace.shape}")
print(f"Centre of k-space (low spatial freq) magnitude: {np.abs(kspace[N//2, N//2]):.1f}")
print(f"Periphery (high spatial freq) magnitude: {np.abs(kspace[0, 0]):.4f}")The centre of k-space encodes contrast (overall signal energy); the periphery encodes spatial resolution (edges, fine detail). This asymmetry is central to modern acceleration techniques like partial Fourier and compressed sensing.
Modern scanners exploit the sparsity of MR images in transformed domains (wavelets, finite differences) to reconstruct from heavily undersampled k-space, directly analogous to ideas in quantum state tomography:
# Conceptual sparse recovery via L1 minimisation (compressed sensing)
# In practice: ADMM, BART, or learned unrolled networks
#
# Minimise: ||F_u · m − y||₂² + λ||Ψm||₁
#
# F_u = undersampled Fourier operator
# y = acquired k-space data
# Ψ = sparsifying transform (e.g. wavelet)
# λ = regularisation weight
#
# This is the MR equivalent of L1-regularised quantum state tomography
# both exploit the prior that the true signal has a sparse representation.
import numpy as np
from scipy.optimize import minimize
def cs_recon_1d(y_undersampled, mask, lam=1e-3, n_iter=200):
"""
Toy 1D compressed sensing reconstruction via proximal gradient.
y_undersampled : observed (masked) k-space
mask : boolean array, True where k-space was sampled
"""
N = len(mask)
x = np.zeros(N, dtype=complex)
for _ in range(n_iter):
# Data consistency gradient step
kx = np.fft.fft(x)
residual = np.zeros(N, dtype=complex)
residual[mask] = kx[mask] - y_undersampled[mask]
grad = np.fft.ifft(residual)
x = x - 0.1 * grad
# Soft-threshold (L1 prox) in image domain
magnitude = np.abs(x)
x = np.where(magnitude > lam, x * (magnitude - lam) / magnitude, 0)
return xThis same mathematical machinery, exploiting structure in a high-dimensional Hilbert space appears in quantum error correction syndrome decoding and in variational quantum eigensolvers, though the physical substrate differs entirely.
The fundamental SNR limit of MRI is set by Johnson-Nyquist thermal noise in the receiver coil, not by shot noise or quantum projection noise (because you’re measuring a classical induced voltage, not individual spin projections). The SNR scales as:
SNR ∝ B₀^(7/4) · √(Vvoxel · NEX)The B₀^(7/4) dependence (slightly better than B₀ due to both signal and coil noise scaling) drives the relentless push to higher field strengths - 7T clinical systems are now approved, and 10.5T and 14T research scanners exist. At these fields, the RF wavelength in tissue (λ ≈ c / (f · √ε_r) ~ 12 cm at 7T) becomes comparable to the FOV, producing B1 inhomogeneity - a wave-optics problem requiring parallel transmit arrays and RF shimming to solve.
| MRI Concept | QC Analogue |
|---|---|
| Single ¹H spin in B₀ | Qubit in |0⟩/|1⟩ basis |
| RF pulse (flip angle α, phase φ) | Single-qubit rotation Ry(α), Rz(φ) |
| Larmor precession | Free evolution under Hamiltonian H = −ℏω₀σz/2 |
| T2 dephasing | Decoherence / pure dephasing channel |
| T1 relaxation | Amplitude damping channel |
| Gradient echo | Refocussing after controlled Z-rotation gradient |
| Spin echo (180° pulse) | Hahn echo / dynamical decoupling (single π pulse) |
| CPMG sequence | XY-8 / CPMG dynamical decoupling |
| Bloch sphere | Qubit state space (same object, different context) |
The Bloch sphere representation that quantum computing textbooks introduce is the same mathematical object: the Bloch sphere was first used to describe NMR spin dynamics in 1946 - quantum computing borrowed it.
The key difference is MRI operates on a mixed-state ensemble at high temperature, where individual spin coherence is irrelevant and only bulk M matters. Gate-model quantum computing requires and maintains pure states on individual qubits, with decoherence as the enemy. In MRI, T2 relaxation is a design constraint to work around and in QC, it is the central engineering challenge.