%matplotlib inline
from collections import namedtuple
import matplotlib.pyplot as plt
import numpy as np
import qutip
from qutip import Qobj, Bloch, basis, ket, tensor
qutip.about()
QuTiP: Quantum Toolbox in Python ================================ Copyright (c) QuTiP team 2011 and later. Original developers: R. J. Johansson & P. D. Nation. Current admin team: Alexander Pitchford, Paul D. Nation, Nathan Shammah, Shahnawaz Ahmed, Neill Lambert, and Eric Giguère. Project Manager: Franco Nori. Currently developed through wide collaboration. See https://github.com/qutip for details. QuTiP Version: 4.4.1 Numpy Version: 1.17.1 Scipy Version: 1.2.0 Cython Version: 0.29.13 Matplotlib Version: 3.1.1 Python Version: 3.6.8 Number of CPUs: 2 BLAS Info: OPENBLAS OPENMP Installed: False INTEL MKL Ext: False Platform Info: Linux (x86_64) Installation path: /home/simon/venvs/qutip-measurements-euroscipy-2019/lib/python3.6/site-packages/qutip ============================================================================== Please cite QuTiP in your publication. ============================================================================== For your convenience a bibtex reference can be easily generated using `qutip.cite()`
Define LaTeX commands:
Why this weird quantum mechanics anyway?
See also:
for more on the history of the Stern-Gerlach experiment.
z = np.array([0, 0, 1])
mu = np.array([0, 1, 1]) / np.sqrt(2)
bloch = Bloch()
bloch.zlabel=("z", "")
bloch.add_vectors([z, mu])
bloch.show()
Simulating a simple classical system.
# Simulation of expected results in the classical case
Direction = namedtuple("Direction", ["theta", "phi"])
def random_direction():
""" Generate a random direction. """
# See http://mathworld.wolfram.com/SpherePointPicking.html
r = 0
while r == 0:
x, y, z = np.random.normal(0, 1, 3)
r = np.sqrt(x**2 + y**2 + z**2)
phi = np.arctan2(y, x)
theta = np.arccos(z / r)
return Direction(theta=theta, phi=phi)
def classical_state(d):
""" Prepare a spin state given a direction. """
x = np.sin(d.theta) * np.cos(d.phi)
y = np.sin(d.theta) * np.sin(d.phi)
z = np.cos(d.theta)
return np.array([x, y, z])
classical_up = np.array([0, 0, 1])
def classical_spin(c):
""" Measure the z-component of the spin. """
return classical_up.dot(c)
def classical_stern_gerlach(n):
""" Simulate the Stern-Gerlach experiment """
directions = [random_direction() for _ in range(n)]
atoms = [classical_state(d) for d in directions]
spins = [classical_spin(c) for c in atoms]
return atoms, spins
def plot_classical_results(atoms, spins):
fig = plt.figure(figsize=(18.0, 8.0))
fig.suptitle("Stern-Gerlach Experiment: Classical Outcome", fontsize="xx-large")
ax1 = plt.subplot(1, 2, 1, projection='3d')
ax2 = plt.subplot(1, 2, 2)
b = Bloch(fig=fig, axes=ax1)
b.vector_width = 1
b.vector_color = ["#ff{:x}0ff".format(i, i) for i in range(10)]
b.zlabel = ["$z$", ""]
b.add_vectors(atoms)
b.render(fig=fig, axes=ax1)
ax2.hist(spins)
ax2.set_xlabel("Z-component of spin")
ax2.set_ylabel("# of atoms")
atoms, spins = classical_stern_gerlach(1000)
plot_classical_results(atoms, spins)
def plot_real_vs_actual(spins):
fig = plt.figure(figsize=(18.0, 8.0))
fig.suptitle("Stern-Gerlach Experiment: Real vs Actual", fontsize="xx-large")
ax1 = plt.subplot(1, 2, 1)
ax2 = plt.subplot(1, 2, 2)
ax1.hist([np.random.choice([1, -1]) for _ in spins])
ax1.set_xlabel("Z-component of spin")
ax1.set_ylabel("# of atoms")
ax2.hist(spins)
ax2.set_xlabel("Z-component of spin")
ax2.set_ylabel("# of atoms")
plot_real_vs_actual(spins)
This really is a system with two outcomes.
Something we're familiar with.
A classical bit (bit) is classical system with two states and two outcomes.
0
and 1
0
and 1
class ClassicalBit:
def __init__(self, outcome):
self.outcome = outcome
b0 = heads = ClassicalBit(outcome=0)
b1 = tails = ClassicalBit(outcome=1)
def measure_cbit(cbit):
return cbit.outcome
print("State:\n", b0)
print("Outcome:", measure_cbit(b0))
State: <__main__.ClassicalBit object at 0x7f8dac87dcc0> Outcome: 0
Putting the puzzle pieces together.
A quantum bit (qubit) is quantum system with two basis states and two outcomes.
|0>
and |1>
0
and 1
Basis means "we will construct more states from these later".
b0 = ket("0") # |0>
b1 = ket("1") # |1>
print("State:\n", b1)
State: Quantum object: dims = [[2], [1]], shape = (2, 1), type = ket Qobj data = [[0.] [1.]]
def measure_qbit(qbit):
if qbit == ket("0"):
return 0
if qbit == ket("1"):
return 1
raise NotImplementedError("No clue yet. :)")
print("Outcome:", measure_qbit(b1))
Outcome: 1
The tricky part. :)
Let's try something simple for the other states:
We need:
def plot_real_a_b():
fig = plt.figure(figsize=(18.0, 8.0))
fig.suptitle("Probabilities: Real $a$ and $b$", fontsize="xx-large")
ax = plt.subplot(1, 1, 1)
ax.plot([0, 1], [1, 0])
ax.set_xlabel("$a$")
ax.set_xlim(-0.5, 1.5)
ax.set_ylabel("$b$")
ax.set_ylim(-0.5, 1.5)
plot_real_a_b()
Let's try something slightly more complicated:
np.abs(a)
np.angle(a)
a|0> + b|1>
0
with probability $\abs{a}^2$1
with probability $\abs{b}^2$This leaves only two real numbers:
|a|
and |b|
.a
and b
.Both of these form circles -- and the together the two circles form a sphere!
Named after Felix Block (also the first director of CERN!)
SU(2)
- the state space of qubits, aka the bloch sphere.SO(3)
- the space of directions in 3D, aka a sphere.b = Bloch()
b.show()
b = Bloch()
up = ket("0")
down = ket("1")
b.add_states([up, down])
b.show()
x = (up + down).unit()
y = (up + (0 + 1j) * down).unit()
z = up
b = Bloch()
b.add_states([x, y, z])
b.show()
def plot_bloch(fig, ax, title, states, color_template):
""" Plot some states on the bloch sphere. """
b = Bloch(fig=fig, axes=ax)
ax.set_title(title, y=-0.01)
b.vector_width = 1
b.vector_color = [color_template.format(i * 10) for i in range(len(states))]
b.add_states(states)
b.render(fig=fig, axes=ax)
def plot_multi_blochs(plots):
""" Plot multiple sets of states on bloch spheres. """
fig = plt.figure(figsize=(18.0, 8.0))
fig.suptitle("Bloch Sphere", fontsize="xx-large")
n = len(plots)
axes = [plt.subplot(1, n, i + 1, projection='3d') for i in range(n)]
for i, (title, states, color_template) in enumerate(plots):
plot_bloch(fig, axes[i], title, states, color_template)
up = ket("0")
down = ket("1")
# magnitude_circle = [Qobj([[a], [np.sqrt(1 - a**2)]]) for a in np.linspace(0, 1, 20)]
magnitude_circle = [
(a * up + np.sqrt(1 - a**2) * down)
for a in np.linspace(0, 1, 20)
]
# angular_circle = [Qobj([[np.sqrt(0.5)], [np.sqrt(0.5) * np.e ** (1j * theta)]]) for theta in np.linspace(0, np.pi, 20)]
angular_circle = [
(np.sqrt(0.5) * up + np.sqrt(0.5) * down * np.e ** (1j * theta))
for theta in np.linspace(0, np.pi, 20)
]
plot_multi_blochs([
["Changing relative magnitude", magnitude_circle, "#ff{:02x}ff"],
["Changing relative angle", angular_circle, "#{:02x}ffff"],
])
def measure_qbit(qbit):
a = qbit.full()[0][0] # a
b = qbit.full()[1][0] # b
if np.random.random() <= np.abs(a) ** 2:
return 0
else:
return 1
# |a|**2 + |b|**2 == 1
a = (1 + 0j) / np.sqrt(2)
b = (0 + 1j) / np.sqrt(2)
qbit = a * ket("0") + b * ket("1")
print("State:\n", qbit)
print("Outcome:", measure_qbit(qbit))
State: Quantum object: dims = [[2], [1]], shape = (2, 1), type = ket Qobj data = [[0.70710678+0.j ] [0. +0.70710678j]] Outcome: 0
qbit = (1 * ket("0")) + (1j * ket("1"))
qbit = qbit.unit()
print("State:\n", qbit)
print("Outcome:", measure_qbit(qbit))
State: Quantum object: dims = [[2], [1]], shape = (2, 1), type = ket Qobj data = [[0.70710678+0.j ] [0. +0.70710678j]] Outcome: 0
def component(qbit, direction):
component_op = direction.dag()
a = component_op * qbit
return a[0][0]
def measure_qbit(qbit, direction):
a = component(qbit, direction)
if np.random.random() <= np.abs(a) ** 2:
return 0
else:
return 1
up, down = ket("0"), ket("1")
x, y, z = (up + down).unit(), (up + (0 + 1j) * down).unit(), up
print("State:\n", x)
print("Outcomes:", [measure_qbit(x, direction=up) for _ in range(10)])
State: Quantum object: dims = [[2], [1]], shape = (2, 1), type = ket Qobj data = [[0.70710678] [0.70710678]] Outcomes: [0, 1, 0, 1, 1, 1, 1, 0, 0, 1]
Simulating a simple classical quantum system.
def quantum_state(d):
""" Prepare a spin state given a direction. """
return np.cos(d.theta / 2) * up + np.exp(1j * d.phi) * np.sin(d.theta / 2) * down
up = ket('0')
def quantum_spin(q):
""" Measurement the z-component of the spin. """
a_up = (up.dag() * q).tr()
prob_up = np.abs(a_up) ** 2
return 1 if np.random.uniform(0, 1) <= prob_up else -1
def quantum_stern_gerlach(n):
""" Simulate the Stern-Gerlach experiment """
directions = [random_direction() for _ in range(n)]
atoms = [quantum_state(d) for d in directions]
spins = [quantum_spin(q) for q in atoms]
return atoms, spins
def plot_quantum_results(atoms, spins):
fig = plt.figure(figsize=(18.0, 8.0))
fig.suptitle("Stern-Gerlach Experiment: Quantum Outcome", fontsize="xx-large")
ax1 = plt.subplot(1, 2, 1, projection='3d')
ax2 = plt.subplot(1, 2, 2)
b = Bloch(fig=fig, axes=ax1)
b.vector_width = 1
b.vector_color = ["#{:x}0{:x}0ff".format(i, i) for i in range(10)]
b.add_states(atoms)
b.render(fig=fig, axes=ax1)
ax2.hist(spins)
ax2.set_xlabel("Z-component of spin")
ax2.set_ylabel("# of atoms")
atoms, spins = quantum_stern_gerlach(1000)
plot_quantum_results(atoms, spins)
QuTiP documentation [ https://qutip.org/ ]
Quantum Computing for the Determined by Michael Nielsen [ https://michaelnielsen.org/blog/quantum-computing-for-the-determined/ ]
Quantum Computing StackExchange [ https://quantumcomputing.stackexchange.com/ ]
History of the Stern-Gerlach experiment [ https://plato.stanford.edu/entries/physics-experiment/app5.html ]
Picking a random point on a sphere [ http://mathworld.wolfram.com/SpherePointPicking.html ]