%run ../code/init_mooc_nb.py import scipy import scipy.linalg as sla from matplotlib import cm import warnings warnings.simplefilter("ignore", UserWarning) sigma0 = np.array([[1, 0], [0, 1]]) sigmax = np.array([[0, 1], [1, 0]]) sigmay = np.array([[0, -1j], [1j, 0]]) sigmaz = np.array([[1, 0], [0, -1]]) def evaluate_on_grid(X, Y, func): """ X, Y should be in np.meshgrid form. It's enough for func to work on floats. """ data = [] for xx, yy in zip(X, Y): row = [] for i,j in zip(xx, yy): row.append(func(i,j)) data.append(row) data = np.array(data) data = [np.array(data[:,:,i]) for i in range(np.shape(data)[2])] return data def get_energy_function(A, B, C): """ Used for plotting of hexagonal warping. """ def func(kx, ky): matrix = A*(kx**2+ky**2)*sigma0\ + B*(kx * sigmay - ky * sigmax) \ + C* 0.5 * ( (kx+1j*ky)**3 + (kx-1j*ky)**3 ) * sigmaz return sla.eigh(matrix)[0] return func # Onsite and hoppings matrices used for building BHZ model def onsite(site, par): A1, A2, B1, B2, D1, D2, C, M = par.A1, par.A2, par.B1, par.B2, par.D1, par.D2, par.C, par.M return (C + 2 * D1 + 4 * D2) * np.kron(sigma0, sigma0) + (M + 2 * B1 + 4 * B2) * np.kron(sigma0, sigmaz) def hopx(site1, site2, par): A1, A2, B1, B2, D1, D2, C, M = par.A1, par.A2, par.B1, par.B2, par.D1, par.D2, par.C, par.M return - D2 * np.kron(sigma0, sigma0) - B2 * np.kron(sigma0, sigmaz) + A2 * 0.5j * np.kron(sigmax, sigmax) def hopy(site1, site2, par): A1, A2, B1, B2, D1, D2, C, M = par.A1, par.A2, par.B1, par.B2, par.D1, par.D2, par.C, par.M return - D2 * np.kron(sigma0, sigma0) - B2 * np.kron(sigma0, sigmaz) + A2 * 0.5j * np.kron(sigmay, sigmax) def hopz(site1, site2, par): A1, A2, B1, B2, D1, D2, C, M = par.A1, par.A2, par.B1, par.B2, par.D1, par.D2, par.C, par.M return - D1 * np.kron(sigma0, sigma0) - B1 * np.kron(sigma0, sigmaz) + A1 * 0.5j * np.kron(sigmaz, sigmax) def bhz_slab(L, W, H): """A cuboid region of BHZ material with two leads attached. parameters for leads and scattering region can be defined separately """ lat = kwant.lattice.general(np.identity(3)) sys = kwant.Builder() sym = kwant.TranslationalSymmetry((1, 0, 0)) lead = kwant.Builder(sym) def shape_lead(pos): (x, y, z) = pos return (0 <= z < H) and (0 <= y < W) def shape(pos): (x, y, z) = pos return (0 <= z < H) and (0 <= y < W) and (0 <= x < L) def hopping_x(site1, site2, par): xt, yt, zt = site1.pos xs, ys, zs = site2.pos return hopx(site1,site2,par) * np.exp(-0.5j * par.Bz * (xt - xs) * (yt + ys)) # scattering system sys[lat.shape(shape, (0,0,0))] = lambda site, par: onsite(site, par.scat) - par.mu * np.eye(4) sys[kwant.HoppingKind((1,0,0), lat)] = lambda site1, site2, par: hopping_x(site1, site2, par.scat) sys[kwant.HoppingKind((0,1,0), lat)] = lambda site1, site2, par: hopy(site1, site2, par.scat) sys[kwant.HoppingKind((0,0,1), lat)] = lambda site1, site2, par: hopz(site1, site2, par.scat) # leads lead[lat.shape(shape_lead, (0,0,0))] = lambda site, par: onsite(site, par.lead) - par.mu_lead * np.eye(4) lead[kwant.HoppingKind((1,0,0), lat)] = lambda site1, site2, par: hopping_x(site1, site2, par.lead) lead[kwant.HoppingKind((0,1,0), lat)] = lambda site1, site2, par: hopy(site1, site2, par.lead) lead[kwant.HoppingKind((0,0,1), lat)] = lambda site1, site2, par: hopz(site1, site2, par.lead) # attaching leads sys.attach_lead(lead) sys.attach_lead(lead.reversed()) return sys.finalized() MoocVideo("62ZObitJ4DM", src_location="6.2-intro") # Making system Bz = 0. mu_lead = 0.7 mu_max = 0.4 L, W, H = 10, 30, 6 sys = bhz_slab(L, W, H) momenta = np.linspace(-np.pi/3, np.pi/3, 100) # calculating bands in scat par_lead = SimpleNamespace(A1=1.0, A2=1.0, B1=1.0, B2=1.0, C=0.0, D1=0., D2=0., M=-1.0, Bz=Bz) par = SimpleNamespace(lead=par_lead, scat=par_lead, mu_lead=0.0) bands = kwant.physics.Bands(sys.leads[0], args=[par]) energies_scat = [bands(k) for k in momenta] # calculating conductance mus = np.linspace(-mu_max, mu_max, 40) data = [] par_lead = SimpleNamespace(A1=1.0, A2=1.0, B1=1.0, B2=1.0, C=0.0, D1=0., D2=0., M=-1.0, Bz=0.0) par_scat = SimpleNamespace(A1=1.0, A2=1.0, B1=1.0, B2=1.0, C=0.0, D1=0., D2=0., M=-1.0, Bz=Bz) par = SimpleNamespace(lead=par_lead, scat=par_scat, mu_lead=mu_lead) sys_leads_fixed = sys.precalculate(energy=0.0, args=[par]) for par.mu in mus: smatrix = kwant.smatrix(sys_leads_fixed, energy=0.0, args=[par]) data.append(smatrix.transmission(1, 0)) # making plot gmax = 8 band_max = 1 def plot(mu): # Figures, axes and matplotlib defaults. fig, ax = plt.subplots(1,2, figsize=([9.5, 4]), tight_layout=True) ax1, ax2 = ax ax2.set_color_cycle(['k']) # Plotting conductance ax1.plot(mus, data, 'r-') ax1.plot([mu, mu], [0, gmax], 'b--') # Plotting bands ax2.plot(momenta, energies_scat) ax2.plot([-4, 4], [mu, mu], 'k--') # Setting labels and titles ax1.set_xlabel(r'$\mu$') ax1.set_ylabel(r'$G\,[e^2/h]$') ax2.set_xlabel('$k$') ax2.set_ylabel('$E$') ax1.set_title('conductance') ax2.set_title('spectrum') # Setting ticks vals = np.arange(-0.4, 0.8, 0.4) ax1.set_xticks(vals) ax1.set_xticklabels(["${0}$".format(i) for i in vals]) vals = range(0, gmax+1, 2) ax1.set_yticks(vals) ax1.set_yticklabels(["${0}$".format(i) for i in vals]) vals = np.arange(-2, 2.5, .5) ax2.set_yticks(vals) ax2.set_yticklabels(["${0}$".format(i) for i in vals]) ax2.set_xticks([-np.pi/3, 0.0, np.pi/3]) ax2.set_xticklabels([r'$-\frac{\pi}{3}$', r'$0$', r'$\frac{\pi}{3}$']) # Setting limits ax1.set_xlim(-mu_max, mu_max) ax1.set_ylim(0, gmax) ax2.set_xlim(-np.pi/3,np.pi/3); ax2.set_ylim(-band_max, band_max) return fig StaticInteract(lambda mu: plot(.1*mu), mu = RangeWidget(-3, 3, name='mu', default=0, show_range=True)) question = ("Which control parameter can remove the 0th plateau in the QHE measurement? ") answers = ["Increasing the magnetic field.", "Gate voltage difference (which controls difference in electron density) between the surfaces.", "Increasing topological mass.", "Adding an in-plane magnetic field."] explanation = ("Gate voltage difference changes the filling of the individual states without shifting the total density of electrons. " "This can therefore shift the plateaus of each surface. Magnetic field and topological mass are part of generating the " " $0^{th}$ plateau to begin with so cannot eliminate it. ") MoocMultipleChoiceAssessment(question=question, answers=answers, correct_answer=1, explanation=explanation) question = ("Why do you think ARPES observes surface states even if there is conductance through the bulk?") answers = ["ARPES can only observe occupied states and therefore bulk conductance is not an issue. ", "Since ARPES measures the spectrum in a momentum resolved way, it can separate out surface and bulk states.", "ARPES does not measure conductance and therefore bulk electronic states are not an issue.", "Since ARPES measures the spectrum in an energy resolved way, it can selectively choose the surface states in the bulk gap."] explanation = ("The surface states live within the energy gap of the bulk. Since ARPES directly measure $E(k)$, it separates out " "surface states from bulk states, which are in different energy ranges. ") MoocMultipleChoiceAssessment(question=question, answers=answers, correct_answer=3, explanation=explanation) PreprintReference("0908.1418", show_abstract=False) def plot_warping(A=0.0, B=1.0, C=1.0, Kmax=0.2): # Plot parameters N = 100 xlim, ylim = 1.2, 1.2 zmin, zmax = -1.0, 3.5 r = np.linspace(0,Kmax,N) p = np.linspace(0,2*np.pi,N) R,P = np.meshgrid(r,p) X,Y = R*np.cos(P),R*np.sin(P) energies = evaluate_on_grid(X, Y, func=get_energy_function(A, B, C)) fig = plt.figure(figsize=(12, 8)) ax = fig.add_subplot(111, projection='3d') vmin = np.array(energies).min() vmax = np.array(energies).max() for z in energies: ax.plot_surface(X, Y, z, rstride=1, cstride=1, cmap=cm.RdBu_r, linewidth=0.05, vmin=vmin, vmax=vmax) ax.set_xlabel('$k_x$')# [\sqrt{v/\lambda}]$') ax.set_ylabel('$k_y$')# [\sqrt{v/\lambda}]$') ax.set_zlabel('$E$')# [E^{*}]$') ax.set_xlim3d(-xlim, xlim) ax.set_ylim3d(-ylim, ylim) ax.set_zlim3d(zmin, zmax) vals = [-xlim,0,xlim] ax.set_xticks(vals) ax.set_xticklabels(["${0}$".format(i) for i in vals]) vals = [-ylim,0,ylim] ax.set_yticks(vals) ax.set_yticklabels(["${0}$".format(i) for i in vals]) vals = [zmin, 0, zmax] ax.set_zticks(vals) ax.set_zticklabels(["${0}$".format(i) for i in vals]) ax.view_init(20,15) return ax plot_warping(A=1.2, B=1.8, C=1.5, Kmax=1.0); MoocVideo("WZmNeEwM1N4", src_location="6.2-summary") MoocDiscussion("Questions", "3DTI materials and signatures")