AlN/GaNデバイスのバンド構造の計算

AlGaN系のバンド計算には、ひずみや分極が組み込まれたBandEngを使うと楽ですが、材料によっては自力で計算する必要がある場合があります。Mathematicaを使って自力でSiのバンド構造の簡易計算をする場合は、下記のようになります。

(*constants*)
q = 1.60218*10^-19; (*elementary charge[C]*)
h = 6.6262*10^-34; (*plank constant[m^2 kg/s]*)
m = 9.1095*10^-31; (*electron weight[kg]*)
d0 = 0.5*q; (*valence band[J]*)
Eg = 1.0*q; (*band gap [J]*)
P = 9*10^-8*q*10^2; (*kp band diagram [J m]*)

hk[k_] := (2*Pi*h*k)^2/(2*0.33*m);

ec[k_] := Eg + hk[k] + P^3*k^2/3*(2/Eg + 1/(Eg + d0)); (*conduction band edge*)
ev1[k_] := hk[k]; (*heavy hole band edge*)
ev2[k_] := hk[k] - 2 P^2 k^2/ (3 Eg); (*light hole band edge*)
ev3[k_] := hk[k] - d0 - P^2 k^2/(3 (Eg + d0)); (*spin split band egde*)
Plot[{ec[x]/q, ev1[x]/q, ev2[x]/q, ev3[x]/q}, {x, -3 10^5,
  3 10^5(*m-1*)}, PlotRange -> {-1 , 1.1}, PlotLabel -> "Energy [eV]"]


上記の計算を発展させて、Mathematicaを使ってGaN/AlNのバンド構造を計算すると下記のようになります。ここまで長くなると、計算が合っているのか不安になってきます。あくまで参考程度でお願いします。

(*constants*)
q = 1.60218*10^-19 (*elementary charge[C]*);
k = 1.38066*10^-23 (*Boltzman constant[J/K]*);
T = 300 (*Absolute temperature[K]*);
h = 6.6262 10 ^-34(*plank constant[Js]*);
me = 0.22 9.10938 10^-31(*electron weight[kg]*);
e0 = 8.85418*10^-14 (*permitivity in vaccum[F/cm]*);
eg = 9.5 e0 (*GaN permitivity*);
ea = 9.0 e0 (*AlN permitivity*);
Nc = 2 (2 Pi me k T/h^2)^1.5*10^-6  (* effective density of state in conduction band[cm-3]*);

(*input parameters*)
tg = 2 10^-6 (*GaN thickness[cm]*);
eb = 0.84(*GaN schottky barrier[eV]*);
dEc = 0.7 (6.05 - 3.42)(*band offset[eV]*);
Cg = eg/tg(*GaN capacitance*);
Ng = 3 10^18(*GaN donor concentration [cm-3]*); (*substrate property*)
ta = 10^-6(*AlN thickness*);
ni = 9.4 10^-34(*intrinsic carrier concentration[cm-3]*);
n0 = 1*10^13 (*AlN donor concentration[cm^-3]*);
p0 = 10 ^5(*ni^2/n0 ? AlN acceptor concentration[cm-3]*);
LD = Sqrt[(k T ea)/(n0 q^2)];(*AlN Debye length*)
Ec = k T/q Log[Nc/n0];(*AlN conduction band edge*)
Print["Nc=", N[Nc], " [cm-3]\tEc=", Ec, " [eV]\tDebye length=", LD, "[cm]"];

(*Semiconductor charge*)
F[x_, y_] := Sign[x] Sqrt[Exp[-x] + x - 1 + y (Exp[x] - x - 1)];
Ex[p_, n0_] := Sqrt[2] k T/q/LD F[q p/(k T), n0/p0];
Qs[p_, n0_] := -ea Ex[p, n0]; 

(*LogPlot[Abs[Qs[p,n0]],{p,-0.5,1.5}]*)

(*GaN Energy profile*)
pg[x_ /; x <= 0] := q Ng (x + tg)^2/(2 eg);
pg[x_ /; x > 0] := eb;
gan = Plot[{eb - pg[x]}, {x, -tg, ta}, PlotRange -> {-0.5, 1}];

(*AlN Energy profile*)
ps = eb - pg[0] + dEc (* [V]*);
eval = NDSolve[{p'[x] + Ex[p[x], n0] == 0, p[0] == ps}, p, {x, 0, ta}];
aln = Plot[{Evaluate[Ec - p[x] /. eval]}, {x, 0, ta}, PlotRange -> {-0.5, 1}];
Show[gan, aln]

(*strain*)
aa = 3.112 10^-10(*lattice constant of AlN*);
ag = 3.189 10^-10(*lattice constant of GaN*);
Psa = -0.081 (*Spontaneous polarization of AlN*);
Psg = -0.029(*Spontaneous polarization of GaN*);
e31a = -0.60(*Dielectric constants of AlN*);
e31g = -0.49(*Dielectric constants of GaN*);
e33a = 1.46(*Dielectric constants of AlN*);
e33g = 0.73(*Dielectric constants of GaN*);
ae[x_] := (aa - ag) x + ag(*AlxGa1-xN lattice constant[m]*);
e31[x_] := (e31a - e31g) x + e31g;(*elastic constant[C m-2]*)
e33[x_] := (e33a - e33g) x + e33g;(*elastic constant[C m-2]*)
c13[x_] := 5 x + 103;(*[GPa]*)
c33[x_] := -32 x + 405;(*[GPa]*)
eb1 = 2.5(*AlN schottky barrier[eV]*);
ppe[x_, y_] := 2 (ae[y] - ae[x])/ae[x] (e31[x] - e33[x] c13[x]/c33[x]) 10^-4;

(*piezo polarization of AlxGa1-xN*)

psp[x_] := ((Psa - Psg) x + Psg) 10^-4;(*spontaneous polarization of AlxGa1-xN*)
qs[x_, y_] := ppe[x, y] + psp[x] 

(*Ppe (AlxGaN)+Psp (AlxGaN:strain)-Psp (AlyGaN:relaxed)*);
Print["Sheet charge density |qs| (=Ppe(GaN)+Psp(GaN)-Psp(AlN)):", Abs[qs[0, 1]/q] "[cm-2]"];

(*AlN cap*)
pa1[x_] := (psp[1] - psp[0]) (x + 2 tg)/ea  ;
aln1 = Plot[{eb1 - 0.1 pa1[z]}, {z, -2 tg, -tg}];

(*GaN carrier*)
pg1[x_] := (psp[0] + ppe[0, 1] - psp[1]) /eg (x + tg);
gan1 = Plot[{eb1 - 0.1 pa1[-tg] - dEc - 23 - pg1[x]}, {x, -tg, 0}];

(*AlN sub*)
ps1 = eb1 - 0.1 pa1[-tg] - dEc - 23 - pg1[0] (* [V]*);
eval1 = NDSolve[{pp'[x] + Ex[pp[x], n0] == 0, pp[ta] == 0}, pp, {x, 0, ta}];
aln2 = Plot[{Evaluate[Ec + ps1 - pp[y] /. eval1]}, {y, 0, ta}];
Show[gan1, aln2, aln1, PlotRange -> {-0.5, 3}, AxesOrigin -> {0, 0}]
(*pa1の0 .1倍とgan1の-23の意味？表面電荷は？２DEGは？*)