Constant energy contours and QPI model¶

Calculation of joint density of states (JDoS) from the DFT-calculated constant energy contours (CECs).

Import data¶

In [1]:

import numpy as np
import matplotlib.pyplot as plt
import os
from tqdm import tqdm

In [2]:

# load bandstructure data
dKGM = np.loadtxt('data/CEC_data/KGM_band_ssc.dat')
dKGM.shape

Out[2]:

(64000, 4)

In [3]:

print('|M| = ', dKGM[:,0].max(), '1/a_B; |K|=', dKGM[:,1].max(), '1/a_B')

|M| =  0.46284401 1/a_B; |K|= 0.53061986 1/a_B

In [4]:

from masci_tools.util.constants import BOHR_A

a = 21.551001
b = 4.177400
c = 51.152020

2*np.pi/(a/BOHR_A) * 3, 2*np.pi/(b/BOHR_A) * 2/3

Out[4]:

(0.4628441819189051, 0.530620078803855)

In [5]:

# import files for constant energy contours
fnames = [i for i in os.listdir('data/CEC_data/') if '_band_ssc.dat' in i and 'KGM' not in i]
d = []
for fname in tqdm(fnames):
    tmp = np.loadtxt('data/CEC_data/'+fname)
    d += list(tmp)
    
# format: kx, ky, energy, weight
d = np.array(d)

100%|██████████| 47/47 [00:03<00:00, 14.53it/s]

In [7]:

def E_k_tf(k, E0, D, V, vF, Delta, B, lmb=0, theta=0):
    """Model from eq. (2) of Zhang et al, Nat. Physics 6, 584-588 (2010)"""
    hbar = 1 #6.5821220e-16 # eV s
    E1 = E0 - D*k**2
    de_plus = np.sqrt((abs(V) + vF*hbar*abs(k))**2 + (Delta/2 - B*k**2)**2 + (lmb*k**3 * np.cos(3*theta))**2)
    de_minus = np.sqrt((abs(V) - vF*hbar*abs(k))**2 + (Delta/2 - B*k**2)**2 + (lmb*k**3 * np.cos(3*theta))**2)
    E4 = E1 - de_minus
    E3 = E1 - de_plus
    E2 = E1 + de_minus
    E1 = E1 + de_plus
    return np.array([E1, E2, E3, E4], dtype=np.float128)

In [ ]:

$\hbar^2/2m_e = 38.09982350\,\mathrm{meV}\,\mathrm{nm^2}$

In [8]:

hbar_by_2me_eVAng = 38.09982350*1e-3*1e-2 # eV/\AA^2

In [9]:

def E(k, theta,
      lmb = 30, #eV A3
      ED = -0.89, #eV
      vD = 3.85, #eV A
      Delta = 0.31, #eV
      hbar_by_2me= 0,
     ):
    vD /= BOHR_A
    lmb /= BOHR_A**3
    
    hbar_by_2me /= BOHR_A**2
        
    return ED + k**2*hbar_by_2me + ((vD*k)**2 + (lmb*k**3 * np.cos(3*theta))**2 + Delta**2)**(1/2.)

In [14]:

import numpy as np
from tqdm import tqdm

def get_d_with_symm(E0, dE):

    # filter energy cut
    dist = abs(d[:,2]-E0)
    dd = d[dist<dE]

    # use symmetries
    d_with_symm = list(dd.copy())
    tmp = dd.copy(); tmp[:,0] *= -1
    d_with_symm += list(tmp)
    tmp = dd.copy(); tmp[:,1] *= -1
    d_with_symm += list(tmp)
    tmp = dd.copy(); tmp[:,0] *= -1; tmp[:,1] *= -1
    d_with_symm += list(tmp)
    d_with_symm = np.array(d_with_symm)

    return d_with_symm

In [15]:

def plot_CEC(E0 = -0.5, dE = 2e-2, s0=1, cmap='binary', levels=None):
    
    if E0 is not None and dE is not None:
        fname_jdos = f'save_jdos_E0={E0:.3f}_dE={dE:.3f}.npz'
        d_with_symm = get_d_with_symm(E0, dE)
        
        d_with_symm = d_with_symm[abs(d_with_symm[:,0])<=0.25]
        d_with_symm = d_with_symm[abs(d_with_symm[:,1])<=0.25]
        print(d_with_symm.shape)
        d_with_symm = np.array(list(d_with_symm)+[[-0.5, 0, 0, 1e-9], [0.5, 0, 0, 1e-9]])
        
        
        c = d_with_symm[:,-1]; c = (c-c.min())/(c.max()-c.min())
        
        if levels is not None:
            plt.contour(ks, ks, E_all[:,:,0], levels=levels, linewidths=2, alpha=0.5, colors='red', zorder=100, linestyles='--')
            plt.contour(ks, ks, E_all[:,:,1], levels=levels, linewidths=2, alpha=0.5, colors='blue', zorder=100, linestyles='--')
    
        plt.scatter(d_with_symm[:,0], d_with_symm[:,1], c=np.sqrt(c), s=s0*c, alpha=np.sqrt(c), cmap=cmap, rasterized=True) #lw=0.5, ec='k',)
        plt.axis('scaled')
        plt.xlabel('$k_x$')
        plt.ylabel('$k_y$')
        
        fac = 1.570129103783706
        plt.ylim(-0.5/fac, 0.5/fac)
        plt.xlim(-0.5, 0.5)
        
        plt.colorbar()

Fig. 10 (Appendix)¶

In [ ]:

In [10]:

Nx = Ny = 1000

ks = np.linspace(-0.5, 0.5, Nx)


vD, lmb, D, B, Delta, E0 = 6.78, 41.07, 28, -29.9, -2*0.51, -1.10
vF = vD / BOHR_A
lmb /= BOHR_A**3


# new values from Jens
E0 = -1.2
vD =7.2
V = 0.05
lmb = 26
Delta = -1.24
D = 7.8
B = -8.4
D /= BOHR_A**2
B /= BOHR_A**2
vF = vD / BOHR_A
lmb /= BOHR_A**3


# our data
E_all = []
for kx in tqdm(ks):
    for ky in ks:
        k = np.sqrt(kx**2+ky**2)
        theta = np.arctan2(ky, kx)
        # E_all.append(E(k, theta,
        #                   lmb = 30, #eV A3
        #                   ED = -0.89, #eV
        #                   vD = 3.85, #eV A
        #                   Delta = 0.31, #eV
        #               ))
        
        E1, E2, E3, E4 = E_k_tf(k, E0, D, V, vF, Delta, B, lmb, theta)
        
        E_all.append([E1, E2])

E_all = np.array(E_all).reshape(Nx, Ny, -1)

100%|██████████| 1000/1000 [00:19<00:00, 50.15it/s]

In [11]:

# data from Kuroda et al., PRL 105, 076802 (2010))
# Bi2Se3

E_all_2 = []
for kx in tqdm(ks):
    for ky in ks:
        k = np.sqrt(kx**2+ky**2)
        theta = np.arctan2(ky, kx)
        E_all_2.append(E(k, theta,
                          lmb = 128, #eV A3
                          ED = 0, #eV
                          vD = 3.55, #eV A
                          Delta = 0., #eV
                        ))
E_all_2 = np.array(E_all_2).reshape(Nx, Ny)

100%|██████████| 1000/1000 [00:07<00:00, 128.47it/s]

In [12]:

# data from Fu, 
# Bi2Te3

E_all_3 = []
for kx in tqdm(ks):
    for ky in ks:
        k = np.sqrt(kx**2+ky**2)
        theta = np.arctan2(ky, kx)
        E_all_3.append(E(k, theta,
                          lmb = 250, #eV A3
                          ED = 0., #eV
                          vD = 2.55, #eV A
                          Delta = 0., #eV
                          hbar_by_2me=0, #hbar_by_2me_eVAng
                        ))
E_all_3 = np.array(E_all_3).reshape(Nx, Ny)

100%|██████████| 1000/1000 [00:07<00:00, 128.20it/s]

In [13]:

# m = np.where(abs(E_all-0.2)<1e-3)
# plt.pcolormesh(E_all)

plt.figure(figsize=(5,8))
# plt.suptitle(r'$E_k=\sqrt{v_F^2k^2+\lambda^2k^6\cos^2(3\theta)+\Delta^2}$', fontsize='large')

plt.subplot(3,1,1)
plt.contour(ks*BOHR_A, ks*BOHR_A, E_all[:,:,0]-E_all[:,:,0].min(), levels=[0.01, 0.2, 0.4, 0.7, 1.2], linewidths=2)
plt.contour(ks*BOHR_A, ks*BOHR_A, E_all[:,:,1]-E_all[:,:,1].min(), levels=[0.01, 0.2, 0.4, 0.7, 1.2], linewidths=2)
# plt.title('this work:\n' + r'$v=3.85\,\mathrm{eV\AA}, \lambda=30\,\mathrm{eV\AA^3}, \Delta=0.31\,\mathrm{eV}$')
plt.colorbar()
plt.axis('equal')
# plt.xlabel(r'$k_x$ (1/$\mathrm{\AA}$)')
plt.ylabel(r'$k_y$ (1/$\mathrm{\AA}$)')
plt.xlim(-0.1, 0.1)
plt.ylim(-0.1, 0.1)

plt.subplot(3,1,2)
plt.contour(ks*BOHR_A, ks*BOHR_A, E_all_2-E_all_2.min(), levels=[0.01, 0.2, 0.4, 0.7, 1.2], linewidths=2)
# plt.title('Bi$_2$Se$_3$ [Kuroda2010]:\n' + r'$v=3.55\,\mathrm{eV\AA}, \lambda=128\,\mathrm{eV\AA^3}, \Delta=0\,\mathrm{eV}$')
plt.colorbar()
plt.axis('equal')
# plt.xlabel(r'$k_x$ (1/$\mathrm{\AA}$)')
plt.ylabel(r'$k_y$ (1/$\mathrm{\AA}$)')
plt.xlim(-0.1, 0.1)
plt.ylim(-0.1, 0.1)

arrkwargs = dict(arrowprops={'arrowstyle': '<-'}, ha='center', va='center')
plt.annotate(xy=(-0.15, -0.08-0.002), xytext=(-0.15, -0.08+0.06), text=r'$\overline{K}$', **arrkwargs)
plt.annotate(xy=(-0.15-0.002, -0.08), xytext=(-0.15+0.05, -0.08), text=r'$\overline{M}$', **arrkwargs)

plt.subplot(3,1,3)
plt.contour(ks*BOHR_A, ks*BOHR_A, E_all_3-E_all_3.min(), levels=[0.01, 0.2, 0.4, 0.7, 1.2], linewidths=2)
# plt.title('Bi$_2$Te$_3$ [Fu2009]:\n' + r'$v=2.55\,\mathrm{eV\AA}, \lambda=250\,\mathrm{eV\AA^3}, \Delta=0\,\mathrm{eV}$')
plt.colorbar()
plt.axis('equal')
plt.xlabel(r'$k_x$ (1/$\mathrm{\AA}$)')
plt.ylabel(r'$k_y$ (1/$\mathrm{\AA}$)')
plt.xlim(-0.1, 0.1)
plt.ylim(-0.1, 0.1)

plt.tight_layout()

plt.annotate(xy=(0.03, 0.960), text=r'(a)', xycoords='figure fraction', fontsize='large')
plt.annotate(xy=(0.03, 0.650), text=r'(b)', xycoords='figure fraction', fontsize='large')
plt.annotate(xy=(0.03, 0.330), text=r'(c)', xycoords='figure fraction', fontsize='large')


plt.savefig('dispersion_comparison.pdf')

plt.show()

No description has been provided for this image

Fig. 3 - bandstructure¶

In [24]:

from matplotlib.gridspec import GridSpec

# look over main figure and inset
for inset in [False, True]:

    fig = plt.figure(figsize=(10,8))

    gs = GridSpec(3,4, figure=fig)

    ax1 = fig.add_subplot(gs[0:2, 1:3])

    c = ddd[:,3]; c = (c-c.min())/(c.max()-c.min())
    plt.scatter(x, y, c=np.sqrt(c), s=c*50, alpha=np.sqrt(c), lw=0.5, ec='k', cmap='binary', rasterized=True)
    cbar = plt.colorbar()
    cbar.set_label('surface weight')
    plt.ylabel('$E-E^\mathrm{exp}_\mathrm{F}$ (eV)', fontsize='large')
    plt.ylim(-2,1)
    plt.xlim(x.min(), x.max())
    plt.xticks([x.min(), 0., x.max()], ['K', '$\Gamma$', 'M'], fontsize='large')

    # for zoom
    if inset:
        plt.ylim(-0.6, -0.35)
        plt.xlim(-0.07, 0.07)

    k = np.linspace(0, 1, 1000)
    
    # fitted values from Jens
    E0 = -1.2
    vD =7.2
    V = 0.05
    lmb = 26
    Delta = -1.24
    D = 7.8
    B = -8.4

    plt.title(fr'$E_D={np.round(E0,2)}$, $v_D={np.round(vD,2)}$, $D={np.round(D,2)}$, $\Delta={np.round(Delta/2,2)}$, $B={np.round(B,2)}$, $V_z={np.round(V,2)}$, $\lambda={np.round(lmb,2)}$')

    vF = vD / BOHR_A
    lmb /= BOHR_A**3
    D /= BOHR_A**2
    B /= BOHR_A**2

    E1, E2, E3, E4 = E_k_tf(k, E0, D, V, vF, Delta, B)
    E1_2, E2_2, E3_2, E4_2 = E_k_tf(-k, E0, D, V, vF, Delta, B, theta=np.pi, lmb=lmb)

    # electron like branches
    m = np.where(E1<=100) #<=0.35)
    plt.plot(k[m], E1[m], color='r', ls='--')
    m = np.where(E2<=100) #<=0.35)
    plt.plot(k[m], E2[m], color='b', ls='--')
    m = np.where(E1_2<=0.35)
    plt.plot(-k[m], E1_2[m], color='r', ls='--')
    m = np.where(E2_2<=0.35)
    plt.plot(-k[m], E2_2[m], color='b', ls='--')

    # hole like branches
    m = np.where(E3<=0.35)
    plt.plot(k[m], E3[m], color='r', ls='--')
    m = np.where(E4<=0.35)
    plt.plot(k[m], E4[m], color='b', ls='--')
    m = np.where(E3_2<=0.35)
    plt.plot(-k[m], E3_2[m], color='r', ls='--')
    m = np.where(E4_2<=0.35)
    plt.plot(-k[m], E4_2[m], color='b', ls='--')


    plt.axhline(0, color='r', ls='--')

    ax2 = fig.add_subplot(gs[2, 0:2])
    p = plot_CEC(E0=-0.2 - 0.2, dE=0.02, s0=5, levels=[-0.2])



    ax3 = fig.add_subplot(gs[2, 2:])
    
    if not inset:
        p = plot_CEC(E0=-0.2 - 0.05, dE=0.02, s0=5, levels=[-0.05])
    else:
        p = plot_CEC(E0=-0.2 + 0.05, dE=0.02, s0=5, levels=[0.05])


    arrkwargs = dict(arrowprops={'arrowstyle': '<-'}, ha='center', va='center')
    plt.annotate(xy=(-0.2, -0.2), xytext=(-0.2+0.15, -0.2), text=r'$\overline{K}$', **arrkwargs)
    plt.annotate(xy=(-0.2, -0.2), xytext=(-0.2, -0.2+0.15), text=r'$\overline{M}$', **arrkwargs)

    # if not inset:
    #     plt.savefig('panel_Fig3_new.svg', dpi=300)
    # else:
    #     plt.savefig('panel_Fig3_inset.svg', dpi=300)

    plt.show()

<>:16: SyntaxWarning: invalid escape sequence '\m'
<>:19: SyntaxWarning: invalid escape sequence '\G'
<>:16: SyntaxWarning: invalid escape sequence '\m'
<>:19: SyntaxWarning: invalid escape sequence '\G'
/tmp/ipykernel_999016/962781318.py:16: SyntaxWarning: invalid escape sequence '\m'
  plt.ylabel('$E-E^\mathrm{exp}_\mathrm{F}$ (eV)', fontsize='large')
/tmp/ipykernel_999016/962781318.py:19: SyntaxWarning: invalid escape sequence '\G'
  plt.xticks([x.min(), 0., x.max()], ['K', '$\Gamma$', 'M'], fontsize='large')

(612, 4)
(964, 4)

(612, 4)
(1440, 4)

In [ ]: