- EuroSciPy 2019
- Simon Cross

- 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.

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

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()
```

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}$

Why this weird quantum mechanics anyway?

See also:

- https://plato.stanford.edu/entries/physics-experiment/app5.html.
- https://en.wikipedia.org/wiki/Stern%E2%80%93Gerlach_experiment#History

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.

\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}

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.

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()
```

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)
```

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)
```

This really is a system with two outcomes.

- 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.

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

Something we're familiar with.

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

- States:
`0`

and`1`

- Outcomes:
`0`

and`1`

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

- We can label these two states
`0`

and`1`

. - Measuring state
`0`

produces an outcome which we will also label`0`

. - Measuring state
`1`

produces an outcome which we will also label`1`

. - These are the only two states.
- Measuring the same state always produces the same outcome.

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:

- A coin can lie on a surface in two different ways. We measure it by looking at it. We call the outcomes
`heads`

and`tails`

. - A digital signal can have one of two different voltages. We measure it by measuring the voltage. We call the outcomes
`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.

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))
```

Putting the puzzle pieces together.

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".

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

- We can label these two basis states
`|0>`

and`|1>`

. - Measuring state
`|0>`

produces an outcome which we will label`0`

. - Measuring state
`|1>`

produces an outcome which we will label`1`

. - There are other states called
*superpositions*which we label`a|0> + b|1>`

.`a`

and`b`

are complex numbers.

- Measuring a state produces outcome
`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:

- A coin can lie on a surface in two different ways. We measure it by looking at it. We call the outcomes
`heads`

and`tails`

. - A digital signal can have one of two different volates. We measure it by measuring the voltage. We call the outcomes
`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.

In [13]:

```
b0 = ket("0") # |0>
b1 = ket("1") # |1>
print("State:\n", b1)
```

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))
```

The tricky part. :)

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()
```

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$

- $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)`

- Length:

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

- 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.

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).

- 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$

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!)

Wow, this looks a lot like a direction in space! Making progress!

- 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.

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"],
])
```

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))
```

In [26]:

```
qbit = (1 * ket("0")) + (1j * ket("1"))
qbit = qbit.unit()
print("State:\n", qbit)
print("Outcome:", measure_qbit(qbit))
```

- 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)])
```

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)
```