In [1]:
import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import brentq as find_root
plt.rc("figure", figsize=(3,2))
%config InlineBackend.figure_formats = ["svg"]
Define Helper Functions¶
In [2]:
adj = lambda u: np.conj(u.T)
In [3]:
def fermi(z):
if(beta*z<-30):
return 1
elif(beta*z>30):
return 0
else:
return 1/(1+np.exp(beta*z))
Effective Hamiltonian and Step¶
In [4]:
def h_eff(k,ns):
"""
k = momentum point in [-pi,pi)
ns = dictionary of current orbital & spin resolved densities, spin-flip susceptibilities
"""
#hopping term
t = t1 + t2*np.exp(-1j*k)
#up-up
E_A_up = m + U*ns["n_A_dn"] #(the other spin densities because these electrons feel the effect of the _other_ spin species)
E_B_up = -m + U*ns["n_B_dn"]
h_uu = np.array([
[E_A_up+hz,t],
[adj(t),E_B_up+hz]
],dtype=complex)
#dn-dns
E_A_dn = m + U*ns["n_A_up"] #(the other spin densities because these electrons feel the effect of the _other_ spin species)
E_B_dn = -m + U*ns["n_B_up"]
h_dd = np.array([
[E_A_dn-hz,t],
[adj(t),E_B_dn-hz]
],dtype=complex)
#up-dn
S_A = -U*ns["s_A"] #(the other spin densities because these electrons feel the effect of the _other_ spin species)
S_B = -U*ns["s_B"]
h_ud = np.array([
[S_A,0],
[0,S_B]
],dtype=complex)
#full hamiltonian
h = np.block([
[h_uu,h_ud],
[adj(h_ud),h_dd]
])
return h
In [5]:
def new_density(ns,mu_guess):
n_A_up_new = 0
n_B_up_new = 0
n_A_dn_new = 0
n_B_dn_new = 0
s_A_new = 0
s_B_new = 0
for k_idx in range(Nk):
k = ks[k_idx]
#calculate eigenstates
Es,us = np.linalg.eigh(h_eff(k,ns))
#calculate densities
for i in range(4): #number of bands is 2*2 spins
n_A_up_new += fermi(Es[i]-mu_guess) * adj(us[0,i])*us[0,i] #0th element is sublattice A, spin up
n_B_up_new += fermi(Es[i]-mu_guess) * adj(us[1,i])*us[1,i] #1st element is sublattice B, spin up
n_A_dn_new += fermi(Es[i]-mu_guess) * adj(us[2,i])*us[2,i] #2nd element is sublattice A, spin dn
n_B_dn_new += fermi(Es[i]-mu_guess) * adj(us[3,i])*us[3,i] #3rd element is sublattice B, spin dn
s_A_new += fermi(Es[i]-mu_guess) * adj(us[2,i])*us[0,i] #s_A flips a down spin on A to an up spin on A
s_B_new += fermi(Es[i]-mu_guess) * adj(us[3,i])*us[1,i] #s_B flips a down spin on B to an up spin on B
#normalize by number of k points
n_A_up_new = np.real(n_A_up_new/Nk)
n_B_up_new = np.real(n_B_up_new/Nk)
n_A_dn_new = np.real(n_A_dn_new/Nk)
n_B_dn_new = np.real(n_B_dn_new/Nk)
s_A_new = s_A_new/Nk #this can be complex
s_B_new = s_B_new/Nk #this can be complex
return [n_A_up_new,n_B_up_new,n_A_dn_new,n_B_dn_new,s_A_new,s_B_new]
def update_density(ns,mu_guess):
n_A_up_new,n_B_up_new,n_A_dn_new,n_B_dn_new,s_A_new,s_B_new = new_density(ns,mu_guess)
#update densities
ns["n_A_up"] = (1-delta)*ns["n_A_up"] + delta*n_A_up_new
ns["n_B_up"] = (1-delta)*ns["n_B_up"] + delta*n_B_up_new
ns["n_A_dn"] = (1-delta)*ns["n_A_dn"] + delta*n_A_dn_new
ns["n_B_dn"] = (1-delta)*ns["n_B_dn"] + delta*n_B_dn_new
ns["s_A"] = (1-delta)*ns["s_A"] + delta*s_A_new
ns["s_B"] = (1-delta)*ns["s_B"] + delta*s_B_new
In [6]:
def total_density(ns,mu_guess):
n_A_up,n_B_up,n_A_dn,n_B_dn,_,_ = new_density(ns,mu_guess)
return n_A_up + n_B_up + n_A_dn + n_B_dn
def diff(mu_guess):
return total_density(ns,mu_guess) - n_target
System Parameters¶
In [7]:
t1 = 1 #first hopping
t2 = 1 #second hopping
m = 0.5 #sublattice mass
U = 1 #Hubbard U
beta = 10**2 #cold limit
Iteration¶
In [ ]:
delta = 0.1 #step size
eps = 1e-5 #initial asymmetry to allow possible magnetic ordering
Nk = 21
ks = np.linspace(-np.pi,np.pi,Nk,endpoint=False)
iterations_0 = 100 #more iterations at each n to
iterations = 10
In [9]:
num_h = 101#21
num_n = 161#31
hz_vals = np.linspace(-4,4,num_h)
n_vals = np.linspace(0.01,3.99,num_n)
n_converged = np.zeros((num_h,num_n,4))
s_converged = np.zeros((num_h,num_n,2),dtype=complex)
mu_converged = np.zeros((num_h,num_n))
for n_idx in range(num_n):
n_target = n_vals[n_idx]
#initial guess at occupations
ns = {
"n_A_up" : n_target/4-eps,
"n_B_up" : n_target/4-eps,
"n_A_dn" : n_target/4+eps,
"n_B_dn" : n_target/4+eps,
"s_A" : eps,
"s_B" : eps
}
hz = 0
mu_guess = 0.0
#iterate
for s in range(iterations_0):
#find mu
if(diff(-10)*diff(10)<0):
pass
else:
raise(ValueError("mu not found"))
mu = find_root(diff,-10,10)
#compute
update_density(ns,mu)
#now vary h's which should be relatively "adiabatic" for fixed n
for h_idx in range(num_h):
hz = hz_vals[h_idx]
#iterate
for s in range(iterations):
#find mu
if(diff(-10)*diff(10)<0):
pass
else:
raise(ValueError("mu not found"))
mu = find_root(diff,-10,10)
#compute
update_density(ns,mu)
n_converged[h_idx,n_idx] = [ns["n_A_up"],ns["n_B_up"],ns["n_A_dn"],ns["n_B_dn"]]
s_converged[h_idx,n_idx] = [ns["s_A"],ns["s_B"]]
mu_converged[h_idx,n_idx] = mu
Plot results¶
In [13]:
fig = plt.figure(figsize=(4,2))
plt.imshow((mu_converged[:,1:-1]-mu_converged[:,0:-2])/((n_vals[-1]-n_vals[0])/num_n),vmin=0,vmax=5,aspect=0.62*(n_vals[-1]-n_vals[0])/(hz_vals[-1]-hz_vals[0]),origin="lower",interpolation=None,extent=[n_vals[0],n_vals[-1],hz_vals[0],hz_vals[-1]],cmap="magma")
plt.colorbar()
plt.xlabel(r"$n$")
plt.ylabel(r"$h_z$")
# plt.axvline(1,c="k")
# plt.axvline(2,c="k")
# plt.axvline(3,c="k")
plt.xticks([0,1,2,3,4])
plt.title(r"Inverse compressibility $d\mu/dn$")
plt.show()
In [15]:
fig = plt.figure(figsize=(4,2))
plt.imshow((n_converged[:,:,0]+n_converged[:,:,1]-n_converged[:,:,2]-n_converged[:,:,3]),vmin=-2,vmax=2,aspect=0.62*(n_vals[-1]-n_vals[0])/(hz_vals[-1]-hz_vals[0]),origin="lower",interpolation=None,extent=[n_vals[0],n_vals[-1],hz_vals[0],hz_vals[-1]],cmap="coolwarm")
plt.colorbar()
plt.xlabel(r"$n$")
plt.ylabel(r"$h_z$")
# plt.axvline(1,c="k")
# plt.axvline(2,c="k")
# plt.axvline(3,c="k")
plt.xticks([0,1,2,3,4])
plt.title(r"Magnetization $n_{A\!\!\uparrow}+n_{B\!\!\uparrow}-n_{A\!\!\downarrow}-n_{B\!\!\downarrow}$")
plt.show()
In [14]:
fig = plt.figure(figsize=(4,2))
plt.imshow((n_converged[:,:,0]-n_converged[:,:,1]+n_converged[:,:,2]-n_converged[:,:,3]),vmin=-0.5,vmax=0.5,aspect=0.62*(n_vals[-1]-n_vals[0])/(hz_vals[-1]-hz_vals[0]),origin="lower",interpolation=None,extent=[n_vals[0],n_vals[-1],hz_vals[0],hz_vals[-1]],cmap="terrain")
plt.colorbar()
plt.xlabel(r"$n$")
plt.ylabel(r"$h_z$")
# plt.axvline(1,c="k")
# plt.axvline(2,c="k")
# plt.axvline(3,c="k")
plt.xticks([0,1,2,3,4])
plt.title(r"Charge Imbalance $n_{A\!\!\uparrow}-n_{B\!\!\uparrow}+n_{A\!\!\downarrow}-n_{B\!\!\downarrow}$")
plt.show()
In [ ]:
In [ ]: