# ID Quantum Ising Model¶

Today we are going to use concepts from quantum computation to solve a real physics problem. We want to solve the one dimensional quantum Ising model.

$H = -J\sum_{i=1}^{N}{\sigma_i^x \sigma_{i+1}^x} - \mu_B B \sum_{i=1}^N{\sigma_i^z}$

The Ising model is a simplified model of how a magnet works. It is useful to study because it is a good example of the physics and computations we need to consider.

First we have to understand what the situation is. We have some material where the electrons are located along a one dimensional chain. We call the location of the electrons "lattice sites". An electron can be spinning up $|\uparrow>$ or spinning down $|\downarrow>$. Also, we have a magnetic field. We can choose any direction, but by convention we say the magnetic field points in the z direction. Then a spin up is aligned with the magnetic field and a spin down is anti-aligned with the magnetic field.

Now let us understand the terms in the model:

1. H is the Hamiltonian. For closed systems this is also the energy
2. J is a number, the units are energy and it represents the interaction energy between two electrons
3. $\mu_B B$ is the energy due to an electron being aligned with the magnetic field
4. The summation runs from 1 to N. N is the number of electrons
5. The electron interaction term reads $\sigma_i \sigma_{i+1}$, this says that an electron at position 1 only interacts with the electron at position 2. But, there are more electrons, we just consider neighbors.
6. We use something called periodic boundary conditions. This is like the game asteroid: when you hit the right edge of the screen the spaceship shows back up on the left side of the screen. For example, if $N=3$ then we will have 3 electron interaction terms: $\sigma_1\sigma_2$, $\sigma_2\sigma_3$ $\sigma_3\sigma_1$. Note that when we reached $\sigma_3$ the next one was $\sigma_1$. These are periodic boundary conditions.

This model is interesting because it is a good example of complications when there are two competing terms in the energy. If $J >> B$ then an electron will want to align with it's neighbor. However, if $J << B$ then the electron will want to align with the magnetic field. Now what if $J$ and $B$ are comparable in magnitude? Well, then that is the heart of the many-body problem. There is no simple solution because the behavior depends on a lot of variables. So what if you just have a few electrons and you want to compute the ground state and ground state energy?

You could do a statistical mechanics treatment of the problem; find the partition function and then compute the energies. If you do this you may quit physics for the rest of your life. I want to show you how to solve this problem in a very clean way using concept from quantum computing. First, let me just change the notation a bit and let $h = \mu_B B$ then the hamiltonian becomes

$H = -J\sum_{i=1}^{N}{\sigma_i^x \sigma_{i+1}^x} - h \sum_{i=1}^N{\sigma_i^z}$

One electron can be described by $|\uparrow>$ or $|\downarrow>$, since these labels are arbitrary we can use $|\uparrow> = |0>$ and $|\downarrow> = |1>$ where $|0>$ and $|1>$ are qubits. Just a reminder

$|0> = \left(\begin{array}{r} 1 \\ 0 \end{array}\right) \qquad |1> = \left(\begin{array}{r} 0 \\ 1 \end{array}\right)$

$\sigma^z$ is the Pauli-Z matrix $Z$ and $\sigma^x$ is the Pauli-X matrix $X$. As a reminder:

$X = \left(\begin{array}{rr} 0 & 1 \\ 1 & 0 \end{array}\right) \qquad Z = \left(\begin{array}{rr} 1 & 0 \\ 0&-1 \end{array}\right)$

This is all for one electron. The Ising model is not very interesting for one electron, so how do we deal with N electrons? The answer is a postulate of quantum mechanics that is overlooked in most intro classes:

### Postulate

The state space of a composite quantum system is the tensor product of the individual systems. For example, suppose we have a two electron system. The first electron is spin up $|0>$ and the second electron is sping down $|1>$. Then then if we bring these two electrons together and let them interact, the composite system will be $|0> \bigotimes |1>$ where $\bigotimes$ is the tensor product. Also note that

$|0> \otimes |1> = |0>|1> = |01> = \left(\begin{array}{r} 0\\1\\0\\0 \end{array}\right)$

Ok, so we have a solution for interacting electrons: just tensor their states. But in the Ising model only neighbors interact. So how do we represent an non interaction between electrons? The answer is to do a fancy version of multiplying by 1. So, if electrons do not interact we tensor with the identity matrix I:

$I = \left(\begin{array}{rr} 1 & 0 \\ 0 & 1 \end{array}\right)$

We are now ready to do a complete example. Suppose we have 3 electrons: N=3 then the Ising model becomes

\begin{equation} \begin{aligned} H &= -J\sum_{i=1}^{3}{\sigma_i^x \sigma_{i+1}^x} - h \sum_{i=1}^3{\sigma_i^z} \\ &= -J\sum_{i=1}^{3}{X_i X_{i+1}} - h \sum_{i=1}^3{Z_i} \\ &= -J[(X_1\otimes X_2 \otimes I) + (I \otimes X_2 \otimes X_3) + (X_1 \otimes I \otimes X_3) ] -h[(Z_1 \otimes I \otimes I) + (I \otimes Z_2 \otimes I) + (I \otimes I \otimes Z_3)] \end{aligned} \end{equation}

Below we use the quantum computation library QuDotPy https://github.com/psakkaris/qudotpy to make our lives easy. We use it to calculate the ground state energy and the ground state. We also verify the ground state energy by explicit calculation of the expectation value $<g|H|g>$ where $|g>$ is the ground state.

In :
from qudotpy import qudot

interaction_mat = (qudot.tensor_gate_list([qudot.X, qudot.X, qudot.I]) +
qudot.tensor_gate_list([qudot.I, qudot.X, qudot.X]) +
qudot.tensor_gate_list([qudot.X, qudot.I, qudot.X]))

z_mat = (qudot.tensor_gate_list([qudot.Z, qudot.I, qudot.I]) +
qudot.tensor_gate_list([qudot.I, qudot.Z, qudot.I]) +
qudot.tensor_gate_list([qudot.I, qudot.I, qudot.Z]))

# for simplicity set J = h = 1 Joule

H = -1 * interaction_mat - z_mat
print H

[[-3.+0.j  0.+0.j  0.+0.j -1.+0.j  0.+0.j -1.+0.j -1.+0.j  0.+0.j]
[ 0.+0.j -1.+0.j -1.+0.j  0.+0.j -1.+0.j  0.+0.j  0.+0.j -1.+0.j]
[ 0.+0.j -1.+0.j -1.+0.j  0.+0.j -1.+0.j  0.+0.j  0.+0.j -1.+0.j]
[-1.+0.j  0.+0.j  0.+0.j  1.+0.j  0.+0.j -1.+0.j -1.+0.j  0.+0.j]
[ 0.+0.j -1.+0.j -1.+0.j  0.+0.j -1.+0.j  0.+0.j  0.+0.j -1.+0.j]
[-1.+0.j  0.+0.j  0.+0.j -1.+0.j  0.+0.j  1.+0.j -1.+0.j  0.+0.j]
[-1.+0.j  0.+0.j  0.+0.j -1.+0.j  0.+0.j -1.+0.j  1.+0.j  0.+0.j]
[ 0.+0.j -1.+0.j -1.+0.j  0.+0.j -1.+0.j  0.+0.j  0.+0.j  3.+0.j]]


In :
import numpy as np

w, v = np.linalg.eigh(H)
ground_state_energy = w
ground_state = qudot.QuState.init_from_vector(np.asarray(v[:,0]))

#verify
print "ground state energy: %s" % ground_state_energy
print "ground state: " + str(ground_state)
energy_check = (ground_state.bra * H * ground_state.ket)[0,0].real
print energy_check
ground_state.possible_measurements()

ground state energy: -4.0
ground state: -0.866025|000> + -0.288675|011> + -0.288675|101> + -0.288675|110>
-3.99999985641


Out:
{u'|000>': 0.75,
u'|011>': 0.083333328,
u'|101>': 0.083333328,
u'|110>': 0.083333328}


Ok, so it is very probable that spins will align with the field but we can also get nearest neighbor alignment. Remember that $|0>$ means aligned with magnetic field and $|1>$ is anti aligned. So just for fun let's crank up the magnetic field. We expect for all spins to align with the field with very high probability.

In :
H = -1 * interaction_mat - 1000 * z_mat
w, v = np.linalg.eigh(H)
ground_state_energy = w
ground_state = qudot.QuState.init_from_vector(np.asarray(v[:,0]))

print "ground state energy: %s" % ground_state_energy
print "ground state: " + str(ground_state)
ground_state.possible_measurements()

ground state energy: -3000.00075038
ground state: -1.0|000> + -2.13412e-23|010> + -0.000250125|011> + -2.13412e-23|100> + -0.000250125|101> + -0.000250125|110> + -7.11609e-24|111>


Out:
{u'|000>': 0.99999976,
u'|010>': 0.0,
u'|011>': 6.2562506e-08,
u'|100>': 0.0,
u'|101>': 6.2562506e-08,
u'|110>': 6.2562506e-08,
u'|111>': 0.0}


Aha! so all spins now want to align with the magnetic field with a probability very close to 1, just like we expect. Note how low the other probabilities are! If you are interested I encourage you to play with J and h to see what happens. Some results may surprise you! Also note that very low probabilities are usually an artifact of floating point leftovers.

### Hamiltonians Extension to QuDotPy

To make things even easier we updated QuDotPy to support building Ising hamiltonians and solving for the ground state. This way you do not need to tensor gates all the time (can get annoying for large N). Also, you will not need to know the details of numpy to solve for the ground state.

In :
from qudotpy import hamiltonians
ising3 = hamiltonians.create_ising_ham(J=1, h=1, N=3)
print ising3

[[-3.+0.j  0.+0.j  0.+0.j -1.+0.j  0.+0.j -1.+0.j -1.+0.j  0.+0.j]
[ 0.+0.j -1.+0.j -1.+0.j  0.+0.j -1.+0.j  0.+0.j  0.+0.j -1.+0.j]
[ 0.+0.j -1.+0.j -1.+0.j  0.+0.j -1.+0.j  0.+0.j  0.+0.j -1.+0.j]
[-1.+0.j  0.+0.j  0.+0.j  1.+0.j  0.+0.j -1.+0.j -1.+0.j  0.+0.j]
[ 0.+0.j -1.+0.j -1.+0.j  0.+0.j -1.+0.j  0.+0.j  0.+0.j -1.+0.j]
[-1.+0.j  0.+0.j  0.+0.j -1.+0.j  0.+0.j  1.+0.j -1.+0.j  0.+0.j]
[-1.+0.j  0.+0.j  0.+0.j -1.+0.j  0.+0.j -1.+0.j  1.+0.j  0.+0.j]
[ 0.+0.j -1.+0.j -1.+0.j  0.+0.j -1.+0.j  0.+0.j  0.+0.j  3.+0.j]]


In :
(energy, state) = hamiltonians.solve_ground_state(ising3)
print "energy: %s \n state: %s" % (energy, state)

energy: -4.0
state: -0.866025|000> + -0.288675|011> + -0.288675|101> + -0.288675|110>