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

adj = lambda u: np.conj(u.T)

def fermi(z):
    if(beta*z<-30):
        return 1
    elif(beta*z>30):
        return 0
    else:
        return 1/(1+np.exp(beta*z))

System Parameters

t1 = 1 #first hopping
t2 = 1 #second hopping
m = 0.5 #sublattice mass
U = 1 #Hubbard U
hz = 0.2 #zeeman splitting term
beta = 10**2 #cold limit
n_target = 1 #3/4 filling

delta = 0.1 #step size
eps = 1e-5 #initial asymmetry to allow possible magnetic ordering
Nk = 101
ks = np.linspace(-np.pi,np.pi,Nk,endpoint=False)
iterations = 100

Effective Hamiltonian and Step

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

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

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

Iteration

#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
}
mu_guess = 0.0

#iterate
densities = np.zeros((iterations,4))
spin_flips = np.zeros((iterations,2),dtype=complex)
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)
    densities[s] = [ns["n_A_up"],ns["n_B_up"],ns["n_A_dn"],ns["n_B_dn"]]
    spin_flips[s] = [ns["s_A"],ns["s_B"]]

Evolution of Densities

plt.plot(densities[:,0],label=r"$n_{A,\!\!\!\uparrow}$")
plt.plot(densities[:,1],label=r"$n_{B,\!\!\!\uparrow}$")
plt.plot(densities[:,2],ls="dashed",label=r"$n_{A,\!\!\!\downarrow}$")
plt.plot(densities[:,3],ls="dashed",label=r"$n_{B,\!\!\!\downarrow}$")
plt.xlabel("iteration")
plt.ylabel("occupation")
plt.legend(fontsize=9)
plt.ylim(0.1,0.405)

(0.1, 0.405)

plt.plot(np.real(spin_flips[:,0]),label=r"Re $s_A$")
plt.plot(np.real(spin_flips[:,1]),ls="dashed",label=r"Re $s_B$")
plt.xlabel("iteration")
plt.ylabel("spin flip coherence")
plt.legend()

<matplotlib.legend.Legend at 0x169291c70>

Free and HF Energy Bands

eigs_free = np.zeros((Nk,4))
eigs_hf = np.zeros((Nk,4))
for k_idx in range(Nk):
    eigs_free[k_idx] = np.linalg.eigvalsh(h_eff(ks[k_idx],{"n_A_up":0,"n_B_up":0,"n_A_dn":0,"n_B_dn":0,"s_A":0,"s_B":0}))
    eigs_hf[k_idx] = np.linalg.eigvalsh(h_eff(ks[k_idx],ns))

plt.plot(ks,eigs_free[:,0],c="r",label=r"$\epsilon_{0}$")
for n in range(1,4): #2 bands * 2 spins
    plt.plot(ks,eigs_free[:,n],c="r")

plt.plot(ks,eigs_hf[:,0]-1/4*U,ls="dashed",c="b",label=r"$\epsilon_{HF}$")
for n in range(1,4): #2 bands * 2 spins
    plt.plot(ks,eigs_hf[:,n]-1/4*U,ls="dashed",c="b")

plt.xticks([-np.pi,0,np.pi],[r"$-\pi$",r"$0$",r"$\pi$"])
plt.xlabel(r"$k$")
plt.ylabel(r"Energy")
plt.legend(loc=1)

<matplotlib.legend.Legend at 0x16935d400>