% Fig_1_17.m % Essential Semiconductor Laser Device Physics % This program calculates the bulk band structure for Si, Ge, AlAs, % GaAs, InP, or InAs, along the L - Gamma - X - U,K - Gamma symmetry % directions. The program prompts the user to enter the material to be % modeled. The method used is the sp3s* model as described by % Vogl et. al. in J. Phys. Chem. Solids 44, 365 (1983). The parameters % used are as given in the paper. clear all clf; %plotting parameters + fontsizes FS = 16; %label fontsize FSN = 14; %number fontsize LW = 2; %linewidth % Change default axes fonts. set(0,'DefaultAxesFontName', 'Times New Roman'); set(0,'DefaultAxesFontSize', FSN); % Change default text fonts. set(0,'DefaultTextFontname', 'Times New Roman'); set(0,'DefaultTextFontSize', FSN); hbarChar=['\fontname{MT Extra}h\fontname{Times New Roman}']; %******************************* % Declare Constants and Arrays * %******************************* % Constants eye = complex(0,1); %Imaginary number, sqrt(-1) d1 = [ 1 1 1]/4; %Nearest neighbor directions d2 = [ 1 -1 -1]/4; d3 = [-1 1 -1]/4; d4 = [-1 -1 1]/4; N = 301; % Number of points calculated for each symmetry direction % Parameters for the sp3s* model % Si Ge AlAs GaAs InP InAs soa = [ 0 0 0.4210 0.4210 0.0670 0.4210]/3; soc = [ 0 0 0.0240 0.1740 0.3920 0.3920]/3; Esa = [-4.2000 -5.8800 -7.5273 -8.3431 -8.5274 -9.5381]; Epa = [ 1.7150 1.6100 0.9833 1.0414 0.8735 0.9099]; Esc = [-4.2000 -5.8800 -1.1627 -2.6569 -1.4826 -2.7219]; Epc = [ 1.7150 1.6100 3.5867 3.6686 4.0465 3.7201]; Esea = [ 6.6850 6.3900 7.4833 8.5914 8.2635 7.4099]; Esec = [ 6.6850 6.3900 6.7267 6.7386 7.0665 6.7401]; Vss = [-8.3000 -6.7800 -6.6642 -6.4513 -5.3614 -5.6052]; Vxx = [ 1.7150 1.6100 1.8780 1.9546 1.8801 1.8398]; Vxy = [ 4.5750 4.9000 4.2919 5.0779 4.2324 4.4693]; Vsapc = [ 5.7292 5.4649 5.1106 4.4800 2.2265 3.0354]; Vpasc = [ 5.7292 5.4649 5.4965 5.7839 5.5825 5.4389]; Vseapc = [ 5.3749 5.2191 4.5216 4.8422 3.4623 3.3744]; Vpasec = [ 5.3749 5.2191 4.9950 4.8077 4.4814 3.9097]; %***************** % Get User Input * %***************** mat = menu('Choose a material.','Si','Ge','AlAs','GaAs','InP','InAs') clc switch mat case 1 matname = 'Si'; case 2 matname = 'Ge'; case 3 matname = 'AlAs'; case 4 matname = 'GaAs'; case 5 matname = 'InP'; case 6 matname = 'InAs'; end uinput = 0; while uinput == 0 direction = input('Type 1 for L-\Gamma-X-U,K-\Gamma symmetry directions, or 2 for arbitrary direction.') if direction == 1 uinput = 1; elseif direction == 2 uinput = 1; else clc disp('Please choose from the options listed.') end end clc if direction == 1 % l, m, and n vectors for Gamma-L, Gamma-X, X-U, % and Gamma-K directions, respectively l = [0.5 1 0 1]; m = [0.5 0 1 1]; n = [0.5 0 1 0]; kmax = [2*pi 2*pi 2*pi 1.5*pi]; kmin = [0 0 1.5*pi 0 ]; else % Input l uinput = 0; while uinput == 0 l = input('Please specify l value'); if isnumeric(l) uinput = 1; else clc disp('Please enter a number') end end clc % Input m uinput = 0; while uinput == 0 m = input('Please specify m value') if isnumeric(m) uinput = 1; else clc disp('Please enter a number') end end clc % Input n uinput = 0; while uinput == 0 n = input('Please specify n value') if isnumeric(n) uinput = 1; else clc disp('Please enter a number') end end clc % Input kmax uinput = 0; while uinput == 0 kmax = input('Specify maximum value of k') if isnumeric(kmax) uinput = 1; else clc disp('Please enter a number') end end clc % Input kmin uinput = 0; while uinput == 0 kmin = input('Specify minimum value of k') if isnumeric(kmin) uinput = 1; else clc disp('Please enter a number') end end clc end %********************** % Execute Main Script * %********************** for d = 1:length(l) % Symmetry direction loop for i = 1:N % Band calculation loop % Calculate k vector based on symmetry direction k = [l(d) m(d) n(d)]*(kmin(d) + (kmax(d)-kmin(d))*(N-i)/(N-1)); p1 = exp(eye*sum(k.*d1)); p2 = exp(eye*sum(k.*d2)); p3 = exp(eye*sum(k.*d3)); p4 = exp(eye*sum(k.*d4)); g0 = (p1+p2+p3+p4)/4; g1 = (p1+p2-p3-p4)/4; g2 = (p1-p2+p3-p4)/4; g3 = (p1-p2-p3+p4)/4; % Create Hamiltonian % |sa> |sc> |pxa> |pya> |pza> |pxc> |pyc> |pzc> |sea> |sec> h = [Esa(mat)/2 Vss(mat)*g0 0 0 0 Vsapc(mat)*g1 Vsapc(mat)*g2 Vsapc(mat)*g3 0 0 ; %|sa> 0 Esc(mat)/2 -Vpasc(mat)*conj(g1) -Vpasc(mat)*conj(g2) -Vpasc(mat)*conj(g3) 0 0 0 0 0 ; %|sc> 0 0 Epa(mat)/2 0 0 Vxx(mat)*g0 Vxy(mat)*g3 Vxy(mat)*g2 0 -Vpasec(mat)*g1; %|pxa> 0 0 0 Epa(mat)/2 0 Vxy(mat)*g3 Vxx(mat)*g0 Vxy(mat)*g1 0 -Vpasec(mat)*g2; %|pya> 0 0 0 0 Epa(mat)/2 Vxy(mat)*g2 Vxy(mat)*g1 Vxx(mat)*g0 0 -Vpasec(mat)*g3; %|pza> 0 0 0 0 0 Epc(mat)/2 0 0 Vseapc(mat)*g1 0 ; %|pxc> 0 0 0 0 0 0 Epc(mat)/2 0 Vseapc(mat)*g2 0 ; %|pyc> 0 0 0 0 0 0 0 Epc(mat)/2 Vseapc(mat)*g3 0 ; %|pzc> 0 0 0 0 0 0 0 0 Esea(mat)/2 0 ; %|sea> 0 0 0 0 0 0 0 0 0 Esec(mat)/2 ]; %|sec> H = [h+h' zeros(10); zeros(10) h+h']; % Compute spin orbital coupling contribution to Hamiltonian hso = zeros(20); hso(3,4) = -eye*soa(mat); hso(3,15) = soa(mat); hso(4,15) = -eye*soa(mat); hso(5,13) = -soa(mat); hso(5,14) = eye*soa(mat); hso(6,7) = -eye*soc(mat); hso(6,18) = soc(mat); hso(7,18) = -eye*soc(mat); hso(8,16) = -soc(mat); hso(8,17) = eye*soc(mat); hso(13,14) = eye*soa(mat); hso(16,17) = eye*soc(mat); Hso = hso + hso'; % Compute eigenvalues of spin-orbital modified Hamiltonian [V,D] = eig(H + Hso); eiglst = sum(D); E(i,:,d) = sort(real(eiglst)); end end %*************** % Plot Results * %*************** GLpt=-sqrt(3)/2;%distance from gamma to L point normalized to gamma-X GXpt=1; XUpt=0.5/sqrt(2); GKpt=0.5*3/sqrt(2); if direction == 1 % Reorder symmetry directions for plotting and create k-axis % L - Gamma symmetry for i = 1:N plotE(i,:) = E(i,:,1); plotK(i) = GLpt*(N-i)/(N-1); end % Gamma - X symmetry for i = 1:N plotE(N+i,:) = E(N+1-i,:,2); plotK(N+i) = GXpt*(i-1)/(N-1); end % X - U symmetry for i = 1:N plotE(2*N+i,:) = E(i,:,3); plotK(2*N+i) = GXpt + XUpt*(i-1)/(N-1); end % K - Gamma symmetry for i = 1:N plotE(3*N+i,:) = E(i,:,4); plotK(3*N+i) = (GXpt + XUpt) + GKpt*(i-1)/(N-1); end p0oset=-20; p1oset=-18;%offset on x-axis labels p2oset=-16; %******************** plot figure ************************ figure(1) clf plot(plotK,plotE,'LineWidth',LW) title(['Energy band diagram for ',matname,... ' along the L-\Gamma-X-U,K-\Gamma symmetry directions']) %text(0,p0oset,'Normalized Wavenumber, kL');%xlabel xlabel({'';'';'';'';'Normalized wave vector, \itkL\rm'}) ylabel('Energy, \itE\rm (eV)') xlim([GLpt (GXpt+XUpt+GKpt)]) grid on set(gca,'XTickLabel',{''}) text(GLpt*0.6,p1oset,'[111]'); text(GLpt,p2oset,'L');text(GLpt/2,p2oset,'\Lambda'); text(0,p2oset,'\Gamma') text(GXpt*0.4,p1oset,'[100]'); text(GXpt,p2oset,'X');text(GXpt/2,p2oset,'\Delta'); text((GXpt+XUpt),p2oset,'U,K') text((GXpt+XUpt+GKpt*0.41),p1oset,'[011]'); text((GXpt+XUpt+GKpt/2),p2oset,'\Sigma'); text((GXpt+XUpt+GKpt),p2oset,'\Gamma') set(gca,'XTick',[GLpt 0 GXpt (GXpt+XUpt) (GXpt+XUpt+GKpt)]) else dK = (kmax-kmin)/(N-1); plotK = kmin:dK:kmax; for i = 1:N plotE(i,:) = E(i,:,1); end plot(plotK,plotE) title(['Energy band diagram for ',matname,' in the [l,m,n]=[',... num2str(l),',',num2str(m),',',num2str(n),... '] direction with k_m_a_x=',num2str(kmax),... ' and k_m_i_n=',num2str(kmin),'.']) xlabel('Wave vector, \itkL\rm/(2\pi)') ylabel('Energy, \itE\rm (eV)') xlim([kmin kmax]) grid on end