function Hamiltonian = Double_Group_Hamiltonian(kx, ky, kz, x)
%function Hamiltonian = Double_Group_Hamiltonian(kx, ky, kz, x)
%
% Returns the block matrix Hamiltonian as a function of the three good
% translational symmetry Bloch vectors (kx ky kz).
%
% Hamiltonian may be 4x4, 6x6, 8x8, 15x15 or 30x30 depending on which final
% expression is uncommented.
%
% material parameters are for Si_{1-x}Ge_x alloys so select x = 0 for Silicon
% and x = 1 for Germanium.
%
% I have also configured the Hamiltonian for arbitrary growth orientations
% under a canonical transformation about angles (theta phi). A good test is
% for the [111] dispersion on a (001) Hamiltonian to match the (001) dispersion
% on a (111) Hamiltonian.
%
% This matlab function was written by Robert M. Ward (19/08/2011);
% robert.ward04@imperial.ac.uk
%
% Ph.D. Thesis - "Modelling of Silicon-Germanium Alloy Heterostructures
% using Double Group Formulation of k.p theory"
% https://spiral.imperial.ac.uk/handle/10044/1/9757
%Rydberg constant
Ryd = 13.6057;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%% This first part uses the single group formulation of k.p theory in a
%%%% 15x15 block matrix using the formulation of Manuel Cardona & Fred H. Pollak;
%%%% "Energy-Band Structure of Germanium & Silicon: The k.p Method"
%%%% Physical Review 142(2):530-543 Feb 1966
G1g1 = -12.700 - 0.180*x;
G5g1 = +00.000 + 0.000*x;
G2u1 = +04.150 - 3.260*x;
G4u1 = +03.335 - 0.222*x;
G1g2 = +08.400 - 1.600*x;
G3u1 = +08.540 + 1.760*x;
G5u2 = +11.700 - 0.340*x;
G2u2 = +15.800 - 1.800*x;
ZoneCenter = blkdiag(G1g1 ,eye(3)*G5g1, G2u1, eye(3)*G4u1, G1g2, eye(2)*G3u1, eye(3)*G5u2, G2u2);
P1 = +16.5990 - 0.4626*x;
P2 = -00.1088 + 1.0612*x - 0.6803*x^2;
P3 = +02.1225 - 0.1102*x;
P4 = +19.3881 - 0.3578*x;
Q1 = +14.5295 + 0.0925*x;
Q2 = -08.9185 - 1.4313*x;
R1 = +07.3838 + 1.2027*x;
R2 = +11.3499 - 0.1714*x;
T1 = +15.8642 - 0.3361*x - 0.5442*x^2;
T2 = +03.9457 + 0.1088*x;
Hk1 = [kx ky kz];
Hk2 = [ -kx -ky 2*kz
-sqrt(3)*kx sqrt(3)*ky 0 ];
Hk3 = [ 0 kz ky
kz 0 kx
ky kx 0 ];
%%%%%%%%%G1+1%%%%%%%%G5+1%%%%%%%%G2-1%%%%%%%%G4-1%%%%%%%%G1+2%%%%%%%%G3-1%%%%%%%%G5+2%%%%%%%%G2-2%%%%%%%%%%%%%%%%
kpHam = [0 zeros(1,3) 0 T2*Hk1 0 zeros(1,2) zeros(1,3) 0 %G1+1
zeros(3,1) zeros(3) P1*Hk1' Q1*Hk3 zeros(3,1) R1*Hk2' zeros(3) P3*Hk1' %G5+1
0 P1*Hk1 0 zeros(1,3) 0 zeros(1,2) P2*Hk1 0 %G2-1
T2*Hk1' Q1*Hk3 zeros(3,1) zeros(3) T1*Hk1' zeros(3,2) Q2*Hk3 zeros(3,1) %G4-1
0 zeros(1,3) 0 T1*Hk1 0 zeros(1,2) zeros(1,3) 0 %G1+2
zeros(2,1) R1*Hk2 zeros(2,1) zeros(2,3) zeros(2,1) zeros(2) R2*Hk2 zeros(2,1) %G3-1
zeros(3,1) zeros(3) P2*Hk1' Q2*Hk3 zeros(3,1) R2*Hk2' zeros(3) P4*Hk1' %G5+2
0 P3*Hk1 zeros(1,1) zeros(1,3) 0 zeros(1,2) P4*Hk1 0 ]; %G2-2
%%Adds hbar^2k^2/2m to the Hamiltonian
Hamiltonian = ZoneCenter + kpHam + Ryd*eye(length(ZoneCenter))*hypot(kx,hypot(ky,kz))^2; %#ok<*NASGU>
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%% This part uses the adapted double group formulation of k.p theory in a
%%%% 30x30 block matrix using the formulation of D. Rideau et al; "Strained
%%%% Si, Ge, and SiGe alloys modeled with a first principled optimized full
%%%% zone k.p method" Phyical Review B 74(195208):1-20 Nov 2006
Delta5g1 = 0.044 + 0.200*x + 0.0520*x^2;
Delta4u1 = 0.033 + 0.157*x;
Delta5g2 = 0.012 + 0.030*x;
HSO = [-1 +1i 0 0 0 +1
-1i -1 0 0 0 +1i
0 0 -1 -1 -1i 0
0 0 -1 -1 -1i 0
0 0 +1i +1i -1 0
+1 -1i 0 0 0 -1]/3;
ZoneCenter = kron(ZoneCenter,eye(2)) + blkdiag(zeros(2), HSO*Delta5g1, zeros(2), HSO*Delta4u1, zeros(6), HSO*Delta5g2, zeros(2));
Delta1 = (0.022 + 0.198*x)/Ryd;
Hk1 = kron(eye(2),Hk1);
Hk2 = kron(eye(2),Hk2);
Hk3 = kron(eye(2),Hk3);
%%%%%%%%%G1+1%%%%%%%%G5+1%%%%%%%%G2-1%%%%%%%%G4-1%%%%%%%%G1+2%%%%%%%%G3-1%%%%%%%%G5+2%%%%%%%%G2-2%%%%%%%%%%%%
kpHam = [zeros(2) zeros(2,6) zeros(2) T2*Hk1 zeros(2) zeros(2,4) zeros(2,6) zeros(2) %G1+1
zeros(6,2) zeros(6) P1*Hk1' Q1*Hk3 zeros(6,2) R1*Hk2' Delta1*HSO P3*Hk1' %G5+1
zeros(2) P1*Hk1 zeros(2) zeros(2,6) zeros(2) zeros(2,4) P2*Hk1 zeros(2) %G2-1
T2*Hk1' Q1*Hk3 zeros(6,2) zeros(6) T1*Hk1' zeros(6,4) Q2*Hk3 zeros(6,2) %G4-1
zeros(2) zeros(2,6) zeros(2) T1*Hk1 zeros(2) zeros(2,4) zeros(2,6) zeros(2) %G1+2
zeros(4,2) R1*Hk2 zeros(4,2) zeros(4,6) zeros(4,2) zeros(4) R2*Hk2 zeros(4,2) %G3-1
zeros(6,2) Delta1*HSO' P2*Hk1' Q2*Hk3 zeros(6,2) R2*Hk2' zeros(6) P4*Hk1' %G5+2
zeros(2) P3*Hk1 zeros(2) zeros(2,6) zeros(2) zeros(2,4) P4*Hk1 zeros(2) ]; %G2-2
%%Adds hbar^2k^2/2m to the Hamiltonian
Hamiltonian = ZoneCenter + kpHam + Ryd*eye(length(ZoneCenter))*hypot(kx,hypot(ky,kz))^2;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%% This part performs a unitary transformation of the 30x30 method of
%%%% Rideau using matrices based on Cracknell & Joshua "Symmetry Adapted
%%%% Functions for Double Point Groups ii. Cubic Point Groups" Proc. Camb.
%%%% Phil. Soc. 67:647-656 1970.
%Valence Band Eigenenergies; Top of Valence Band = 0eV
G6g1 = -12.700 - 0.180*x;
G7g1 = -00.044 - 0.200*x - 0.052*x^2;
G8g1 = +00.000 + 0.000;
%Conduction Band Eigenenergies
G7u1 = +04.1500 - 3.2600*x;
G6u1 = +03.3020 - 0.3790*x;
G8u1 = +03.3350 - 0.2220*x;
G6g2 = +08.4000 - 1.6000*x;
G8u2 = +08.5400 + 1.7600*x;
G7g2 = +11.6988 - 0.3700*x;
G8g2 = +11.7000 - 0.3400*x;
G7u2 = +15.8000 - 1.8000*x;
ZoneCenter = blkdiag(eye(2)*G6g1 ,eye(2)*G7g1, eye(4)*G8g1, eye(2)*G7u1, ...
eye(2)*G6u1, eye(4)*G8u1, eye(2)*G6g2, eye(4)*G8u2, ...
eye(2)*G7g2, eye(4)*G8g2, eye(2)*G7u2);
%G5+ x G6+ = G7+ + G8+ Basis transformation
CM1 = [ 0 0 -1i/sqrt(3) -1i/sqrt(3) +1/sqrt(3) 0
-1i/sqrt(3) -1/sqrt(3) 0 0 0 +1i/sqrt(3)
+1i/sqrt(6) +1/sqrt(6) 0 0 0 +1i*sqrt(2/3)
0 0 0 -1i/sqrt(2) -1/sqrt(2) 0
+1i/sqrt(2) -1/sqrt(2) 0 0 0 0
0 0 +1i*sqrt(2/3) -1i/sqrt(6) +1/sqrt(6) 0];
%G4- x G6+ = G7- + G8- Basis transformation
CM2 = [ 0 0 -1/sqrt(3) -1/sqrt(3) -1i/sqrt(3) 0
-1/sqrt(3) +1i/sqrt(3) 0 0 0 +1/sqrt(3)
-1/sqrt(2) -1i/sqrt(2) 0 0 0 0
0 0 +sqrt(2/3) -1/sqrt(6) -1i/sqrt(6) 0
+1/sqrt(6) -1i/sqrt(6) 0 0 0 +sqrt(2/3)
0 0 0 +1/sqrt(2) -1i/sqrt(2) 0];
%G3- x G6+ = G8- Basis transformation
CM3 = [0 0 -1 0
0 +1 0 0
0 0 0 -1
+1 0 0 0];
Hk1 = Hk1*CM1';
Hk2 = CM3*Hk2*CM1';
Hk3 = CM2*Hk3*CM1';
Hk4 = kron(eye(2),[kx ky kz])*CM2';
%%%%%%%%%G1+1xG6+%%%%G5+1xG6+%%%%G2-1xG6+%%%%G4-1xG6+%%%%G1+2xG6+%%%%G3-1xG6+%%%%G5+2xG6+%%%%G2-2xG6+%%%%%%%%
kpHam = [zeros(2) zeros(2,6) zeros(2) T2*Hk4 zeros(2) zeros(2,4) zeros(2,6) zeros(2) %G1+1xG6+
zeros(6,2) zeros(6) P1*Hk1' Q1*Hk3' zeros(6,2) R1*Hk2' zeros(6) P3*Hk1' %G5+1xG6+
zeros(2) P1*Hk1 zeros(2) zeros(2,6) zeros(2) zeros(2,4) P2*Hk1 zeros(2) %G2-1xG6+
T2*Hk4' Q1*Hk3 zeros(6,2) zeros(6) T1*Hk4' zeros(6,4) Q2*Hk3 zeros(6,2) %G4-1xG6+
zeros(2) zeros(2,6) zeros(2) T1*Hk4 zeros(2) zeros(2,4) zeros(2,6) zeros(2) %G1+2xG6+
zeros(4,2) R1*Hk2 zeros(4,2) zeros(4,6) zeros(4,2) zeros(4) R2*Hk2 zeros(4,2) %G3-1xG6+
zeros(6,2) zeros(6) P2*Hk1' Q2*Hk3' zeros(6,2) R2*Hk2' zeros(6) P4*Hk1' %G5+2xG6+
zeros(2) P3*Hk1 zeros(2) zeros(2,6) zeros(2) zeros(2,4) P4*Hk1 zeros(2) ]; %G2-2xG6+
%%Adds hbar^2k^2/2m to the Hamiltonian
Hamiltonian = ZoneCenter + kpHam + Ryd*eye(length(ZoneCenter))*hypot(kx,hypot(ky,kz))^2;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%% In this part I enter the k.p block matrices in accordance with the
%%%% double group formulation presented by Elder-Ward-Zhang "Double Group
%%%% Formulation of k.p theory for Cubic Crystals" Physical Review B,
%%%% 83(165210):1-23 Apr 2011
%%%% puts kx and ky into spherical basis for use in Wigner-Eckart theorem
kplus = (kx + 1i*ky);
kminus = (kx - 1i*ky);
%These are all the old parameters with a scaling of sqrt(6) for normalisation purposes;
P1 = +09.584 - 0.267*x;
P2 = -00.063 + 0.613*x - 0.393*x^2;
P3 = +01.226 - 0.064*x;
P4 = +11.194 - 0.207*x;
Q1 = -08.389 - 0.054*x;
Q2 = +05.149 + 0.826*x;
R1 = -04.263 - 0.694*x;
R2 = -06.553 + 0.099*x;
T1 = -09.159 + 0.194*x + 0.314*x^2;
T2 = -02.278 - 0.063*x;
%This is +A(k.pi)(G6-G6+)
Hk1 = [+kz +kplus
+kminus -kz];
%This is +A(k.pi)(G8-G6+)
Hk2 = -[-sqrt(3)*kplus 0
+2*kz -kplus
+kminus +2*kz
0 +sqrt(3)*kminus]/sqrt(2);
%This is +A(k.pi)(G8-G7+)
Hk3 = [+kminus +2*kz
0 -sqrt(3)*kminus
+sqrt(3)*kplus 0
+2*kz -kplus]/sqrt(2);
%This is A(k.pi,1)(G8-G8+)
Hk4 = [+kz 0 0 +kminus
0 -kz +kplus 0
0 +kminus +kz 0
+kplus 0 0 -kz ];
%This is A(k.pi,2)(G8-G8+)
Hk5 = [ 0 -sqrt(3)*kplus 0 +3*kminus
-sqrt(3)*kminus -4*kz +kplus 0
0 +kminus +4*kz -sqrt(3)*kplus
+3*kplus 0 -sqrt(3)*kminus 0 ]/sqrt(6);
%%%% This returns the Rotation matrices for appropriate Irreps used in the near set
[DG5, DG6, DG8] = Rotation_Matrices(0, 0); %[001] Hamiltonian
% [DG5, DG6, DG8] = Rotation_Matrices(0, pi/4); %[011] Hamiltonian
% [DG5, DG6, DG8] = Rotation_Matrices(pi/4, atan(sqrt(2))); %[111] Hamiltonian
%%%% This returns the matrices above using the Clebsch-Gordan procedure function
function Hk = Linear_Matrices(m1, m2)
A1 = zeros(length(m1),length(m2));
A2 = zeros(length(m1),length(m2));
A3 = zeros(length(m1),length(m2));
for p = 1:length(m1)
for q = 1:length(m2)
A1(p,q) = +Clebsch_Gordan(1,-1,max(m1),m1(p),max(m2),m2(q));
A2(p,q) = -Clebsch_Gordan(1,+0,max(m1),m1(p),max(m2),m2(q));
A3(p,q) = +Clebsch_Gordan(1,+1,max(m1),m1(p),max(m2),m2(q));
end
end
Hk = (DG5(1,1)*A1 + DG5(2,1)*A2 + DG5(3,1)*A3)*kplus/sqrt(2) - ...
(DG5(1,2)*A1 + DG5(2,2)*A2 + DG5(3,2)*A3)*kz - ...
(DG5(1,3)*A1 + DG5(2,3)*A2 + DG5(3,3)*A3)*kminus/sqrt(2);
end
Hk1 = -sqrt(3)*Linear_Matrices(+1/2:-1:-1/2, +1/2:-1:-1/2);
Hk2 = +sqrt(6)*Linear_Matrices(+3/2:-1:-3/2, +1/2:-1:-1/2);
U = [ -sqrt(1/6) 0 0 0 +sqrt(5/6) 0
0 +sqrt(5/6) 0 0 0 -sqrt(1/6)
0 -sqrt(1/6) 0 0 0 -sqrt(5/6)
0 0 +1 0 0 0
0 0 0 -1 0 0
+sqrt(5/6) 0 0 0 +sqrt(1/6) 0 ];
Hkb = -sqrt(6)*Linear_Matrices(+3/2:-1:-3/2, +5/2:-1:-5/2)*U';
Hk3 = -Hkb(1:4,1:2);
Hka = -Linear_Matrices(+3/2:-1:-3/2, +3/2:-1:-3/2);
Hk4 = Hkb(1:4,3:6)*sqrt(2/5) + Hka*sqrt(3/5);
Hk5 = Hkb(1:4,3:6)*sqrt(3/5) - Hka*sqrt(2/5);
%%%% Try switching the comments and see if there is a difference???
%%%% The T-matrix was presented in the publication by Elder-Ward-Zhang
T = [ [ 0 0 +1i 0 ]
[ 0 0 0 -1i ]
[ -1i 0 0 0 ]
[ 0 +1i 0 0 ] ];
%Hk2 = -Hkb(1:4,1:2);
%Hk3 = +sqrt(6)*Linear_Matrices(+3/2:-1:-3/2, +1/2:-1:-1/2);
%Hk5 = T*(Hkb(1:4,3:6)*sqrt(3/5) - Hka*sqrt(2/5))*T';
%%%%%%%%%%%%G6+1%%%%%%%%G7+1%%%%%%%%%%%%G8+1%%%%%%%%%%%%%G7-1%%%%%%%%G6-1%%%%%%%%G8-1%%%%%%%%%G6+2%%%%%%%%G8-2%%%%%%%%%%%%%G7+2%%%%%%%%%%%%G8+2%%%%%%%%%%%%%G7-2%%%%%%%%%%%%%%%%%%
kpHam = [ [ zeros(2) zeros(2) zeros(2,4) zeros(2) T2*Hk1 T2*Hk3' zeros(2) zeros(2,4) zeros(2) zeros(2,4) zeros(2) ] %G6+1
[ zeros(2) zeros(2) zeros(2,4) -P1*Hk1 zeros(2) -Q1*Hk2' zeros(2) R1*sqrt(2)*Hk2' zeros(2) zeros(2,4) -P3*Hk1 ] %G7+1
[ zeros(4,2) zeros(4,2) zeros(4) P1*Hk2 Q1*Hk3 Q1*Hk4 zeros(4,2) R1*sqrt(3)*Hk5 zeros(4,2) zeros(4) P3*Hk2 ] %G8+1
[ zeros(2) -P1*Hk1' P1*Hk2' zeros(2) zeros(2) zeros(2,4) zeros(2) zeros(2,4) -P2*Hk1 P2*Hk2' zeros(2) ] %G7-1
[ T2*Hk1' zeros(2) Q1*Hk3' zeros(2) zeros(2) zeros(2,4) T1*Hk1 zeros(2,4) zeros(2) Q2*Hk3' zeros(2) ] %G6-1
[ T2*Hk3 -Q1*Hk2 Q1*Hk4' zeros(4,2) zeros(4,2) zeros(4) T1*Hk3 zeros(4) -Q2*Hk2 Q2*Hk4 zeros(4,2) ] %G8-1
[ zeros(2) zeros(2) zeros(2,4) zeros(2) T1*Hk1' T1*Hk3' zeros(2) zeros(2,4) zeros(2) zeros(2,4) zeros(2) ] %G6+2
[ zeros(4,2) R1*sqrt(2)*Hk2 R1*sqrt(3)*Hk5' zeros(4,2) zeros(4,2) zeros(4) zeros(4,2) zeros(4) R2*sqrt(2)*Hk2 R2*sqrt(3)*Hk5 zeros(4,2) ] %G8-2
[ zeros(2) zeros(2) zeros(2,4) -P2*Hk1' zeros(2) -Q2*Hk2' zeros(2) R2*sqrt(2)*Hk2' zeros(2) zeros(2,4) -P4*Hk1 ] %G7+2
[ zeros(4,2) zeros(4,2) zeros(4) P2*Hk2 Q2*Hk3 Q2*Hk4' zeros(4,2) R2*sqrt(3)*Hk5' zeros(4,2) zeros(4) P4*Hk2 ] %G8+2
[ zeros(2) -P3*Hk1' P3*Hk2' zeros(2) zeros(2) zeros(2,4) zeros(2) zeros(2,4) -P4*Hk1' P4*Hk2' zeros(2) ] ]; %G7-2
%%Adds hbar^2k^2/2m to the Hamiltonian
Hamiltonian = ZoneCenter + kpHam + Ryd*eye(length(ZoneCenter))*hypot(kx,hypot(ky,kz))^2;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%% In this part I use the Lowdin renormalisation procedure of second
%%%% order degenerate perturbation theory to generate 6-band k.p effective
%%%% mass Hamiltonians. Per-Olov-Lowdin "A note on the quantum mechanical
%%%% perturbation theory" The Journal of Chemical Physics, 19(11):1396-1401
%%%% Nov 1951
%%%%%%%%%%%%%%%%%%%%%%G7-1 + G7-2%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%G6-1%%%%%%%%%%%%%%%%%%%%%%G8-1%%%%%%%%%%%%%%%%%%%%%%%%G8-2%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Projector = [ DG6'*[ -sqrt(P1^2/(G7u1 - G7g1) + P3^2/(G7u2 - G7g1))*Hk1 zeros(2) -Q1*Hk2'/sqrt(G8u1 - G7g1) +R1*sqrt(2)*Hk2'/sqrt(G8u2 - G7g1) ] %G7+1
DG8'*[ +sqrt(P1^2/(G7u1 - G8g1) + P3^2/(G7u2 - G8g1))*Hk2 Q1*Hk3/sqrt(G6u1 - G8g1) Q1*Hk4/sqrt(G8u1 - G8g1) +R1*sqrt(3)*Hk5/sqrt(G8u2 - G8g1) ] ]; %G8+1
%Adds hbar^2k^2/2m to the Hamiltonian
%Hamiltonian = ZoneCenter(3:8,3:8) - (Projector*Projector') + Ryd*eye(6)*hypot(kx,hypot(ky,kz))^2;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%% This is what one returns from using second order degenerate
%%%% perturbation theory. I have expressed the scaling parameters in terms
%%%% of Luttinger invariants J. M. Luttinger "Quantum Theory of Cyclotron
%%%% Resonance in Semiconductors: General Theory" Physical Review,
%%%% 102(4):1030-1041 May 1956
%%%% One should compare expressions for these with the Foreman parameters,
%%%% (sigma delta pi) given in Bradley A. Foreman "Effective-mass
%%%% Hamiltonian and Boundary conditions for the valence bands of
%%%% semiconductor microstructures" Physical Review B Rapid Communications,
%%%% 48(7):4964-4969, Aug 1993
gamma1 = P1^2/(G7u1 - G8g1) + P3^2/(G7u2 - G8g1) + Q1^2*(1/(G6u1 - G8g1) + 1/(G8u1 - G8g1)) + 4*R1^2/(G8u2 - G8g1) - Ryd;
gamma2 = (+1/2)*(P1^2/(G7u1 - G8g1) + P3^2/(G7u2 - G8g1) - Q1^2/(G6u1 - G8g1) + 4*R1^2/(G8u2 - G8g1));
gamma3 = (+1/2)*(P1^2/(G7u1 - G8g1) + P3^2/(G7u2 - G8g1) + Q1^2/(G6u1 - G8g1) - 2*R1^2/(G8u2 - G8g1));
mu = (gamma3 - gamma2)/2;
gamma = (gamma3 + gamma2)/2;
gamma1pr = P1^2/(G7u1 - G7g1) + P3^2/(G7u2 - G7g1) + 2*Q1^2/(G8u1 - G7g1) + 4*R1^2/(G8u2 - G7g1) - Ryd;
gamma2pr = (+1/2)*(P1^2/sqrt((G7u1 - G7g1)*(G7u1 - G8g1)) + ...
P3^2/sqrt((G7u2 - G7g1)*(G7u2 - G8g1)) - ...
Q1^2/sqrt((G8u1 - G7g1)*(G8u1 - G8g1)) + ...
4*R1^2/sqrt((G8u2 - G7g1)*(G8u2 - G8g1)));
gamma3pr = (+1/2)*(P1^2/sqrt((G7u1 - G7g1)*(G7u1 - G8g1)) + ...
P3^2/sqrt((G7u2 - G7g1)*(G7u2 - G8g1)) + ...
Q1^2/sqrt((G8u1 - G7g1)*(G8u1 - G8g1)) - ...
2*R1^2/sqrt((G8u2 - G7g1)*(G8u2 - G8g1)));
mupr = (gamma3pr - gamma2pr)/2;
gammapr = (gamma3pr + gamma2pr)/2;
%%%% These are expressions for different growth directions under a
%%%% canonical transformation. One should consult Jian-Bai Xia "Effective
%%%% mass Theory for supperlattices grown on (11N)-oriented substrates"
%%%% Physical Review B 43(12):9856-9864 Apr 1991
%%%%%%%%%%% [001] Terms
P = +gamma1*(kplus*kminus + kz^2);
Ppr = +gamma1pr*(kplus*kminus + kz^2);
Q = +gamma2*(kplus*kminus - 2*kz^2);
Qpr = +gamma2pr*(kplus*kminus - 2*kz^2);
R = +sqrt(3)*(mu*kminus^2 - gamma*kplus^2);
Rpr = +sqrt(3)*(mupr*kminus^2 - gammapr*kplus^2);
S = -2*sqrt(3)*gamma3*kz*kplus;
Spr = -2*sqrt(3)*gamma3pr*kz*kplus;
% %%%%%%%%%% [011] Terms
% P = +gamma1*(kplus*kminus + kz^2);
% Ppr = +gamma1pr*(kplus*kminus + kz^2);
%
% Q = +gamma2*kx^2 - (1/2)*(gamma2 - 3*gamma3)*ky^2 - (1/2)*(gamma2 + 3*gamma3)*kz^2;
% Qpr = +gamma2pr*kx^2 - (1/2)*(gamma2pr - 3*gamma3pr)*ky^2 - (1/2)*(gamma2pr + 3*gamma3pr)*kz^2;
%
% R = -sqrt(3)*(gamma2*kx^2 - gamma*ky^2 + 2i*gamma3*kx*ky + mu*kz^2);
% Rpr = -sqrt(3)*(gamma2pr*kx^2 - gammapr*ky^2 + 2i*gamma3pr*kx*ky + mupr*kz^2);
%
% S = -2*sqrt(3)*kz*(gamma2*ky + 1i*gamma3*kx);
% Spr = -2*sqrt(3)*kz*(gamma2pr*ky + 1i*gamma3pr*kx);
% %%%%%%%%%%% [111] Terms
% P = +gamma1*(kplus*kminus + kz^2);
% Ppr = +gamma1pr*(kplus*kminus + kz^2);
%
% Q = +gamma3*(kplus*kminus - 2*kz^2);
% Qpr = +gamma3pr*(kplus*kminus - 2*kz^2);
%
% R = +(2/sqrt(3))*( 2i*sqrt(2)*mu*kminus*kz + (gamma + (1/2)*gamma3)*kminus^2);
% Rpr = +(2/sqrt(3))*( 2i*sqrt(2)*mupr*kminus*kz + (gammapr + (1/2)*gamma3pr)*kminus^2);
%
% S = -(2/sqrt(3))*((gamma2 + 2*gamma)*kplus*kz - 1i*sqrt(2)*mu*kminus^2);
% Spr = -(2/sqrt(3))*((gamma2pr + 2*gammapr)*kplus*kz - 1i*sqrt(2)*mupr*kminus^2);
HSOSO = [Ppr 0
0 Ppr ];
HVBVB = [ P + Q +S +R 0
+conj(S) P - Q 0 +R
+conj(R) 0 P - Q -S
0 +conj(R) -conj(S) P + Q ];
HSOVB = [ +Spr/2 +Rpr
-Qpr -sqrt(3/4)*Spr
-sqrt(3/4)*conj(Spr) +Qpr
-conj(Rpr) +conj(Spr)/2 ]*sqrt(2);
Lowdin = [HSOSO HSOVB'
HSOVB HVBVB];
%Hamiltonian = ZoneCenter(3:8,3:8) - Lowdin;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%% In this part I include the conduction band minima into the near set.
%%%% Note that this is only strictly valid for high Germanium compositions
%%%% when the direct gap minima is G7- symmetry. This is similar to the
%%%% method given by Thomas B. Bahder; "Eight-Band k.p model of strained
%%%% zinc blende crystals" Physical Review B, 41(17):11992-12001, Jun 1990
%%%%%%%%%%G7+1%%%%%%%%%%%%%%%G8+1%%%%%%%%%%%%%%%G7-1%%%%%%%%%%
kpHam = [ zeros(2) zeros(2,4) -P1*DG6'*Hk1*DG6 %G7+1
zeros(4,2) zeros(4) +P1*DG8'*Hk2*DG6 %G8+1
-P1*DG6'*Hk1'*DG6 +P1*DG6'*Hk2'*DG8 zeros(2) ]; %G7-1
%%%%%%%%%%%%%%%%%%%%%%G6-1%%%%%%%%%%%%%%%%%%%%%%%G8-1%%%%%%%%%%%%%%%%%%%%%%%%G8-2%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%G7+2%%%%%%%%%%%%%%%%%%%%%%G8+2%%%%%%%%%%%%%%%%%%%%%%G7-2%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Projector = [ DG6'*[ zeros(2) -Q1*Hk2'/sqrt(G8u1 - G7g1) R1*sqrt(2)*Hk2'/sqrt(G8u2 - G7g1) zeros(2) zeros(2,4) -P3*Hk1/sqrt(G7u2 - G7g1) ] %G7+1
DG8'*[ Q1*Hk3/sqrt(G6u1 - G8g1) Q1*Hk4/sqrt(G8u1 - G8g1) R1*sqrt(3)*Hk5/sqrt(G8u2 - G8g1) zeros(4,2) zeros(4) P3*Hk2/sqrt(G7u2 - G8g1) ] %G8+1
DG6'*[ zeros(2) zeros(2,4) zeros(2,4) -P2*Hk1/sqrt(G7g2 - G7u1) P2*Hk2'/sqrt(G8g2 - G7u1) zeros(2) ] ]; %G7-1
%%Adds hbar^2k^2/2m to the Hamiltonian
%Hamiltonian = ZoneCenter(3:10,3:10) + kpHam - (Projector*Projector') + Ryd*eye(8)*hypot(kx,hypot(ky,kz))^2;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
gamma1 = gamma1 - P1^2/(G7u1 - G8g1);
gamma2 = gamma2 - (+1/2)*P1^2/(G7u1 - G8g1);
gamma3 = gamma3 - (+1/2)*P1^2/(G7u1 - G8g1);
mu = (gamma3 - gamma2)/2;
gamma = (gamma3 + gamma2)/2;
gamma1pr = gamma1pr - P1^2/(G7u1 - G7g1);
gamma2pr = gamma2pr - (+1/2)*P1^2/sqrt((G7u1 - G7g1)*(G7u1 - G8g1));
gamma3pr = gamma3pr - (+1/2)*P1^2/sqrt((G7u1 - G7g1)*(G7u1 - G8g1));
mupr = (gamma3pr - gamma2pr)/2;
gammapr = (gamma3pr + gamma2pr)/2;
me = P2^2*(1/(G7g2 - G7u1) + 2/(G8g2 - G7u1)) - Ryd;
%%%%%%%%%%% [001] Terms
P = +gamma1*(kplus*kminus + kz^2);
Ppr = +gamma1pr*(kplus*kminus + kz^2);
Q = +gamma2*(kplus*kminus - 2*kz^2);
Qpr = +gamma2pr*(kplus*kminus - 2*kz^2);
R = +sqrt(3)*(mu*kminus^2 - gamma*kplus^2);
Rpr = +sqrt(3)*(mupr*kminus^2 - gammapr*kplus^2);
S = -2*sqrt(3)*gamma3*kz*kplus;
Spr = -2*sqrt(3)*gamma3pr*kz*kplus;
A = +me*(kplus*kminus + kz^2);
% %%%%%%%%%% [011] Terms
% P = +gamma1*(kplus*kminus + kz^2);
% Ppr = +gamma1pr*(kplus*kminus + kz^2);
%
% Q = +gamma2*kx^2 - (1/2)*(gamma2 - 3*gamma3)*ky^2 - (1/2)*(gamma2 + 3*gamma3)*kz^2;
% Qpr = +gamma2pr*kx^2 - (1/2)*(gamma2pr - 3*gamma3pr)*ky^2 - (1/2)*(gamma2pr + 3*gamma3pr)*kz^2;
%
% R = -sqrt(3)*(gamma2*kx^2 - gamma*ky^2 + 2i*gamma3*kx*ky + mu*kz^2);
% Rpr = -sqrt(3)*(gamma2pr*kx^2 - gammapr*ky^2 + 2i*gamma3pr*kx*ky + mupr*kz^2);
%
% S = -2*sqrt(3)*kz*(gamma2*ky + 1i*gamma3*kx);
% Spr = -2*sqrt(3)*kz*(gamma2pr*ky + 1i*gamma3pr*kx);
%
% A = +me*(kplus*kminus + kz^2);
% %%%%%%%%%%% [111] Terms
% P = +gamma1*(kplus*kminus + kz^2);
% Ppr = +gamma1pr*(kplus*kminus + kz^2);
%
% Q = +gamma3*(kplus*kminus - 2*kz^2);
% Qpr = +gamma3pr*(kplus*kminus - 2*kz^2);
%
% R = +(2/sqrt(3))*( 2i*sqrt(2)*mu*kminus*kz + (gamma + (1/2)*gamma3)*kminus^2);
% Rpr = +(2/sqrt(3))*( 2i*sqrt(2)*mupr*kminus*kz + (gammapr + (1/2)*gamma3pr)*kminus^2);
%
% S = -(2/sqrt(3))*((gamma2 + 2*gamma)*kplus*kz - 1i*sqrt(2)*mu*kminus^2);
% Spr = -(2/sqrt(3))*((gamma2pr + 2*gammapr)*kplus*kz - 1i*sqrt(2)*mupr*kminus^2);
%
% A = +me*(kplus*kminus + kz^2);
HSOSO = [Ppr 0
0 Ppr ];
HVBVB = [ P + Q +S R 0
+conj(S) P - Q 0 R
+conj(R) 0 P - Q -S
0 +conj(R) -conj(S) P + Q ];
HSOVB = [ +Spr/2 +Rpr
-Qpr -sqrt(3/4)*Spr
-sqrt(3/4)*conj(Spr) +Qpr
-conj(Rpr) +conj(Spr)/2 ]*sqrt(2);
HCBCB = [ A 0
0 A ];
HSOCB = zeros(2);
HVBCB = zeros(2,4);
Lowdin = [HSOSO HSOVB' HSOCB'
HSOVB HVBVB HVBCB'
HSOCB HVBCB HCBCB];
%Hamiltonian = ZoneCenter(3:10,3:10) + kpHam - Lowdin;
end
%%%%robert.ward04@imperial.ac.uk 19/08/2011