{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# HGSFP Graduate Days SS2019\n",
"\n",
"## Programming exercise 1: Spin model toolbox and Ising model"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# load standard libraries\n",
"\n",
"import numpy as np # standard numerics library\n",
"import numpy.linalg as LA\n",
"\n",
"import matplotlib.pyplot as plt # for making plots\n",
"\n",
"import time as time\n",
"\n",
"import scipy.sparse as sparse\n",
"import scipy.sparse.linalg as sLA\n",
"\n",
"from numpy import (array, pi, cos, sin, ones, size, sqrt, real, mod, append, arange)\n",
"\n",
"#from qutip import (tensor, basis) \n",
"\n",
"%matplotlib inline"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Building general spin 1/2 Hamiltonians\n",
"\n",
"We want to consider a model with nearest neighbor spin interactions, the one-dimensional transverse-field Ising model (TFIM). This model is analytically solvable (see e.g. https://www.sciencedirect.com/science/article/pii/0003491670902708, note the slightly different definition of the model parameters in this reference. Their $\\Gamma$ is $2B$ and their $J$ is $4J$ in our convention). For the sake of learning how to build general spin models, we will solve this model here by numerical means and compare the result to the exact solution as a check. The tools you develop here, you can then apply to build any spin Hamiltonian straight forwardly. The Hamiltonian of the TFIM reads\n",
"$$\n",
"H=\\sum_{i=0}^{N-1} -J\\sigma_z^{(i)}\\sigma_z^{(i+1)} - B \\sigma_x^{(i)}\n",
"$$\n",
"where we want to use periodic boundary conditions, i.e. the Nth spin is identified with the 0th spin.\n",
"\n",
"To solve it we first set up a framework for general spin 1/2 models.\n",
"We will work in the canonical product basis of states $|i_1,...,i_N\\rangle$, where $i_k \\in {0,1}$. This basis of course has the problem of exponential scaling of the number of basis states with particle number N. But we still want to use it to keep our approach as general as possible, which then just allows us to use up to 12 or so spins with reasonable computation time.\n",
"\n",
"The recipe for conctructing general spin models is the following:\n",
"\n",
"1) Build the single-spin operators $\\sigma_x$, $\\sigma_y$, $\\sigma_z$, and the 2x2 identity, $\\mathbb{1}$ as numpy arrays. You can also work with sparse matrices since the Hamiltonian operator will contain a lot of zeros this should significantly speed up your code (csr_matrix should work best).\n",
"\n",
"2) Build a list of all single-spin operators acting on spin $i=0...N-1$ in the $N$-spin Hilbert space.\n",
"$$\n",
"\\sigma_\\alpha^{(i)} = \\mathbb{1}^{\\otimes i}\\otimes \\sigma_\\alpha \\otimes \\mathbb{1}^{\\otimes (N-i-1)}\n",
"$$\n",
"where $\\alpha \\in \\{x,y,z\\}$. Use the kronecker product np.kron() (for sparse matrices sparse.kron()) to do this. (Make sure you understand how the indexing/ordering of the matrix elements is done by kron()!)\n",
"\n",
"3) Build the Hamiltonian of the transverse field Ising model by adding up its individual terms, using the dot product for interaction terms.\n",
"\n",
"Test your implementation by calculating the ground state energy of the transverse field Ising model and comparing it to the analytical result (for even $N$)\n",
"$$\n",
"E_0^{analyt} = -\\sum_{k=-(N-1)/2}^{(N-1)/2}\\sqrt{1+B^2 + 2B\\cos(2\\pi k/N)}\n",
"$$\n",
"where the convention $J=1$ was chosen.\n",
"For N=10, scan $B$ from 0 to 2 and plot the numerical and analytical result for the ground state energy.\n",
"\n",
"Also think about what the ground state should be at large $B$ and at $B=0$. \n",
"\n",
"Thoughts on diagonalization routines:\n",
"- If you work with full matrices, you can use LA.eig() which will give you the eigenvalues and eigenfunctions. Here the eigenvalues will not be sorted by their magnitude, so the the ground state is not the first eigenvlaue.\n",
"- Since we are dealing with Hermitian matrices, which have real eigenvalues, it makes sense to use eigh() instead, which returns the eigenvalues already sorted.\n",
"- If you use sparse matrices you should use scipy.sparse.linalg.eigsh(). If you are only interested in the ground state or low lying states, this can be very efficient. You have to specify how many states you want to calculate where in the spectrum (in our case 'SA').\n",
"\n",
"If you are already done with the above...\n",
"\n",
"4) Now you can calculate the gap between the ground and excited states as a function of the field strength to see the closing of the gap at the phase transition point. Note that the model also has the $Z_2$ symmetry (invariance under flipping of all spins) so you would look at the gap of ground and second excited state.\n",
"What other properties could you use to see the quantum phase transition?\n",
"\n",
"5) Dynamics: If you calculate all eigenenergies and eigenfunction you can now also calculate the time evolution of the state after a quench. E.g. choose an initial state where all spins are in the state $(|0\\rangle+|1\\rangle)/\\sqrt{2}$ (ground state for infinite B) and evolve it in time under the Hamiltonian with some value of J and B to see how the magnetization changes in time."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.7.0"
}
},
"nbformat": 4,
"nbformat_minor": 2
}