%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.
Why did they decide to do this?
It had become generally accepted that atoms could behave as tiny magnets.
In the classical model of the atom, electrons orbit the nucleus. A charge orbiting in a circle generates a magnetic field -- so atoms should act like tiny magnets and be deflected by an inhomogenous magnetic field.
Stern had the idea that this deflection would help better measure and understand the magnetic moment of the atoms & Gerlach brought the experimental expertise.
Here $\hat \mu$ is the magnetic moment (strength and direction of the atom's magnetic field) and $\hat z$ is the direction in which the magnetic field varies (and which the measurement is made in).
The assumptions that the magnitude of the the atoms magnetic moment $\abs{\mu}$ and the inhomogeneity of the magnetic field are constant allow us to make the last statement that $F_z$ is proportional to the dot product of $\hat \mu$ and $\hat z$.
Now for fun, let's display the $\hat \mu$ and $\hat z$ vectors using QuTiP's Bloch
class. For now consider this class just a handy way to plot vectors -- we will be learning more about what it is later on.
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
A classical bit (bit) is classical system with two states and two outcomes.
0
and 1
.0
produces an outcome which we will also label 0
.1
produces an outcome which we will also label 1
.The state is like a Python object. Measurement is an operation that acts on the state and returns an outcome.
Bit is a portmanteau of binary digit.
Examples:
heads
and tails
.high
and low
(or 1
and 0
).With two bits there are four outcomes. With three bits there are eight outcomes. With n
bits there are 2^n
outcomes.
We need n
bits to represent n
bits.
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".
A quantum bit (qubit) is quantum system with two two basis states and two outcomes.
|0>
and |1>
.|0>
produces an outcome which we will label 0
.|1>
produces an outcome which we will label 1
.a|0> + b|1>
.a
and b
are complex numbers.0
with probability |a|^2
and outcome 1
with probability |b|^2
.Qubit is a portmanteau of quantum bit (it's portmanatees all the way down).
Examples:
heads
and tails
.high
and low
(or 1
and 0
).With two qubits there are four outcomes. With three qubits there are eight outcomes. With n
qubits there are 2^n
outcomes.
We 2^n
(minus 1) complex numbers to represent n
qubits.
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)
This means we can rotate the two vectors on the complex plane as long as we don't change the angle between them. Note that this doesn't change the magnitude of a or b (and thus doesn't change the probabilities).
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!)
Wow, this looks a lot like a direction in space! Making progress!
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 ]