{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# HGSFP Graduate Days SS2019\n",
"\n",
"## Programming exercise 3: Spin model toolbox and Ising model"
]
},
{
"cell_type": "code",
"execution_count": 2,
"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": [
"### Reduced density matrix and entanglement entropy\n",
"\n",
"The TFIM has a quantum phase transition at $B/J=1$ from a ferromagnetic to a paramagnetic phase (see exercise 1). We saw in exercise 1 that the gap is mimimal around this value. An interesting feature is also that at the critical point the half-chain entanglement entropy diverges logarithmically as $N$ goes to infinity (see for example https://arxiv.org/abs/0803.3610 Fig. 3). For finite $N$ we should already see that the entanglement entropy becomes maximal approximately at the critical point.\n",
"\n",
"To calculate the half-chain entanglement entropy in the ground state of the transverse Ising chain, you should proceed as follows:\n",
"\n",
"1) Calculate the partial trace over $N/2$ spins. The matrix elements of the reduced density of a state $|\\psi\\rangle = \\sum_i c_i |i\\rangle$ (where $|i\\rangle$ are the $N$-spin basis states) can be expressed as\n",
"$$\n",
"(\\rho_{red,N/2})_{ij} = \\sum_{k=0}^{2^{N/2}-1} c_{i+k*2^{N/2}}c^*_{j+k*2^{N/2}}\n",
"$$\n",
"(This might look different depending on how the states are ordered in your basis. The form given here corresponds to adding up the $2^{N/2}\\times 2^{N/2}$ diagonal blocks, as sketched in the lecture notes.)\n",
"\n",
"2) Calculate the eigenvalues $p_i$ of $\\rho_{red,N/2}$, sort them from largest to smallest. This is the entanglement spectrum. The rank of $\\rho_{red,N/2}$ is called entanglement dimension (or Schmidt rank). For a separable state the entanglement dimension is 1, i.e. only one eigenvalue is non-zero.\n",
"\n",
"3) Calculate the von-Neumann entanglement entropy $S_E = -\\sum_i p_i \\log(p_i)$. Consider that some of the eigenvalues can be zero. Due to the finite numerical precision they can even turn negative. These you want to exclude from the sum. (The limit $\\lim_{p\\rightarrow 0} p\\log(p) =0$ so we are not actually modifying the outcome significantly by neglecting small p's.)\n",
"\n",
"Test your code for some cases where you know the answer, for example a separable state or a Bell pair state of 2 qubits. The qutip module also has the capability of calculating partial traces. If you are interested you can compare to this for testing.\n",
"\n",
"Then calculate the entanglement spectrum and entropy for the same parameter scan as in exercise 1. Interpret your observations.\n",
"\n",
"If you still have time:\n",
"\n",
"Use singular value decomposition to calculate the entanglement spectrum: Write the coefficients of your bipartite state $c_{ij}$ ($i$ labels the basis of subsystem A and j the basis states of subsystem B) as a matrix using the reshape command. Then calculate the singular values of this using svd() from numpy.linalg, which correspond to the Schmidt coefficients, i.e. the square roots of the p's from above. Is this numerically more efficient than calculating the reduced density matrix and diagonalizing it?\n",
"\n",
"Calculate the time-evolution for all spins being in the state $|+\\rangle = (|0\\rangle+|1\\rangle)/\\sqrt{2}$ initially.\n",
"If you evolve this state with the critical Ising Hamiltonian, the enatnglement entropy should increase very rapidly in time and saturate due to the finite system size. The saturation value is extensive in system size (volume law). This makes simulating the time evolution for quenches near the critical point challenging for techniques based on matrix prduct states. "
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [],
"source": [
"# useful definitions for building the TFIM Hamiltonian\n",
"\n",
"# spin operators\n",
"sx = np.array([[0,1],[1,0]])\n",
"sy = np.array([[0,1j],[-1j,0]])\n",
"sz = np.array([[1,0],[0,-1]])\n",
"\n",
"# sparse spin operators\n",
"sxSp = sparse.csr_matrix(sx)\n",
"sySp = sparse.csr_matrix(sy)\n",
"szSp = sparse.csr_matrix(sz)\n",
"oneSp = sparse.identity(2)\n",
"\n",
"# function that builds all single spin operators\n",
"def BuildSingleSpinOps(N):\n",
" sxi=[]\n",
" syi=[]\n",
" szi=[]\n",
" # build the single spin operators\n",
" for i in range(N):\n",
" sxi.append( sparse.kron(sparse.kron(sparse.identity(2**(N-i-1)),sxSp),sparse.identity(2**i)) )\n",
" syi.append( sparse.kron(sparse.kron(sparse.identity(2**(N-i-1)),sySp),sparse.identity(2**i)) )\n",
" szi.append( sparse.kron(sparse.kron(sparse.identity(2**(N-i-1)),szSp),sparse.identity(2**i)) )\n",
" \n",
" return sxi, syi, szi\n",
"\n",
"# Build the transverse field ising Hamiltonian with periodic boundaries\n",
"def Build_H_TFIM(N,B):\n",
" sxis, syis, szis = BuildSingleSpinOps(N)\n",
" dim=2**N\n",
" Hmat = sparse.csr_matrix((dim,dim))\n",
" for i in range(N):\n",
" Hmat = Hmat - szis[i] @ szis[(i+1)%N] # interaction term\n",
" Hmat = Hmat - B*sxis[i] # field term\n",
" return Hmat\n",
"\n"
]
},
{
"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
}