In [1]:
import numpy as np
import ipywidgets as ipw
from bokeh.io import push_notebook, output_notebook, show
from bokeh.layouts import row, column
from bokeh.plotting import figure
from bokeh.models import ColumnDataSource
output_notebook()
from ipywidgets import interact
from collections import OrderedDict
old_settings = np.seterr(over = 'ignore') #Ignore warnings about overflow data points
Start by setting some parameters including an (artificial) analytical expression of the Fermi surface along kz direction¶
In [2]:
kz = np.linspace(-np.pi, np.pi);
Ez = -2*np.cos(kz)-3.0;
B = 0.2; # Magnetic field
hbar = 1.0;
e = 1.0; # Electron charge
m = 1.0; # (free) Electron mass
c = 1.0; # Speed of light
wc = e*B/(m*c); # Cyclotron frequency
In [3]:
# Define the Fermi surface in zero magnetic field
freeE = figure(height=400, width=400, title="Dispersion relationship of a free electron",
tools="pan,reset,save,wheel_zoom", x_range = [1.1*np.amin(Ez), -1.1*np.amin(Ez)],
y_range = [np.amin(kz), np.amax(kz)]);
source1 = ColumnDataSource(data = {'xVal':-Ez, 'yVal':kz, 'xxVal':Ez});
freeE.line('xVal', 'yVal', source = source1, line_width=2, line_alpha=1);
freeE.line('xxVal', 'yVal', source = source1, line_width=2, line_alpha=1);
freeE.xaxis.axis_label = 'kx & ky';
freeE.yaxis.axis_label = 'kz';
show(freeE)
Define some expressions according to Landau quantization, such as the energy spacing and projection area in k-space¶
In [4]:
dE = hbar*wc;
dA = 2*m*np.pi*wc/hbar;
A = np.arange(np.amin(Ez)**2*np.pi, np.amax(Ez)**2*np.pi,-dA);
kk = np.sqrt(A/np.pi);
kF = np.array([np.flip(-kk),kk]);
kF = kF.flatten();
Ef = np.arccos((kk-3.0)/2.0);
Ef = np.array([np.flip(-Ef),Ef]);
Ef = Ef.flatten();
gExy = m/(np.pi*hbar**2);
source2 = ColumnDataSource(data = {'xVal':kF, 'yVal': Ef, 'yyVal': -Ef});
freeE.circle('xVal','yVal', source = source2, line_width = 3, line_alpha = 0.2);
freeE.circle('xVal','yyVal', source = source2, line_width = 3, line_alpha = 0.2);
freeE.xaxis.axis_label = 'kx & ky';
freeE.yaxis.axis_label = 'kz';
Set up a histogram to count interceptions between the Fermi surface and "Landau tubes", compute density of states¶
In [5]:
# Plot the histogram of states per energy interval
# # Option 1: use adaptive bin width
# dEf = np.diff(Ef);
# bWidth = np.amax(np.abs(dEf)); # Define the bin width of the histogram so minimum it contains one kxy tube
# # Option 2: use fixed bin size to slice
bWidth = (np.amax(kz)-np.amin(kz))/100.0;
histoN = (2*np.amax(kz)/bWidth).astype(int); # calculate the number of bars in the histogram
kCount = np.empty(histoN, dtype = np.int8); # Declare an array to store histogram information
gEs = np.empty(histoN, dtype = np.double); # Declare an array to store density of stats
kIntv = np.linspace(kz.min(),kz.max(), histoN);
zero = np.zeros(histoN);
for idx in range(histoN):
kCount[idx] = len(Ef[ (Ef >= kIntv[idx-1]) & (Ef < kIntv[idx]) ])
gEs[idx] = kCount[idx]*gExy/(0.001+np.abs(np.sin(0.5*(kIntv[idx-1]+kIntv[idx]))));
histo = figure(height=400, width=300, title="allowed states per Energy interval",
tools="pan,reset,save,wheel_zoom", y_range = [np.amin(kz), np.amax(kz)],
x_range = [0, 1.1*np.amax(gEs)]);
histo.xaxis[0].axis_label='Number of states';
histo.yaxis[0].axis_label='Energy (eV)';
source3 = ColumnDataSource(data = {'bottom': kIntv, 'top': kIntv + bWidth, 'right':kCount, 'left': zero});
source4 = ColumnDataSource(data = {'bottom': kIntv, 'top': kIntv + bWidth, 'right':gEs, 'left': zero})
histo.quad(bottom = 'bottom', top = 'top', right = 'right', left = 'left', source = source3,
fill_color="navy", line_color="white", alpha=1);
histo2 = figure(height=400, width=300, title = "extrema orbitals",
tools="pan,reset,save,wheel_zoom", y_range = [np.amin(kz), np.amax(kz)],
x_range = [0, 1.1*np.amax(gEs)]);
histo2.quad(bottom = 'bottom', top = 'top', right = 'right', left = 'left', source = source4,
fill_color="navy", line_color="white", alpha=1);
show(row(freeE,histo,histo2))
Haas-van Alphen effect¶
Next, we use the above method to compute the density of states at different external fields to visualize the de Haas-van Alphen effect¶
In [ ]:
# Remove variables to avoid conflicts
%reset
In [ ]:
import numpy as np
import ipywidgets as ipw
from bokeh.io import push_notebook, output_notebook, show
from bokeh.layouts import row, column
from bokeh.plotting import figure
from bokeh.models import ColumnDataSource
output_notebook()
from ipywidgets import interact
from collections import OrderedDict
old_settings = np.seterr(over = 'ignore') #Ignore warnings about overflow data points
Wrap everything inside a function that takes a value of external field and computes the density of states¶
In [ ]:
# Wrap everything in a function to compute densities of states at different external magnetic fields
def LandauDen(B):
kz = np.linspace(-np.pi, np.pi);
Ez = -2*np.cos(kz)-3.0;
hbar = 1.0;
e = 1.0; # Electron charge
m = 1.0; # (free) Electron mass
c = 1.0; # Speed of lightDensity of States in Magnetic Fields Using Landau Levels
wc = e*B/(m*c); # Cyclotron frequency
dE = hbar*wc;
A = np.arange(np.amin(Ez)**2*np.pi, np.amax(Ez)**2*np.pi,-2*m*np.pi*wc/hbar);
kk = np.sqrt(A/np.pi);
kF = np.array([np.flip(-kk),kk]);
kF = kF.flatten();
Ef = np.arccos((kk-3.0)/2.0);
Ef = np.array([np.flip(-Ef),Ef]);
Ef = Ef.flatten();
gExy = m/(np.pi*hbar**2);
bWidth = (np.amax(kz)-np.amin(kz))/100.0; # Fix the bin width\
histoN = (2*np.amax(kz)/bWidth).astype(int); # calculate the number of bars in the histogram
gEs = np.empty(histoN, dtype = np.double); # Declare an array to store density of stats
haas = np.empty(histoN, dtype = np.double);
kIntv = np.linspace(kz.min(),kz.max(), histoN);
zero = np.zeros(histoN);
for idx in range(histoN-1):
temp = len(Ef[ (Ef >= kIntv[idx]) & (Ef <= kIntv[idx+1]) ]);
gEs[idx] = temp*gExy/(0.001+np.abs(np.sin(0.5*(kIntv[idx+1]+kIntv[idx]))));
haas = abs(np.sin(Ef**-2));
return np.sum(gEs)
Generate an array of sampling points in magnetic field¶
In [ ]:
invB = np.arange(0.1, 1, 0.001);
B = 1/invB;
densStat = np.empty(np.size(B),dtype = np.double);
for ii in range(np.size(densStat)):
densStat[ii] = LandauDen(B[ii]);
In [ ]:
source5 = ColumnDataSource(data = {'xVal':B, 'yVal': densStat});
source6 = ColumnDataSource(data = {'xVal':invB, 'yVal': densStat});
dHvA = figure(height=350, width=450, title="Density of states of free electrons in magnetic field");
dHvA.line('xVal','yVal', source = source5, line_width = 2, line_alpha = 1);
dHvA.xaxis.axis_label = 'Magnetic Field (B)';
dHvA.yaxis.axis_label = 'Density of States';
dHvA2 = figure(height=350, width=450, title="Density of states of free electrons in magnetic field");
dHvA2.line('xVal','yVal', source = source6, line_width = 2, line_alpha = 1);
dHvA2.xaxis.axis_label = 'Inverse of Magnetic Field (1/B)';
dHvA2.yaxis.axis_label = 'Density of States';
show(row(dHvA,dHvA2))
In [ ]:
In [ ]: