Performing Quantum Measurements in QuTiP

  • EuroSciPy 2019
  • Simon Cross

Outline

  • Why this weird quantum mechanics anyway?
  • Simulating a simple classical system.
  • What is a classical bit?
  • What is a qubit?
  • Measuring a qubit.
  • Simulating a simple quantum experiment in QuTiP.

Meta goals

  • Use QuTiP to try things out.
  • Try understand what we're doing!

Imports: Our tools

In [1]:
%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
In [2]:
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:

  • $\newcommand{\ket}[1]{\left|{#1}\right\rangle}$ $\ket{0}$
  • $\newcommand{\bra}[1]{\left\langle{#1}\right|}$ $\bra{1}$
  • $\newcommand{\abs}[1]{\lvert{#1}\rvert}$ $\abs{x}$

The Stern-Gerlach Experiment

Why this weird quantum mechanics anyway?

\begin{align} F_z & = \mu_z \frac{\partial B}{\partial z} \\ & = ( \hat{\mu} \cdot \hat{z} ) |\mu| \frac{\partial B}{\partial z} \\ & \propto \hat{\mu} \cdot \hat{z} \end{align}
In [3]:
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()

Let's simulate the Stern-Gerlach in Python!

Simulating a simple classical system.

In [4]:
# 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)
In [5]:
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])
In [6]:
classical_up = np.array([0, 0, 1])

def classical_spin(c):
    """ Measure the z-component of the spin. """
    return classical_up.dot(c)
In [7]:
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
In [8]:
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")
In [9]:
atoms, spins = classical_stern_gerlach(1000)
plot_classical_results(atoms, spins)

The actual results

In [10]:
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")
In [11]:
plot_real_vs_actual(spins)

Uh-oh.

This really is a system with two outcomes.

If we align all the spins perpendicular to the z-axis:

  • Our classical simulator would show no force on any of the atoms.
  • The actual result would be the same as for the random spin case.

If we align all the spins with the z-axis:

  • Our classical simulator would show always 1.
  • The actual result would match that of the classical simulator.

It's like tossing a biased coin.

Classical bit

Something we're familiar with.

Classical bit

A classical bit (bit) is classical system with two states and two outcomes.

  • States: 0 and 1
  • Outcomes: 0 and 1
In [12]:
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

Quantum bit

Putting the puzzle pieces together.

Quantum bit

A quantum bit (qubit) is quantum system with two basis states and two outcomes.

  • Basis states: |0> and |1>
  • Outcomes: 0 and 1

Basis means "we will construct more states from these later".

In [13]:
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.]]
In [14]:
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

What about the "in between" states?

The tricky part. :)

We want probabilistic outcomes

Let's try something simple for the other states:

  • $a \ket{0} + b \ket{1}$
  • Probability of outcome $0$: $a$
  • Probability of outcome $1$: $b$

We need:

  • $a + b = 1$
  • $1 >= a >= 0$
  • $1 >= b >= 0$
In [15]:
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)
In [16]:
plot_real_a_b()

Hmm. That didn't work.

Let's try something slightly more complicated:

  • $a \ket{0} + b \ket{1}$
  • $a, b \in \mathbb{C}$
  • Probability of outcome $0$: $\abs{a}^2$
  • Probability of outcome $1$: $\abs{b}^2$

Complex numbers summary

  • $a = x + y \cdot j$
    • Length: $\sqrt{x^2 + y^2}$
    • Angle ($\theta$): $arctan(y / x)$
  • Python:
    • Length: np.abs(a)
    • Angle ($\theta$): np.angle(a)

Complex numbers in terms of length and angle

  • $a = m \cdot cos(t) + m \cdot sin(t) \cdot j$
  • $a = m \cdot e ^ {t \cdot j}$
    • Length: $m$
    • Angle: $t$

Global phase

  • 4 real parameters, minus one from $\abs{a}^2 + \abs{b}^2 = 1$
  • Still one too many.
  • Let's remove one!
  • Declare that only the relative angles of $a$ and $b$ are important.

The qubit states

  • General state: a|0> + b|1>
    • where $a$ and $b$ are complex numbers
    • $\abs{a}^2 + \abs{b}^2 = 1$
    • global phase is irrelevant:
      • $a = a \cdot e^{t \cdot j}, b = b \cdot e^{t \cdot j}$
  • Outcome:
    • 0 with probability $\abs{a}^2$
    • 1 with probability $\abs{b}^2$

Bloch sphere

This leaves only two real numbers:

  • The relative magnitudes of |a| and |b|.
  • The relative angle between 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!)

How to sound smart at parties

  • SU(2) is isomorphic to SO(3)
  • SU(2) - the state space of qubits, aka the bloch sphere.
  • SO(3) - the space of directions in 3D, aka a sphere.
  • The state space of possible qubits is a sphere.

QuTiP has nice tools for visualising Bloch spheres

In [17]:
b = Bloch()
b.show()
In [18]:
b = Bloch()
up = ket("0")
down = ket("1")
b.add_states([up, down])
b.show()
In [19]:
x = (up + down).unit()
y = (up + (0 + 1j) * down).unit()
z = up
b = Bloch()
b.add_states([x, y, z])
b.show()
In [20]:
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)
In [21]:
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)
In [22]:
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)
]
In [23]:
plot_multi_blochs([
    ["Changing relative magnitude", magnitude_circle, "#ff{:02x}ff"],
    ["Changing relative angle", angular_circle, "#{:02x}ffff"],
])

Measuring a general state

In [24]:
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
In [25]:
# |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
In [26]:
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

Other measures

  • Was there really anything special about $\ket{0}$ and $\ket{1}$?
  • No! :D
    • Except they pointed in opposite directions.
  • Can we make other measurements?
    • Yes!
In [27]:
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
In [28]:
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]

Let's simulate the Stern-Gerlach with QuTiP!

Simulating a simple classical quantum system.

In [29]:
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
In [30]:
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
In [31]:
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
In [32]:
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")
In [33]:
atoms, spins = quantum_stern_gerlach(1000)
plot_quantum_results(atoms, spins)