Having seen some of the typical maths you need in part I and part II we’ll turn now to some Python showing how to simulate quantum systems.

I’ll use below the QuTip Python library but there are other very nice libs out there. Also, if Python is not your cup of tea you can also take a look at this list of QC simulators here, the Liquid F# lib in particular.


Creating qubits, bras and kets goes typically likes this:

from qutip import *
import numpy as np
sqrt = np.sqrt
import matplotlib.pyplot as plt
plot = plt.plot
%matplotlib inline

which returns the $|0\rangle$ state in a 4-dimensional space:

$\begin{pmatrix}1.0 \\ 0.0\\0.0\\0.0 \end{pmatrix}$

The tangled EPR state is hence

epr = (basis(4,0) + basis(4,3))/sqrt(2)

$\begin{pmatrix}0.707 \\ 0.0\\0.0\\0.707 \end{pmatrix}$

We saw earlier that qubits can be represented on a Bloch sphere and the library gives you easy access to visualizations:

bloch = Bloch()
state = sigmaz()*sigmax()*basis(2,0)


The sigmax, sigmay, sigmaz above are the Pauli spin matrices. Turning the states into a density matrix can be done with the ket2dm command:


$\begin{pmatrix}0.5 & 0.0 & 0.0 & 0.5\\ 0.0 & 0.0 & 0.0 & 0.0\\0.0 & 0.0 & 0.0 & 0.0\\0.5 & 0.0 & 0.0 & 0.5 \end{pmatrix}$

Remember that this is nothing but turning an arbitrary ket $|\psi\rangle$ into a matrix $|\psi\rangle\langle\psi|$.

Tensor products are obtained using the tensor method. For example, creating the tensor product of $\sigma_x\otimes\sigma_y$:


$\begin{pmatrix}0.0 & 0.0 & 0.0 & -1.0j\\ 0.0 & 0.0 & 1.0j & 0.0\\0.0 & -1.0j & 0.0 & 0.0\\1.0j & 0.0 & 0.0 & 0.0 \end{pmatrix}$

so, $(\sigma_x\otimes\sigma_y)(|1\rangle\otimes|0\rangle)$ is hence

t = tensor(sigmax(),sigmay())
t*tensor(basis(2,1), basis(2,0))

$\begin{pmatrix}0.0 \\ 1.0j\\0.0\\0.0 \end{pmatrix}$.


Operators are created via the Qobj

op = Qobj([[0,5j], [-5j,0]])

and applied via multiplication on a state

op*(basis(2,0) + basis(2,1))

$\begin{pmatrix}5.0j \\ -5.0j \end{pmatrix}$.

Eigenstates and eigenvalues are fetched via the eigenstates method on the operator

eig = sy.eigenstates()
print("The eingenvalues %s has eigenvector"%eig[0][0])

$\begin{pmatrix}-0.707j \\ -0.707j \end{pmatrix}$.

Fock space mechanics works as should be and you can check e.g. that $a^\dagger|0\rangle = |1\rangle$

create(2)*basis(2,0) == basis(2,1)

Similarly, $a^\dagger a^\dagger|0\rangle = 0$ since ${a^\dagger, a^\dagger}=0$:


$\begin{pmatrix}0.0 \\ 0.0 \end{pmatrix}$.


The framework allows you to simulate the evolution of kets or observables using the Lindblad or master equation. For example, the following simple Hamiltonian

$\hat{H} = \sigma_x\sigma_z$

evolution on the $|0\rangle$ ket gives:

H = sigmax()*sigmaz()
psi0 = basis(2, 0)
# list of times for which the solver should store the state vector
tlist = np.linspace(0, 10, 100)
result = mesolve(H, psi0, tlist, [], [])
plot([np.absolute(x[0][0]) for x in result.states])