import numpy as np
import matplotlib.pyplot as plt
%config InlineBackend.figure_formats = ["svg"]
Spenser Talkington
Assignment 1
Physics 662: Computational Many-Body Physics
26 January 2023
I found this resource helpful for the Jordan-Wigner transformation: (https://learn.microsoft.com/en-us/azure/quantum/user-guide/libraries/chemistry/concepts/jordan-wigner)
s0 = np.array([[1,0],[0,1]],dtype=float)
s1 = np.array([[0,1],[1,0]],dtype=float)
s2 = np.array([[0,-1j],[1j,0]],dtype=complex)
s3 = np.array([[1,0],[0,-1]],dtype=float)
sp = np.array([[0,1],[0,0]],dtype=float)
sm = np.array([[0,0],[1,0]],dtype=float)
def cD_(n, L):
if(n==L-1):
result = sm
else:
result = s0
for i in range(L-2,n,-1):
result = np.kron(s0,result)
result = np.kron(sm,result)
for i in range(n,0,-1):
result = np.kron(s3,result)
return result.astype(float)
def c_(n, L):
if(n==L-1):
result = sp
else:
result = s0
for i in range(L-2,n,-1):
result = np.kron(s0,result)
result = np.kron(sp,result)
for i in range(n,0,-1):
result = np.kron(s3,result)
return result.astype(float)
def vacuum(L):
vac = np.zeros(2**L)
vac[0] = 1
return vac
Indexing scheme:
$m\sigma=\{d_{xz}\uparrow,d_{yz}\uparrow,d_{xy}\uparrow,d_{x^2-y^2}\uparrow,d_{z^2}\uparrow,d_{xz}\downarrow,d_{yz}\downarrow,d_{xy}\downarrow,d_{x^2-y^2}\downarrow,d_{z^2}\downarrow\}$ $m\sigma=\{0\quad\ \ ,1\quad\ \ ,2\quad\ \ ,3\qquad\ \ ,4\quad\ \ ,5\quad\ \ ,6\quad\ \ ,7\quad\ \ ,8\qquad\ \ ,9\quad\ \ \}$
def CrystalFieldHamiltonian(Dq):
H_cf = np.zeros((2**10,2**10),dtype=float)
for s in [0,5]:
H_cf += 6*Dq *(cD_(s+3,10) @ c_(s+3,10) + cD_(s+4,10) @ c_(s+4,10)) - 4*Dq*(cD_(s+2,10) @ c_(s+2,10) + cD_(s+0,10) @ c_(s+0,10) + cD_(s+1,10) @ c_(s+1,10))
return H_cf
We express the Coulomb interation term in terms of Slater-Condon parameters. Indexing scheme:
$m=\{d_{xz},d_{yz},d_{xy},d_{x^2-y^2},d_{z^2}\}$
$m=\{0\ \ \ \, ,1\ \ \ ,2\ \ \ \, ,3\quad \ \ \ ,4\ \ \ \}$
def SlaterU(F0, F2, F4):
U = np.zeros((5,5,5,5))
#direct Coulomb
for m in range(0,5,1):
U[m,m,m,m] = F0 + 4*F2 + 36*F4
#density repulsion
U[3,0,0,3] = U[0,3,3,0] = U[3,1,1,3] = U[1,3,3,1] = U[2,0,0,2] = U[0,2,2,0] = U[2,1,1,2] = U[1,2,2,1] = U[0,1,1,0] = U[1,0,0,1] = F0 - 2*F2 - 4*F4
U[4,0,0,4] = U[0,4,4,0] = U[4,1,1,4] = U[1,4,4,1] = F0 + 2*F2 - 24*F4
U[4,2,2,4] = U[2,4,4,2] = U[4,3,3,4] = U[3,4,4,3] = F0 - 4*F2 + 6*F4
U[3,2,2,3] = U[2,3,3,2] = F0 + 4*F2 - 34*F4
#exchange coefficients
U[2,0,2,0] = U[2,2,0,0] = U[0,0,2,2] = U[0,2,0,2] = U[2,1,2,1] = U[2,2,1,1] = U[1,1,2,2] = U[1,2,1,2] = U[0,1,0,1] = U[0,0,1,1] = U[1,1,0,0] = U[1,0,1,0] = U[3,0,3,0] = U[3,3,0,0] = U[0,0,3,3] = U[0,3,0,3] = U[3,1,3,1] = U[3,3,1,1] = U[1,1,3,3] = U[1,3,1,3] = 3*F2 + 20*F4
U[4,3,4,3] = U[4,4,3,3] = U[3,3,4,4] = U[3,4,3,4] = U[4,2,4,2] = U[4,4,2,2] = U[2,2,4,4] = U[2,4,2,4] = 4*F2 + 15*F4
U[4,0,4,0] = U[4,4,0,0] = U[0,0,4,4] = U[0,4,0,4] = U[4,1,4,1] = U[4,4,1,1] = U[1,1,4,4] = U[1,4,1,4] = F2 + 30*F4
U[3,2,3,2] = U[3,3,2,2] = U[2,2,3,3] = U[2,3,2,3] = 35*F4
#multi-orbital interactions
init = np.array([[0,4,3,0],[1,4,3,1],[0,0,3,4],[1,1,3,4],[4,2,1,0],[4,2,0,1],[4,0,1,2],[3,2,1,0],[3,2,0,1]],dtype=int)
vals = np.array([-2*np.sqrt(3)*F2 + 10*np.sqrt(3)*F4,
2*np.sqrt(3)*F2 - 10*np.sqrt(3)*F4,
np.sqrt(3)*F2 - 5*np.sqrt(3)*F4,
-np.sqrt(3)*F2 + 5*np.sqrt(3)*F4,
np.sqrt(3)*F2 - 5*np.sqrt(3)*F4,
np.sqrt(3)*F2 - 5*np.sqrt(3)*F4,
-2*np.sqrt(3)*F2 + 10*np.sqrt(3)*F4,
3*F2 - 15*F4,
-3*F2 + 15*F4])
for i in range(0,9,1):
i1,i2,i3,i4 = init[i]
val = vals[i]
U[i1,i2,i3,i4] = U[i1,i3,i2,i4] = U[i4,i2,i3,i1] = U[i4,i3,i2,i1] = U[i2,i1,i4,i3] = U[i2,i4,i1,i3] = U[i3,i1,i4,i2] = U[i3,i4,i1,i2] = val
return U
def SlaterHamiltonian(F0, F2, F4):
H_int = np.zeros((2**10,2**10),dtype=float)
U = SlaterU(F0,F2,F4)
for s1 in [0,5]:
for s2 in [0,5]:
for m1 in range(0,5,1):
for m2 in range(0,5,1):
if(m2!=m1 or s2!=s1):
for m3 in range(0,5,1):
for m4 in range(0,5,1):
if(m3!=m4 or s2!=s1):
H_int += (1/2) * U[m1,m2,m3,m4] * cD_(m1+s1,10) @ cD_(m2+s2,10) @ c_(m3+s2,10) @ c_(m4+s1,10)
return H_int
Rotated into the fixed particle number subspace to save on the computational complexity.
def indicesFixedN(N, L):
#this sequence is given by: https://oeis.org/A000120
tot = np.zeros((2**L,2**L))
for i in range(0,L,1):
tot += cD_(i,L) @ c_(i,L)
#print(tot)
my_list = []
for i in range(0,2**L,1):
if(tot[i,i]==N):
my_list.append(i)
return np.array(my_list)
def dn_Hamiltonian(F0, F2, F4, Dq, n):
H_total = CrystalFieldHamiltonian(Dq) + SlaterHamiltonian(F0, F2, F4)
#rotate to diagonal block in upper left
indices = indicesFixedN(n,10)
num = np.shape(indices)[0]
rot = np.zeros((2**10,2**10),dtype=float)
for i in range(0,num,1):
rot[i,indices[i]] = 1
H_reduced = rot @ H_total @ rot.T
#H_reduced = H_total[indices,:][:,indices]
return H_reduced[0:num,0:num]
For
%%time
for n in [2,3,4]:
H_cf = CrystalFieldHamiltonian(1) #Dq=1, but it is proportional to Dq
H_int = SlaterHamiltonian(0, 1/(1-5*0.078), 0.078/(1-5*0.078)) #units where B=1
#rotate to diagonal block in upper left
indices = indicesFixedN(n,10)
num = np.shape(indices)[0]
rot = np.zeros((2**10,2**10),dtype=float)
for i in range(0,num,1):
rot[i,indices[i]] = 1
H_cf_r = (rot @ H_cf @ rot.T)[0:num,0:num]
H_int_r = (rot @ H_int @ rot.T)[0:num,0:num]
#calculate many-body energy spectrum
N = 100
energies = np.zeros((N,int(np.math.factorial(10)/(np.math.factorial(n)*np.math.factorial(10-n)))))
for i in range(0,N,1):
energies[i] = np.linalg.eigvalsh(H_int_r+H_cf_r*i*4/N)
#plot
Dqs = np.linspace(0,4,np.shape(energies)[0])
for i in range(0,int(np.math.factorial(10)/(np.math.factorial(n)*np.math.factorial(10-n))),1):
plt.plot(Dqs,energies[:,i]-energies[:,0],c="k")
plt.ylim(-1,40)
plt.xlabel("Dq/B")
plt.ylabel("E/B")
plt.title(r"Tanabe-Sugano Diagram for $d^{a}$ Transition Metals".format(a=n))
plt.show()
plt.close()
CPU times: user 49min 5s, sys: 9min 27s, total: 58min 32s Wall time: 7min 26s