Put the module file in the same directory as your working script of notebook files.
# import for Python 3
from quantumising_sparse import *
# import for Python 2
from quantumising_sparse_p2 import *
Let's write the Hamiltonian of the Ising in a transverse field (ITF) model on a chain of length $L$ $$ H = - J \sum_{j=1}^L S_j^zS_{j+1}^z -h\sum_{j=1}^L S_j^x$$ with periodic boundary conditions $S_{L+1} = S_1$.
1/ In the $J = 0$ limit, recall the ground-state $|\text{GS}\rangle$. The first excitation states correspond to flipping a spin at position $j$ along the $x$-axis $|\text{j}\rangle =|\text{+}\rangle_1|\text{+}\rangle_2\ldots|\text{-}\rangle_j\ldots|\text{+}\rangle_L $. Apply the Hamiltonian on this state. What is the energy of such a state ?
We introduce the following $q$-states using translational invariance (beware that $|0\rangle$ is not the ground-state): $$ \vert{q}\rangle = \frac{1}{\sqrt{L}}\sum_{j=1}^L e^{iqj}|\text{j}\rangle$$
2/ To first order perturbation theory in $J\ll 1$, compute the energies of the ground-state $E_0^{(1)}$ and the energies of the $q$-states $E_q^{(1)}$. Compute the dispersion relation $\omega(q) = E_q^{(1)}-E_0^{(1)}$.
3/ Infer that in first-order perturbation theory, the gap to the first excited state reads $$ \Delta = h -h^c\quad \text{ for } h \geq h^c$$ What is the corresponding prediction for the critical point $h^c$ ? This perturbative result is actually true from $h=h^c$ to $h=\infty$ in the thermodynamical limit. For $h\leq h_c$, we give that the gap reads $$ \Delta = 2(h^c-h)\quad \text{ for } h \leq h^c$$
4/ Compute numerically the gap $\Delta$ as a function of $h$ for $L=8,10,12$ spins. For $h$, take 16 values between 0 and $1.5J$. Plot it against the exact result above.
from quantumising_sparse import *
def IsingHam(L=4,h=0.0,J=1.0):
HamH = GeneralOp([ [ -h, Sx(i)] for i in range(L) ])
HamIsing = GeneralOp([ [ -J, SzSz(i,(i+1)%L)] for i in range(L) ])
return HamIsing + HamH
def Simulation(L,hs = [0.0], num=???,J=1):
hilbert = IsingHilbert(L=L)
res = []
for h in hs:
Ham = IsingHam(L,h,J)
res.append(Ham.Energies(hilbert,num=num,algorithm="Arnoldi"))
return res
hs = np.linspace(0.0,1.5,16)
...
...
5/ Compute numerically the local magnetizations $\langle{S_j^{x,z}}\rangle$ for a $L=10$ spin chain as a function of $h$ (take 21 points between 0 and 1). Does it behave the way you expect from the quantum phase transition lecture ? Why ? Using the np.gradient() function to compute the derivative, show that the transverse magnetization has a cusp close to the critical point.
6/ Compute the correlation function $C^{z}(r) = \langle{S_0^{z}S_{r}^{z}}\rangle$ for $r$ in $[0,L/2]$. How does it behave in each phase ? How to compute correctly the order parameter $\langle{S_j^{z}}\rangle$ from these correlations ? Check that it has the appropriate behavior with this approach.
def SzFromCorrelation(L,h):
hilbert = IsingHilbert(L=L)
Ham = IsingHam(L,h)
gs = Ham.Groundstate(hilbert,algorithm="Arnoldi")
return ????
hs = np.linspace(0.0,1.0,51)
mzs = []
for L in [8,10,12,14]: # calculations are long, possibly add L=16 if you have time
mzs.append([ SzFromCorrelation(L,h) for h in hs ])
7/ Show graphically that the finite size scaling of the order parameter obtained as in the previous question follows reasonably well a scaling form $$ m(h,L) \propto L^{-\beta/\nu} \mathcal{M}(L^{1/\nu}(h-h_c)) $$ with $\mathcal{M}$ a scaling function, by performing a collapse of the data using the critical exponents of the 2D classical Ising model.
8/ From a few points in $h$ from each side of the critical point, show that the correlation length $\xi$ - obtained from fitting the correlations $C^z$ as a decaying exponential $Ae^{-r/\xi}+B$ - satisfy tends to diverge at the critical point. Use the script for fitting curve below and remove the first and last points of the correlations.
from scipy.optimize import curve_fit
# pos contains positions of range(L/2+1), Cr the correlations data
fitfunc = lambda x, a, b, c: a*np.exp(-x/b)+ c
p, pcov = curve_fit(fitfunc,pos[1:-1],Cr[1:-1],p0 = [1.0,1.0,0.5])
# optimized results for a is in p[0], b in p[1] and c in p[2]
scatter(pos,Cr)
x = np.linspace(0,L/2,101)
semilogy(x,fitfunc(x,p[0],p[1],p[2]))
show()
# when c is large, maybe subtract it from the data...
We would like to use the exact diagonalization technique to simulate the dynamics of a quantum system. As for ground-state calculations, we would like to handle Hilbert spaces as large as possible.
1/ Write the evolution operator for small time step $dt$ as a basic series expansion. Numerically, we need a simple criteria to truncate the expansion. What is your guess for the criteria ?
2/ Show that one can obtain $|{\psi(t+dt)}\rangle$ from $|{\psi(t)}\rangle$ using sparse-matrix vector multiplication and only a few temporary vectors. Implement a function timestep(psi) which propagates the wave-vector from $t$ to $t+dt$. Complete the following code
def timestep(Ham,psi,dt):
v1 = deepcopy(psi)
...
res = deepcopy(psi)
for n in range(1,100):
...
...
...
if ??criteria?? < 1e-15:
break
return res
In this part, we use the conservation of total $S_{tot}^z$ but not translational symmetries.
3/ Start from a state with one spin up in the middle of an odd chain and all spin-down everywhere else. Evolve in time with the Heisenberg Hamiltonian. What is the approximate speed of propagation ?
Code sample to be completed :
from heisenberg import *
def HeisenbergHam(L=4):
HamZZ = GeneralOp([ [ 1.0, SzSz(i,(i+1)%L)] for i in range(L) ])
HamXY = GeneralOp([ [ 1.0, Spinflip(i,(i+1)%L)] for i in range(L) ])
return HamZZ + HamXY
def magnetization(psi):
return np.array([ ??? for r in range(L)])
def simulation(psi0,H,times):
psi = zeroWf(psi0.Hilbert,dtype=np.complex128)
psi[:] = psi0[:]
res = [ magnetization(psi) ]
for dt in np.diff(times):
psi = ???
res.append(???)
return res
L = 25
sz = ??? # do not use Sz, already defined as an operator
hilbert = HeisenbergHilbert(L=L,Sz=sz)
Ham = HeisenbergHam(L)
psi0 = zeroWf(hilbert)
initialconf = [0]*(L//2)+[1]+[0]*(L//2)
print(initialconf)
psi0.setvalue(Conf.conf2int(initialconf), 1.0)
print(psi0.readable())
Ham.createMatrix(hilbert,sparse=True)
times = np.linspace(0,10,101)
mx = simulation(psi0,Ham,times)
from pylab import *
matshow(mx, aspect=L/float(len(times)))
show()
4/ Compute the ground-state of an antiferromagnetic Heisenberg chain of length $L=16$ in the sector $S_{tot}^z=0$ with periodic boundary conditions. Using an operator, flip the central spin to obtain a new (normalized...) state which will be the initial state for the time evolution (notice that flipping the central spin change the total spin by a amount of 1). WARNING: debug using small chains at the begining
5/ Propagate this out-of-equilibrium initial state and measure the magnetization profile as a function of time. Interpret the results as two counter-propagating excitations carrying a spin-$1/2$, the spinons. Make sketches of the propagation using Néel-like configurations.
L, sz = 16, 0
hilbert = HeisenbergHilbert(L=L,Sz=sz)
hilbertP = HeisenbergHilbert(L=L,Sz=???)
Ham = HeisenbergHam(L)
gs = Ham.Groundstate(hilbert)
psi0 = ???? * gs
psi0.normalize()
# needs to clear the matrix and reconstruct it on the new hilbert space.
Ham.having_matrix = False
Ham.createMatrix(???,sparse=True)
times = np.linspace(0,6,61)
mx = simulation(psi,Ham,times)
matshow(mx, aspect=L/float(len(times)))
show()