DJ’s Mathematica File for Quantum Physics of Semiconductors and Nanostructures
Debdeep Jena (djena@cornell.edu)

Module I: Fundamentals

1. And off we go! [includes Basics of Mathematica]

Basics of Mathematica I

q=1.6*(10^(-19))  (* Electron charge in Coulombs, how to Define a constant and comment *)

f[x_]:=(x^2)-x+1;  (* How to define a function *)

Solve[f[x]==0,x]  (* How to solve a polynomial equation *)

Solve[f[x]==0,x][[1]][[1]][[2]] (* How to get the 1st solution of the equation *)

Solve[f[x]==0,x][[2]][[1]][[2]] (* How to get the 2nd solution of the equation *)

Plot[{x,Cos[x]},{x,0,2*π},Frame->True,FrameLabel->{"x","x and cos[x]"}, BaseStyle->{FontSize->15}]  (* How to plot and label *)

FindRoot[Cos[x]==x,{x,0}]  (* How to find numerical roots of difficult equations *)
FindRoot[Cos[x]==x,{x,0}][[1]][[2]]  (* How to get the numerical value of the root *)

∫(Sqrt[Tan[x]])\[DifferentialD]x  (* How to evaluate an indefinite integral analytically *)
(∫_0^∞)(x^2)*Exp[-((((x/a)))^2)]\[DifferentialD]x  (* How to evaluate a definite integral analytically *)
NIntegrate[(x^2)*Exp[-(x^2)],{x,0,1000}] (* How to evaluate any integral numerically *)

Basics of Mathematica 2

A={{1,2},{0,1}}  (* How to define a simple matrix *)
MatrixForm[A]  (* How to display a matrix in the normal format *)

B={{2,1},{1,2}}  (* How to define a simple matrix *)
MatrixForm[A.B] (* Product of two matrices *)
MatrixForm[B.A] (* Product of two matrices; note A.B≠B.A *)

MatrixForm[MatrixPower[A,20]]  (* How to take powers of matrices *)
MatrixForm[N[MatrixFunction[Sin,A]]]  (* How to take numerical functions of matrices *)

Eigenvalues[A]  (* How to find eigenvalues of matrices *)
Eigenvectors[A]  (* How to find eigenvectors of matrices *)

MatrixForm[IdentityMatrix[10]] (* How to display a 10x10 identity matrix *)

tightbinding[center_,left_,right_]=center*IdentityMatrix[10]+SparseArray[{Band[{1,2}]->left,Band[{2,1}]->right},{10,10}] (* How to create a 10x10 tight-binding matrix *)
MatrixForm[tightbinding[1.0,-0.1,-0.1]] (* How to display the 10x10 tight-binding matrix *)
MatrixPlot[tightbinding[1.0,-0.1,-0.1]] (* How to "Plot" the 10x10 tight-binding matrix *)
ListPlot[Eigenvalues[tightbinding[1.0,-0.1,-0.1]]]  (* How to Plot the eigenvalues of the 10x10 tight-binding matrix *)

tightbindingBig[center_,left_,right_]=center*IdentityMatrix[100]+SparseArray[{Band[{1,2}]->left,Band[{2,1}]->right},{100,100}] ;(* How to create a 100x100 tight-binding matrix *)
Manipulate[Quiet[ListPlot[Eigenvalues[tightbindingBig[1.0,left,left]],PlotRange->{0,2},Frame->True,FrameLabel->{"State number","Eigenvalue"},PlotStyle->{Dashed,Thick},BaseStyle->{FontSize->15},AspectRatio->GoldenRatio,GridLines->Automatic]]  ,{left,0,-0.3}](* How to Plot the eigenvalues of the 100x100 tight-binding matrix with a tunable hopping parameter.  This is exactly how electrons develop a bandstructure in solids! *)

tightbindingBigger[center_,left1_,right1_,left2_,right2_]=center*IdentityMatrix[200]+SparseArray[{Band[{1,2}]->left1,Band[{2,1}]->right1},{200,200}]+SparseArray[{Band[{1,3}]->left2,Band[{3,1}]->right2},{200,200}] ;;
(* How to create a 200x200 tight-binding matrix *)
Manipulate[Quiet[ListPlot[Eigenvalues[tightbindingBigger[1.0,left1,left1,left2,left2]],PlotRange->{0,2},Frame->True,FrameLabel->{"State number","Eigenvalue"},PlotStyle->{Dashed,Thick},BaseStyle->{FontSize->15},AspectRatio->GoldenRatio,GridLines->Automatic]]  ,{left1,0,-0.3},{left2,0,-0.3}](* How to Plot the eigenvalues of the 200x200 tight-binding matrix with a tunable hopping parameter.  This is exactly how electrons develop a bandstructure in solids! *)

2. Secrets of the Classical Electron

Optical conductivity of metals

sigma[ω_,τ_]:=(τ/(1-I*ω*τ));
Manipulate[LogLinearPlot[{Re[sigma[ω,τ]],Im[
sigma[ω,τ]]},{ω,0.001,11},PlotRange->{0,15},Frame->True,PlotStyle->Thick,FrameLabel->{"ω ((10^12) /s)","σ(ω) (a.u.)"},BaseStyle->{FontSize->15}],{τ,10,0.1}]

3. Quantum Mechanics in a Nutshell

Planck’s resolution of the Ultraviolet Catastrophe

ℏ:=((6.63*(10^(-34)))/(2π));  (* Planck's constant *)
c:=3*(10^8); (* Speed of Light *)
kb:=1.38*(10^(-23)); (* Boltzmann's constant *)
q:=1.6*(10^(-19)); (* Coulomb, Electron Charge *)
u[ω_,T_]:=(1/(q*(π^2)*(ℏ^2)*(c^3)))*((((ℏ*ω*(10^12)))^3)/(Exp[((ℏ*(ω*(10^12)))/(kb*T))]-1)); (* Spectral Density of Blackbody Radiation *)
uλ[λ_,T_]:=(1/(q*(π^2)*(ℏ^2)*(c^3)))*((((ℏ*((2*π*c)/(λ*(10^(-9))))))^3)/(Exp[((ℏ*(((2*π*c)/(λ*(10^(-9))))))/(kb*T))]-1)); (* Spectral Density of Blackbody Radiation *)
Manipulate[Plot[{u[ω,T],u[ω,3000],u[ω,4500],u[ω,5800]},{ω,0,6000},AspectRatio->GoldenRatio,Frame->True,FrameLabel->{"ω ((10^12) /s)","u[ω] (eV.s/(m^3))"},PlotStyle->{{Gray,Thick},{Blue,Thick},{Green,Thick},{Red,Thick}},BaseStyle->{FontSize->15},PlotRange->{-1,1800}],{T,4,6000}]
Manipulate[Quiet[Plot[{uλ[λ,T],uλ[λ,3000],uλ[λ,4500],uλ[λ,5800]},{λ,0,3000},AspectRatio->GoldenRatio,Frame->True,FrameLabel->{"λ (nm)","u[ω] (eV.s/(m^3))"},PlotStyle->{{Gray,Thick},{Blue,Thick},{Green,Thick},{Red,Thick}},BaseStyle->{FontSize->15},PlotRange->{-1,1800}]],{T,4,6000}]

Hydrogen Atom

λ[n1_,n2_]:=(1240/(13.6((1/(n1^2))-(1/(n2^2)))));  (* Hydrogen atom spectral lines in nm *)
sharpfunction[x_,x0_,Δ_]:=(1/(π*Δ))*(1/(1+(((((x-x0)/Δ)))^2))); (* A Lorentzian function to mimic spectral lines of width Δ nm and center wavelength x0 *)
Plot[{sharpfunction[x,λ[1,2],0.1],sharpfunction[x,λ[1,3],0.1],sharpfunction[x,λ[1,4],0.1],sharpfunction[x,λ[1,5],0.1],sharpfunction[x,λ[1,6],0.1],sharpfunction[x,λ[1,7],0.1],sharpfunction[x,λ[1,8],0.1],sharpfunction[x,λ[2,3],0.1],sharpfunction[x,λ[2,4],0.1],sharpfunction[x,λ[2,5],0.1],sharpfunction[x,λ[2,6],0.1],sharpfunction[x,λ[2,7],0.1],sharpfunction[x,λ[2,8],0.1],sharpfunction[x,λ[2,9],0.1],sharpfunction[x,λ[3,4],0.1],sharpfunction[x,λ[3,5],0.1],sharpfunction[x,λ[3,6],0.1],sharpfunction[x,λ[3,7],0.1],sharpfunction[x,λ[3,8],0.1],sharpfunction[x,λ[3,9],0.1],sharpfunction[x,λ[3,10],0.1]},{x,0,2000},PlotRange->{0,0.2},Frame->True,FrameLabel->{"λ (nm)","Hydrogen Atom Specta"},PlotStyle->Thick,BaseStyle->{FontSize->15}]

Free Electron Bandstructure and Compton Wavelength

(* Free Electron Energy Bandstructure *)
h:=6.63*(10^(-34)); (* J.s, Planck's constant *)
hbar:=((6.63*(10^(-34)))/(2π)); (* J.s, Reduced Planck's constant *)
q:=1.6*(10^(-19 ));  (* Coulomb, electron charge *)
m0:=9.1*(10^(-31)); (* kg, free electron rest mass *)
freeenergy[k_]:=(((hbar^2)*(((k*(10^9)))^2))/(2*m0*q));  (* particle in a box energies, in eV units. Note k in (nm^(-1)) *)
Plot[{freeenergy[k]},{k,-10,10},Frame->True,FrameLabel->{"k ((nm^(-1)))","E(k) (eV)"},PlotStyle->Thick,BaseStyle->{FontSize->15},AspectRatio->GoldenRatio]  (* Note the energy and wavevector scales *)

(* Compton wavelength, and the de Broglie relation for free electrons *)
c=3*(10^8);  (* Speed of light, m/s *)
λc=((6.63*(10^(-34)))/(9.1*(10^(-31))*c));  (* Compton wavelength of electron *)
v[λ_]:=(c/(Sqrt[1+((((λ/(λc))))^2)]));  (* relativistic velocity *)
LogLogPlot[{v[0],v[λ]},{λ,(10^(-15)),(10^(-6))},Frame->True,FrameLabel->{"λ (m)","v (m/s)"},PlotStyle->{Dashed,Thick},BaseStyle->{FontSize->15},AspectRatio->GoldenRatio,GridLines->Automatic]

Schrodinger Equation for a linear potential, V(x)=x

soln=DSolve[{y''[x]==x*y[x],y[-3]==-1,y[+10]==0},y,x]
Plot[Evaluate[y[x]/.soln],{x,-20,5},Frame->True,FrameLabel->{"x","y(x)"},PlotStyle->Thick,BaseStyle->{FontSize->15}]

(* Linear Potential, Constant Electric Field *)
solnConstField:=NDSolve[{y''[x]==x*y[x],y[-3]==-1,y[+5]==0},y,{x,-40,10}];
Plot[Evaluate[y[x]/.solnConstField],{x,-20,5},Frame->True,FrameLabel->{"x","y(x)"},PlotStyle->Thick,BaseStyle->{FontSize->15}]
Plot[Evaluate[((y[x])^2)/.solnConstField],{x,-20,5},Frame->True,FrameLabel->{"x","|y(x)((|)^2)"},PlotStyle->Thick,BaseStyle->{FontSize->15}]

Particle in a Box with Infinite Barriers

ψ[n_,L_,x_]:=(Sqrt[(2/L)])*Sin[((n*π)/L)*x]*(UnitStep[x-L]-UnitStep[x]);  (* Eigenfunctions *)
Manipulate[Plot[ψ[n,2,x],{x,-1,3},Frame->True,PlotRange->{-1.1,1.1},FrameLabel->{"x","ψ(x)"},PlotStyle->Thick,BaseStyle->{FontSize->15}],{n,1,10,1}]

h:=6.63*(10^(-34)); (* J.s, Planck's constant *)
q:=1.6*(10^(-19 ));  (* Coulomb, electron charge *)
m0:=9.1*(10^(-31)); (* kg, free electron rest mass *)
boxenergy[n_,L_]:=(n^2)*((h^2)/(8*m0*(((L*(10^(-9))))^2)*q));  (* particle in a box energies, in eV units. Note the box length is in nm *)
Plot[{boxenergy[1,L],boxenergy[2,L],boxenergy[3,L]},{L,0,10},Frame->True,FrameLabel->{"L (nm)","(Subscript[E, 1])(eV)"},PlotStyle->Thick,BaseStyle->{FontSize->15},AspectRatio->GoldenRatio]  (* 3 lowest eigenvalues as a function of the well width *)

ListPlot[Table[boxenergy[n,2],{n,1,10,1}],Frame->True,FrameLabel->{"n","(Subscript[E, n])(eV)"},PlotStyle->Thick,BaseStyle->{FontSize->15},AspectRatio->GoldenRatio]  (* Eigenvalues are discrete,and grow as (n^2) *)

Particle in a Box with Finite Barriers

(* Finite well, number of bound states *)
hbar:=((6.63*(10^(-34)))/(2*π)); (* reduced Planck's constant, J.s *)
m0:=9.1*(10^(-31)); (* Free electron mass, kg *)
q:=1.6*(10^(-19)); (* electron charge, C *)
Eqw[meff_,Lw_]:=((((2*π*hbar))^2)/(q*8*m0*meff*(((Lw*10*(10^(-9))))^2)));
θ0[U0_,meff_,Lw_]:=(Sqrt[(((q*U0)*m0*meff*(((Lw*(10^(-9))))^2))/(2*(hbar^2)))]); (* U0 in eV, Lw in nm *)
lhs[U0_,meff_,Lw_,θ_]:=(Sqrt[(((((θ0[U0,meff,Lw])/θ)))^2)-1]);
rhs1[θ_]:=Tan[θ];
rhs2[θ_]:=-Cot[θ];
Manipulate[Plot[{lhs[U0,0.067,Lw,θ],rhs1[θ],rhs2[θ]},{θ,0,2π},Frame->True,PlotRange->{-6,6},FrameLabel->{"((Subscript[kL, w])/2)","LHS and RHS"},AspectRatio->GoldenRatio],{U0,0.1,1},{Lw,10,1}]

eigenθ1[U0_,meff_,Lw_]:=Quiet[Re[FindRoot[lhs[U0,meff,Lw,θ]==rhs1[θ],{θ,{0.01,θ0[U0,meff,Lw]}}][[1]][[2]]]];
eigenθ2[U0_,meff_,Lw_]:=Quiet[Re[FindRoot[lhs[U0,meff,Lw,θ]==rhs2[θ],{θ,{0.01,θ0[U0,meff,Lw]}}][[1]][[2]]]];
eigenEn1[U0_,meff_,Lw_]:=Sort[((2*(hbar^2)*((eigenθ1[U0,meff,Lw])^2))/(q*m0*meff*(((Lw*(10^(-9))))^2)))];
eigenEn2[U0_,meff_,Lw_]:=Sort[((2*(hbar^2)*((eigenθ2[U0,meff,Lw])^2))/(q*m0*meff*(((Lw*(10^(-9))))^2)))];
eigenEn1[0.3,0.067,5]
Quiet[Plot[{eigenEn1[0.3,0.067,Lw][[1]],eigenEn1[0.3,0.067,Lw][[2]],eigenEn1[0.3,0.067,Lw][[1]],eigenEn2[0.3,0.067,Lw][[1]],eigenEn2[0.3,0.067,Lw][[2]],eigenEn2[0.3,0.067,Lw][[3]]},{Lw,1,20},Frame->True,PlotRange->{0,0.3},FrameLabel->{"(Subscript[L, w])(nm)","Eigenvalue (eV)"}]]

Quiet[Plot[{eigenEn1[U0,0.067,10][[1]],eigenEn1[U0,0.067,10][[2]],eigenEn1[U0,0.067,10][[1]],eigenEn2[U0,0.067,10][[1]],eigenEn2[U0,0.067,10][[2]],eigenEn2[U0,0.067,10][[3]]},{U0,0.1,1.0},Frame->True,PlotRange->{0,1.0},FrameLabel->{"(Subscript[U, 0])(eV)","Eigenvalue (eV)"}]]

4. Damned Lies, and Statistics

Fermi-Dirac and Bose-Einstein Occupation Functions

kb:=1.38*(10^(-23));  (* Boltzmann constant, in J/K *)
q:=1.6*(10^(-19)); (* Electron charge, in Coulomb *) 
mb[en_,μ_,T_]:=Exp[-((q*(en-μ))/(kb*T))] ;(* The Maxwell-Boltzmann Function *)
fd[en_,μ_,T_]:=(1/(Exp[((q*(en-μ))/(kb*T))]+1)) ;(* The Fermi-Dirac Function *)
be[en_,μ_,T_]:=(1/(Exp[((q*(en-μ))/(kb*T))]-1)) ;(* The Bose-Einstein Function *)
Manipulate[Quiet[Plot[{mb[en,μ,T],fd[en,μ,T],fd[en,μ,1],be[en,μ,T]},{en,0,10},PlotRange->{0,2},Frame->True,AspectRatio->1,FrameLabel->{"Energy (eV)","Occupation Probability"},BaseStyle->{FontSize->15},PlotStyle->{{Red},{Blue},{Dashed,Gray},{Green}}]],{T,3000,4},{μ,5,1}]
(* Red: Maxwell-Boltzmann, Green: Bose-Einstein, Blue: Fermi-Dirac *)

Fermi-Dirac Integrals

(* Some useful constants *)
Print["Γ[-2]=",Gamma[-2],"  ","ζ[-2]=",Zeta[-2]]
Print["Γ[-3/2]=",Gamma[-3/2],"  ","ζ[-3/2]=",Zeta[-3/2]+0.0]
Print["Γ[-1]=",Gamma[-1],"  ","ζ[-1]=",Zeta[-1]]
Print["Γ[-1/2]=",Gamma[-1/2],"  ","ζ[-1/2]=",Zeta[-1/2]+0.0]
Print["Γ[0]=",Gamma[0],"  ","ζ[0]=",Zeta[0]]
Print["Γ[+1/2]=",Gamma[+1/2],"  ","ζ[+1/2]=",Zeta[+1/2]+0.0]
Print["Γ[+1]=",Gamma[+1],"  ","ζ[+1]=",Zeta[+1]]
Print["Γ[+3/2]=",Gamma[+3/2],"  ","ζ[+3/2]=",Zeta[+3/2]+0.0]
Print["Γ[+2]=",Gamma[+2],"  ","ζ[+2]=",Zeta[+2]]
Print["Γ[+5/2]=",Gamma[+5/2],"  ","ζ[+5/2]=",Zeta[+5/2]+0.0]
Print["Γ[+3]=",Gamma[+3],"  ","ζ[+3]=",Zeta[+3]+0.0]
Print["Γ[+7/2]=",Gamma[+7/2],"  ","ζ[+7/2]=",Zeta[+7/2]+0.0]
Print["Γ[+4]=",Gamma[+4],"  ","ζ[+4]=",Zeta[+4]]
Print["Γ[+9/2]=",Gamma[+9/2],"  ","ζ[+9/2]=",Zeta[+9/2]+0.0]
Print["Γ[+5]=",Gamma[+5],"  ","ζ[+5]=",Zeta[+5]+0.0]

(∫_0^∞)((x^j)/(Exp[x-μ]+1))\[DifferentialD]x  (* Integrating the Fermi-Dirac Function *)
(∫_0^∞)((x^j)/(Exp[x-μ]-1))\[DifferentialD]x   (* Integrating the Bose-Einstein Function *)

((-1)/(Gamma[1+j]))(∫_0^∞)((x^j)/(Exp[x-μ]+1))\[DifferentialD]x  (* Integrating the Fermi-Dirac Function *)

(∫_0^∞)((x*(Sqrt[x]))/(Exp[x-η]+1))\[DifferentialD]x
Series[(∫_0^∞)((x*(Sqrt[x]))/(Exp[x-η]+1))\[DifferentialD]x,{η,0,4}]

Series[-(1/2) (Sqrt[π]) PolyLog[(3/2),-(E^η)],{η,∞,0}]
Series[-(3/4) (Sqrt[π]) PolyLog[(5/2),-(E^η)],{η,∞,0}]

(* Defining the Fermi-Dirac Integral *)
FD[j_,η_]:=(1/(Gamma[j+1]))(∫_0^∞)((x^j)/(1+Exp[x-η]))\[DifferentialD]x;
Print[" FD[-2,η]= ",D[FD[0,η],{η,2}]]  (* Fermi-Dirac integral for j=-1 *)
Print[" FD[-1,η]= ",D[FD[0,η],{η,1}]]  (* Fermi-Dirac integral for j=-1 *)
Print[" FD[-(1/2),η]= ",FD[-(1/2),η]]  (* Fermi-Dirac integral for j=-(1/2) *)
Print[" FD[0,η]= ",FD[0,η]]  (* Fermi-Dirac integral for j=0 *)
Print[" FD[(1/2),η]= ",FD[(1/2),η]]  (* Fermi-Dirac integral for j=(1/2) *)
Print[" FD[1,η]= ",FD[1,η]]  (* Fermi-Dirac integral for j=1 *)
Print[" FD[(3/2),η]= ",FD[(3/2),η]]  (* Fermi-Dirac integral for j=(3/2) *)
Print[" FD[2,η]= ",FD[2,η]]  (* Fermi-Dirac integral for j=2 *)

-PolyLog[(1/2),-(E^0)]
-PolyLog[(1/2),-(E^0)]+0.0
((Sqrt[π])/2)*(-PolyLog[(1/2),-(E^0)]+0.0)
Log[1+(E^0)]
-PolyLog[(3/2),-(E^0)]

(* Plot of The Fermi-Dirac Integrals above *)
LogPlot[{-PolyLog[(3/2),-(E^η)],Log[1+(E^η)],-PolyLog[(1/2),-(E^η)],-PolyLog[2,-(E^η)]},{η,-8,8},Frame->True, AspectRatio->GoldenRatio,FrameLabel->{"η","Fermi-Dirac Integrals"},BaseStyle->{FontSize->15},PlotStyle->Thick,GridLines->Automatic]

(* Plot of The Fermi-Dirac Integrals above *)
LogPlot[{(E^η),-((E^(2 η))/(((1+(E^η)))^2))+((E^η)/(1+(E^η))),((E^η)/(1+(E^η))),-PolyLog[(1/2),-(E^η)],Log[1+(E^η)],-PolyLog[(3/2),-(E^η)],-PolyLog[2,-(E^η)],-PolyLog[(5/2),-(E^η)],-PolyLog[3,-(E^η)]},{η,-8,8},Frame->True, AspectRatio->GoldenRatio,FrameLabel->{"η","Fermi-Dirac Integrals"},BaseStyle->{FontSize->15},PlotRange->{(10^(-3)),(10^3)},GridLines->{{-8,-7,-6,-5,-4,-3,-2,-1,1,2,3,4,5,6,7,8},{(10^(-3)),2*(10^(-3)),3*(10^(-3)),4*(10^(-3)),5*(10^(-3)),6*(10^(-3)),7*(10^(-3)),8*(10^(-3)),9*(10^(-3)),(10^(-2)),2*(10^(-2)),3*(10^(-2)),4*(10^(-2)),5*(10^(-2)),6*(10^(-2)),7*(10^(-2)),8*(10^(-2)),9*(10^(-2)),(10^(-1)),2*(10^(-1)),3*(10^(-1)),4*(10^(-1)),5*(10^(-1)),6*(10^(-1)),7*(10^(-1)),8*(10^(-1)),9*(10^(-1)),(10^0),1,2,3,4,5,6,7,8,9,10,2*(10^1),3*(10^1),4*(10^1),5*(10^1),6*(10^1),7*(10^1),8*(10^1),9*(10^1),(10^2),2*(10^2),3*(10^2),4*(10^2),5*(10^2),6*(10^2),7*(10^2),8*(10^2),9*(10^2),(10^3)}}]

(* Plot of The Fermi-Dirac Integrals above *)
LogPlot[{(E^η),-((E^(2 η))/(((1+(E^η)))^2))+((E^η)/(1+(E^η))),((E^η)/(1+(E^η))),-PolyLog[(1/2),-(E^η)],Log[1+(E^η)],-PolyLog[(3/2),-(E^η)],-PolyLog[2,-(E^η)],-PolyLog[(5/2),-(E^η)],-PolyLog[3,-(E^η)]},{η,-5,15},Frame->True, AspectRatio->1,FrameLabel->{"η","Fermi-Dirac Integrals"},BaseStyle->{FontSize->15},PlotStyle->{{Dashed,Black},{Black},{Black},{Black},{Black},{Black},{Black},{Black},{Black}},PlotRange->{(10^(-2)),(10^3)},GridLines->{{-5,-4,-3,-2,-1,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15},{(10^(-2)),2*(10^(-2)),3*(10^(-2)),4*(10^(-2)),5*(10^(-2)),6*(10^(-2)),7*(10^(-2)),8*(10^(-2)),9*(10^(-2)),(10^(-1)),2*(10^(-1)),3*(10^(-1)),4*(10^(-1)),5*(10^(-1)),6*(10^(-1)),7*(10^(-1)),8*(10^(-1)),9*(10^(-1)),(10^0),1,2,3,4,5,6,7,8,9,10,2*(10^1),3*(10^1),4*(10^1),5*(10^1),6*(10^1),7*(10^1),8*(10^1),9*(10^1),(10^2),2*(10^2),3*(10^2),4*(10^2),5*(10^2),6*(10^2),7*(10^2),8*(10^2),9*(10^2),(10^3)}}]

(* Plot of The Fermi-Dirac Integrals above *)
LogPlot[{(E^η),-((E^(2 η))/(((1+(E^η)))^2))+((E^η)/(1+(E^η))),((E^η)/(1+(E^η))),-PolyLog[(1/2),-(E^η)],Log[1+(E^η)],-PolyLog[(3/2),-(E^η)],-PolyLog[2,-(E^η)],-PolyLog[(5/2),-(E^η)],-PolyLog[3,-(E^η)]},{η,-5,15},Frame->True, AspectRatio->1,FrameLabel->{"η","Fermi-Dirac Integrals"},BaseStyle->{FontSize->15},PlotStyle->{{Dashed,Thick,Black},{Thick,Black},{Thick,Black},{Thick,Black},{Thick,Black},{Thick,Black},{Thick,Black},{Thick,Black},{Thick,Black}},PlotRange->{(10^(-2)),(10^3)},GridLines->{{-5,-4,-3,-2,-1,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15},{(10^(-2)),2*(10^(-2)),3*(10^(-2)),4*(10^(-2)),5*(10^(-2)),6*(10^(-2)),7*(10^(-2)),8*(10^(-2)),9*(10^(-2)),(10^(-1)),2*(10^(-1)),3*(10^(-1)),4*(10^(-1)),5*(10^(-1)),6*(10^(-1)),7*(10^(-1)),8*(10^(-1)),9*(10^(-1)),(10^0),1,2,3,4,5,6,7,8,9,10,2*(10^1),3*(10^1),4*(10^1),5*(10^1),6*(10^1),7*(10^1),8*(10^1),9*(10^1),(10^2),2*(10^2),3*(10^2),4*(10^2),5*(10^2),6*(10^2),7*(10^2),8*(10^2),9*(10^2),(10^3)}}]

(* Plot of The Fermi-Dirac Integrals above *)
LogPlot[{(E^η),-PolyLog[(3/2),-(E^η)],(4/(3*(Sqrt[π])))*(η^(3/2))},{η,-8,8},Frame->True, AspectRatio->GoldenRatio,FrameLabel->{"η","Fermi-Dirac Integrals"},BaseStyle->{FontSize->15},PlotStyle->{{Dashed,Thick,Black},{Thick,Black},{Dashed,Thick,Black}},PlotRange->{(10^(-4)),(10^2)},GridLines->Automatic]

Some Statistical Integrals and Related Results

(((∫_(-∞)^∞)(1/2)*m*(vx^2)*Exp[-(((1/2)*m*(vx^2))/(k*T))]\[DifferentialD]vx)/((∫_(-∞)^∞)Exp[-(((1/2)*m*(vx^2))/(k*T))]\[DifferentialD]vx))  (* Classical 1D Equipartition of Energy *)

(((∫_(-∞)^∞)(∫_(-∞)^∞)(∫_(-∞)^∞)(1/2)*m*((vx^2)+(vy^2)+(vz^2))*Exp[-(((1/2)*m*((vx^2)+(vy^2)+(vz^2)))/(k*T))]\[DifferentialD]vx\[DifferentialD]vy\[DifferentialD]vz)/((∫_(-∞)^∞)(∫_(-∞)^∞)(∫_(-∞)^∞)Exp[-(((1/2)*m*((vx^2)+(vy^2)+(vz^2)))/(k*T))]\[DifferentialD]vx\[DifferentialD]vy\[DifferentialD]vz))  (* Classical 3D Equipartition of Energy *)

(((∫_(-∞)^∞)(1/2)*m*(vx^2)*(1/(1+Exp[(((1/2)*m*(vx^2)-μ)/(k*T))]))\[DifferentialD]vx)/((∫_(-∞)^∞)(1/(1+Exp[(((1/2)*m*(vx^2)-μ)/(k*T))]))\[DifferentialD]vx))  (* Quantum 1D Equipartition of Energy for Fermions *)
(((∫_(-∞)^∞)(∫_(-∞)^∞)(1/2)*m*((vx^2)+(vy^2))*(1/(Exp[(((1/2)*m*((vx^2)+(vy^2))-μ)/(k*T))]+1))\[DifferentialD]vx\[DifferentialD]vy)/((∫_(-∞)^∞)(∫_(-∞)^∞)(1/(Exp[(((1/2)*m*((vx^2)+(vy^2))-μ)/(k*T))]+1))\[DifferentialD]vx\[DifferentialD]vy))  (* Quantum 2D Equipartition of Energy for Fermions *)
(((∫_(-∞)^∞)(∫_(-∞)^∞)(∫_(-∞)^∞)(1/2)*m*((vx^2)+(vy^2)+(vz^2))*(1/(Exp[(((1/2)*m*((vx^2)+(vy^2)+(vz^2))-μ)/(k*T))]+1))\[DifferentialD]vx\[DifferentialD]vy\[DifferentialD]vz)/((∫_(-∞)^∞)(∫_(-∞)^∞)(∫_(-∞)^∞)(1/(Exp[(((1/2)*m*((vx^2)+(vy^2)+(vz^2))-μ)/(k*T))]+1))\[DifferentialD]vx\[DifferentialD]vy\[DifferentialD]vz))  (* Quantum 3D Equipartition of Energy for Fermions *)

QuantFermionEquipartEnergy1D[q_,kb_,T_,μ_]:=((kb*T)/(2*q))*((PolyLog[(3/2),-(E^((q*μ)/(kb*T)))])/(PolyLog[(1/2),-(E^((q*μ)/(kb*T)))]));
QuantFermionEquipartEnergy2D[q_,kb_,T_,μ_]:=((kb*T)/q)*((-PolyLog[(4/2),-(E^((q*μ)/(kb*T)))])/(Log[1+(E^((q*μ)/(kb*T)))]));
QuantFermionEquipartEnergy3D[q_,kb_,T_,μ_]:=((3 *kb*T)/(2*q))*((PolyLog[(5/2),-(E^((q*μ)/(kb*T)))])/(PolyLog[(3/2),-(E^((q*μ)/(kb*T)))]));
Quiet[LogPlot[{QuantFermionEquipartEnergy1D[1.6*(10^(-19)),1.38*(10^(-23)),1000,μ],QuantFermionEquipartEnergy1D[1.6*(10^(-19)),1.38*(10^(-23)),300,μ],QuantFermionEquipartEnergy1D[1.6*(10^(-19)),1.38*(10^(-23)),4,μ]},{μ,-0.5,1},PlotRange->{(10^(-5)),(10^0)},Frame->True,AspectRatio->GoldenRatio, FrameLabel->{"μ (eV)","1D Fermionic Equipartition Energy (eV)"},BaseStyle->{FontSize->15},PlotStyle->{{Red,Thick},{Green,Thick},{Blue,Thick}},GridLines->Automatic]]
Quiet[LogPlot[{QuantFermionEquipartEnergy2D[1.6*(10^(-19)),1.38*(10^(-23)),1000,μ],QuantFermionEquipartEnergy2D[1.6*(10^(-19)),1.38*(10^(-23)),300,μ],Chop[QuantFermionEquipartEnergy2D[1.6*(10^(-19)),1.38*(10^(-23)),4,μ]]},{μ,-0.5,1},PlotRange->{(10^(-5)),(10^0)},Frame->True,AspectRatio->GoldenRatio, FrameLabel->{"μ (eV)","2D Fermionic Equipartition Energy (eV)"},BaseStyle->{FontSize->15},PlotStyle->{{Red,Thick},{Green,Thick},{Blue,Thick}},GridLines->Automatic]]
Quiet[LogPlot[{QuantFermionEquipartEnergy3D[1.6*(10^(-19)),1.38*(10^(-23)),1000,μ],QuantFermionEquipartEnergy3D[1.6*(10^(-19)),1.38*(10^(-23)),300,μ],QuantFermionEquipartEnergy3D[1.6*(10^(-19)),1.38*(10^(-23)),4,μ]},{μ,-0.5,1},PlotRange->{(10^(-5)),(10^0)},Frame->True,AspectRatio->GoldenRatio, FrameLabel->{"μ (eV)","3D Fermionic Equipartition Energy (eV)"},BaseStyle->{FontSize->15},PlotStyle->{{Red,Thick},{Green,Thick},{Blue,Thick}},GridLines->Automatic]]

(((∫_(-∞)^∞)(1/2)*m*(vx^2)*(1/(Exp[(((1/2)*m*(vx^2)-μ)/(k*T))]-1))\[DifferentialD]vx)/((∫_(-∞)^∞)(1/(Exp[(((1/2)*m*(vx^2)-μ)/(k*T))]-1))\[DifferentialD]vx))  (* Quantum 1D Equipartition of Energy for Bosons *)
(((∫_(-∞)^∞)(∫_(-∞)^∞)(1/2)*m*((vx^2)+(vy^2))*(1/(Exp[(((1/2)*m*((vx^2)+(vy^2))-μ)/(k*T))]-1))\[DifferentialD]vx\[DifferentialD]vy)/((∫_(-∞)^∞)(∫_(-∞)^∞)(1/(Exp[(((1/2)*m*((vx^2)+(vy^2))-μ)/(k*T))]-1))\[DifferentialD]vx\[DifferentialD]vy))  (* Quantum 2D Equipartition of Energy for Bosons *)
(((∫_(-∞)^∞)(∫_(-∞)^∞)(∫_(-∞)^∞)(1/2)*m*((vx^2)+(vy^2)+(vz^2))*(1/(Exp[(((1/2)*m*((vx^2)+(vy^2)+(vz^2))-μ)/(k*T))]-1))\[DifferentialD]vx\[DifferentialD]vy\[DifferentialD]vz)/((∫_(-∞)^∞)(∫_(-∞)^∞)(∫_(-∞)^∞)(1/(Exp[(((1/2)*m*((vx^2)+(vy^2)+(vz^2))-μ)/(k*T))]-1))\[DifferentialD]vx\[DifferentialD]vy\[DifferentialD]vz))  (* Quantum 3D Equipartition of Energy for Bosons *)

QuantBosonEquipartEnergy3D[q_,kb_,T_,μ_]:=((3 *kb*T)/(2*q))*((PolyLog[(5/2),(E^((q*μ)/(kb*T)))])/(PolyLog[(3/2),(E^((q*μ)/(kb*T)))]));
LogPlot[{QuantBosonEquipartEnergy3D[1.6*(10^(-19)),1.38*(10^(-23)),1000,μ],QuantBosonEquipartEnergy3D[1.6*(10^(-19)),1.38*(10^(-23)),300,μ],QuantBosonEquipartEnergy3D[1.6*(10^(-19)),1.38*(10^(-23)),4,μ]},{μ,-0.2,1},PlotRange->{(10^(-5)),(10^0)},Frame->True,AspectRatio->GoldenRatio, FrameLabel->{"μ (eV)","Bosonic Equipartition Energy (eV)"},BaseStyle->{FontSize->15},PlotStyle->{{Red,Thick},{Green,Thick},{Blue,Thick}},GridLines->Automatic]

Print["(((∫_0^∞)(x^2)*Exp[-β*x] \[DifferentialD]x)/((∫_0^∞)x*Exp[-β*x] \[DifferentialD]x))= ",(((∫_0^∞)(x^2)*Exp[-β*x]\[DifferentialD]x)/((∫_0^∞)x*Exp[-β*x]\[DifferentialD]x))]
Print["(∑_(n = 0)^∞)(n*h*ω)*Exp[-((n*h*ω)/(k*T))]= ",(∑_(n=0)^∞)(n*h*ω)*Exp[-((n*h*ω)/(k*T))]]
Print["(((∑_(n = 0)^∞)(n*h*ω)*Exp[-((n*h*ω)/(k*T))])/((∑_(n = 0)^∞)Exp[-((n*h*ω)/(k*T))]))= ",(((∑_(n=0)^∞)(n*h*ω)*Exp[-((n*h*ω)/(k*T))])/((∑_(n=0)^∞)Exp[-((n*h*ω)/(k*T))]))]
Print["(∑_(n = 0)^∞)((ω^2)/((π^2)*(c^3)))*(n*h*ω)*Exp[-((n*h*ω)/(k*T))]= ",(∑_(n=0)^∞)((ω^2)/((π^2)*(c^3)))*(n*h*ω)*Exp[-((n*h*ω)/(k*T))]]
Print["(((∑_(n = 0)^∞)((ω^2)/(h*(π^2)*(c^3)))*(n*h*ω)*Exp[-((n*h*ω)/(k*T))])/((∑_(n = 0)^∞)((ω^2)/(h*(π^2)*(c^3)))*Exp[-((n*h*ω)/(k*T))]))= ",(((∑_(n=0)^∞)((ω^2)/(h*(π^2)*(c^3)))*(n*h*ω)*Exp[-((n*h*ω)/(k*T))])/((∑_(n=0)^∞)((ω^2)/(h*(π^2)*(c^3)))*Exp[-((n*h*ω)/(k*T))]))](* Planck's Quantum Hypothesis *)
Print["(∫_0^∞)(((x^2)*(E^x))/((((E^x) - 1))^2))\[DifferentialD]x= ",(∫_0^∞)(((x^2)*(E^x))/((((E^x)-1))^2))\[DifferentialD]x]
Print["(∫_0^∞)(((x^3)*(E^x))/((((E^x) - 1))^2))\[DifferentialD]x= ",(∫_0^∞)(((x^3)*(E^x))/((((E^x)-1))^2))\[DifferentialD]x]
Print["(∫_0^∞)(((x^4)*(E^x))/((((E^x) - 1))^2))\[DifferentialD]x= ",(∫_0^∞)(((x^4)*(E^x))/((((E^x)-1))^2))\[DifferentialD]x]
Print["(∫_0^∞)((x^3)/((E^x) - 1))\[DifferentialD]x= ",(∫_0^∞)((x^3)/((E^x)-1))\[DifferentialD]x]
Print["(∫_0^∞)(x^4)*Csch[\!\(\*FractionBox[\(x\), \(\(2\)\(\\ \)\)]\)]^4*Sinh[x]\[DifferentialD]x= ",(∫_0^∞)(x^4)*((Csch[(x/(2 ))])^4)*Sinh[x]\[DifferentialD]x]
Print["D[\!\(\*FractionBox[\(h\\ ω\), \((\(-1\) + \*SuperscriptBox[\(E\), FractionBox[\(h\\ ω\), \(k\\ T\)]])\)]\),T]= ",D[((h ω)/((-1+(E^((h ω)/(k T)))))),T]]
Print["D[((T/(ΘD)))^4*(∫_0^((ΘD)/T))((x^3)/((E^x) - 1))\[DifferentialD]x,T]= ",FullSimplify[D[((((T/(ΘD))))^4)*(∫_0^((ΘD)/T))((x^3)/((E^x)-1))\[DifferentialD]x,T]]]

(* Specific heat capacity: Dulong-Petit vs Einstein vs Debye *)
Debye[ΘD_,T_]:=3*((1/(-T+(E^(-((ΘD)/T))) T))-((4 (π^4) (T^3))/(15 ((ΘD)^4)))+((4 Log[1-(E^((ΘD)/T))])/(ΘD))+(1/((ΘD)^4))12 T (((ΘD)^2) PolyLog[2,(E^((ΘD)/T))]+2 T (-ΘD PolyLog[3,(E^((ΘD)/T))]+T PolyLog[4,(E^((ΘD)/T))])));
Einstein[ΘD_,T_]:=(((((((ΘD)/T)))^2)*(E^((ΘD)/T)) )/((((E^((ΘD)/T))-1))^2));
Plot[{1,Einstein[1,T],Chop[N[Debye[1,T]]]},{T,0,1},Frame->True, FrameLabel->{"T [(Subscript[ℏω, D])/(Subscript[k, B])]","(Subscript[c, v]) [3(Subscript[nk, B])]"},BaseStyle->{FontSize->15},PlotRange->{-0.1,1.2},PlotStyle->{{Black,Dashed},{Red,Thick},{Blue,Thick}},GridLines->Automatic]
Plot[{1,Einstein[1,T],Chop[N[Debye[1,T]]]},{T,0,3},Frame->True, FrameLabel->{"T [(Subscript[ℏω, D])/(Subscript[k, B])]","(Subscript[c, v]) [3(Subscript[nk, B])]"},BaseStyle->{FontSize->15},PlotRange->{-0.1,1.2},PlotStyle->{{Black,Dashed},{Red,Thick},{Blue,Thick}},GridLines->Automatic]

5. Electrons in the Quantum World

1D Carrier Statistics, Currents, Energies, and Specific Heat

(* Carrier Statistics And Quantum Currents in 1, 2, and 3 Dimensions *)
ClearAll
h:=6.63*(10^(-34)); (* Planck's constant in J.s *)
q:=1.6*(10^(-19));               (* Electron charge, Coulomb *)
hbar:=(6.63/(2π))*(10^(-34));   (* Reduced Planck's constant, J.s *)
kb:=1.38*(10^(-23));           (* Boltzmann constant, J/K *)
me:=9.1*(10^(-31));             (* Electron rest mass, Kg *)
FermiDiracIntegral[d_,η_]:=Switch[d,-1,((E^η)/(1+(E^η))),-(1/2),-PolyLog[(1/2),-(E^η)],0,Log[1+(E^η)],(1/2),-PolyLog[(3/2),-(E^η)],1,-PolyLog[2,-(E^η)],(3/2),-PolyLog[(5/2),-(E^η)],2,-PolyLog[3,-(E^η)],(5/2),-PolyLog[(7/2),-(E^η)],3,-PolyLog[4,-(E^η)]]; (* Define Fermi-Dirac Integrals as PolyLog functions for fast numerical evaluation *)
Nc[d_,gs_,gv_,meff_,T_]:=gs*gv*(((((2*π*meff*me*kb*T)/(h^2))))^(d/2)); (* Band-edge DOS in 1/(m^d) for d dimensions *)
n[d_,gs_,gv_,meff_,T_,Ef_]:=Nc[d,gs,gv,meff,T]*FermiDiracIntegral[((d-2)/2),(Ef/((kb*T)/q))];  (* Carrier density in: 1d in 1/(m^1), 2d in 1/(m^2), and 3d in 1/(m^3) *)
LogPlot[{(10^(-2))*n[1,2,1,1,77,Ef],(10^(-2))*n[1,2,1,1,200,Ef],(10^(-2))*n[1,2,1,1,300,Ef]},{Ef,-0.4,1},Frame->True, FrameLabel->{"(Subscript[E, F]) (eV)","(Subscript[n, 1  d])(1/cm)"},BaseStyle->{FontSize->15},PlotStyle->{{Thick,Blue},{Thick,Green},{Thick,Red}},AspectRatio->GoldenRatio,PlotRange->{(10^(-1)),(10^8)}]
(* Carrier density in: 1d in 1/(cm^1), LogScale *)
Plot[{(10^(-2))*n[1,2,1,1,77,Ef],(10^(-2))*n[1,2,1,1,200,Ef],(10^(-2))*n[1,2,1,1,300,Ef]},{Ef,-0.4,1},Frame->True, FrameLabel->{"(Subscript[E, F]) (eV)","(Subscript[n, 1  d])(1/cm)"},BaseStyle->{FontSize->15},PlotStyle->{{Thick,Blue},{Thick,Green},{Thick,Red}},AspectRatio->GoldenRatio,PlotRange->{-(10^6),4*(10^7)}]
(* Carrier density in: 1d in 1/(cm^1), LinearScale *)
EF[d_,nd_,gs_,gv_,meff_,T_]:=Chop[Quiet[FindRoot[FermiDiracIntegral[((d-2)/2),(Efs/((kb*T)/q))]+FermiDiracIntegral[((d-2)/2),(Efs/((kb*T)/q))]==((2*nd)/(Nc[d,gs,gv,meff,T]*(10^(-2*d)))),{Efs,0}]][[1,2]]]; (* nd is in 1/(cm^d) *)
EF[1,1*(10^7),2,1,1,300]
(*Plot[{EF[1,n1d*(10^7),2,1,1,77],EF[1,n1d*(10^7),2,1,1,200],EF[1,n1d*(10^7),2,1,1,300]},{n1d,0,10},Frame->True, FrameLabel->{"(Subscript[n, 1 d])(1/cm)","(Subscript[E, F]) (eV)"},BaseStyle->{FontSize->15},PlotStyle->{{Thick,Blue},{Thick,Green},{Thick,Red}},AspectRatio->GoldenRatio,PlotRange->All]
LogLinearPlot[{EF[1,n1d*(10^7),2,1,1,77],EF[1,n1d*(10^7),2,1,1,200],EF[1,n1d*(10^7),2,1,1,300]},{n1d,(10^(-5)),(10^1)},Frame->True, FrameLabel->{"(Subscript[n, 1 d])(1/cm)","(Subscript[E, F]) (eV)"},BaseStyle->{FontSize->15},PlotStyle->{{Thick,Blue},{Thick,Green},{Thick,Red}},AspectRatio->GoldenRatio,PlotRange->All]*)
ud[d_,T_,η_]:=(d/2)*((kb*T)/q)*((FermiDiracIntegral[(d/2),η])/(FermiDiracIntegral[((d-2)/2),η])); (* Quantum Equipartition of Energy in eV *)
LogPlot[{ud[1,4,((q*Ef)/(kb*4))],ud[1,300,((q*Ef)/(kb*300))],ud[1,1000,((q*Ef)/(kb*1000))]},{Ef,-0.5,1},Frame->True, FrameLabel->{"(Subscript[E, F])(eV)","(Subscript[u, 1  d])(eV)"},BaseStyle->{FontSize->15},PlotStyle->{{Thick,Blue},{Thick,Green},{Thick,Red}},AspectRatio->GoldenRatio,PlotRange->{(10^(-5)),1}]

(* Nc[d_,gs_,gv_,meff_,T_] *)
LogLogPlot[{(10^(-2))*Nc[1,2,1,1,T],(10^(-4))*Nc[2,2,1,1,T],(10^(-6))*Nc[3,2,1,1,T]},{T,1,1100},Frame->True, FrameLabel->{"T(K)","(Subscript[N, c]^d)(1/(cm^d))"},BaseStyle->{FontSize->15},PlotStyle->{{Thick,Blue},{Thick,Green},{Thick,Red}},AspectRatio->GoldenRatio,PlotRange->{(10^5),(10^21)},GridLines->Automatic]

(10^(-2))*Nc[1,2,1,1,300]
(10^(-4))*Nc[2,2,1,1,300]
(10^(-6))*Nc[3,2,1,1,300]

ud[d_,T_,η_]:=(d/2)*kb*T*((FermiDiracIntegral[(d/2),η])/(FermiDiracIntegral[((d-2)/2),η])); (* Quantum Equipartition of Energy *)
q:=1.6*(10^(-19));               (* Electron charge, Coulomb *)
kb:=1.38*(10^(-23));           (* Boltzmann constant, J/K *)
cv[d_,T_,Ef_]=((D[ud[d,T,((q*Ef)/(kb*T))],T])/kb);Plot[{(1/2),cv[1,T,0.1],cv[1,T,10],((2*1)/2),cv[2,T,0.1],cv[2,T,10],((3*1)/2),cv[3,T,0.1],cv[3,T,10]},{T,0,200000}]

cv[1,200000,0.1]
N[Log[2]]

(* Carrier Statistics And Quantum Currents in 1, 2, and 3 Dimensions *)
ClearAll
h:=6.63*(10^(-34)); (* Planck's constant in J.s *)
q:=1.6*(10^(-19));               (* Electron charge, Coulomb *)
hbar:=(6.63/(2π))*(10^(-34));   (* Reduced Planck's constant, J.s *)
kb:=1.38*(10^(-23));           (* Boltzmann constant, J/K *)
me:=9.1*(10^(-31));             (* Electron rest mass, Kg *)
n[d_,gs_,gv_,meff_,T_,Ef_]:=Nc[d,gs,gv,meff,T]*Switch[d,1,-PolyLog[(1/2),-(E^η)],2,Log[1+(E^η)],3,-PolyLog[(3/2),-(E^η)]];  (* Carrier density in 1d, 2d, and 3d in 1/(m^3) *)
Nc[d_,gs_,gv_,meff_,T_]:=gs*gv*(((((2*π*meff*me*kb*T)/(h^2))))^(d/2)); (* Band-edge DOS in 1/(m^d) for d dimensions *)
FermiDiracIntegral[d_,η_]:=Switch[d,-1,((E^η)/(1+(E^η))),-(1/2),-PolyLog[(1/2),-(E^η)],0,Log[1+(E^η)],(1/2),-PolyLog[(3/2),-(E^η)],1,-PolyLog[2,-(E^η)],(3/2),-PolyLog[(5/2),-(E^η)],2,-PolyLog[3,-(E^η)],(5/2),-PolyLog[(7/2),-(E^η)],3,-PolyLog[4,-(E^η)]]; (* Define Fermi-Dirac Integrals as PolyLog functions for fast numerical evaluation *)
EFs[d_,nd_,gs_,gv_,meff_,T_,v_]:=Chop[Quiet[FindRoot[FermiDiracIntegral[((d-2)/2),(Efs/((kb*T)/q))]+FermiDiracIntegral[((d-2)/2),((Efs-v)/((kb*T)/q))]==((2*nd)/(Nc[d,gs,gv,meff,T]*(10^(-2*d)))),{Efs,0}]][[1,2]]]; (* nd is in 1/(cm^d) *)
CurrentPreFactor[d_]:=Switch[d,1,2,2,(1/2),3,((3(Sqrt[π]))/4)]  (* Prefactor for the current density *)
jd[d_,nd_,gs_,gv_,meff_,T_,v_]:=Chop[((q^2)/h)*Nc[d-1,gs,gv,meff,T]*((kb*T)/q)*(FermiDiracIntegral[((d-1)/2),((EFs[d,nd,gs,gv,meff,T,v])/((kb*T)/q))]-FermiDiracIntegral[((d-1)/2),((EFs[d,nd,gs,gv,meff,T,v]-v)/((kb*T)/q))])];
(* Current Density in A/(m^(d-1)) *)

(* n[d_,gs_,gv_,meff_,T_,Ef_] *)
(* Nc[d_,gs_,gv_,meff_,T_] *)
Print["Nc[1,2,1,0.2,300]*(10^(-2))=",Nc[1,2,1,0.2,300]*(10^(-2)),"(cm^(-1))"] (* 1D Band Edge DOS in 1/cm *)
Print["Nc[2,2,1,0.2,300]*(10^(-2))=",Nc[2,2,1,0.2,300]*(10^(-4)),"(cm^(-2))"] (* 2D Band Edge DOS in 1/(cm^2) *)
Print["Nc[3,2,1,0.2,300]*(10^(-2))=",Nc[3,2,1,0.2,300]*(10^(-6)),"(cm^(-3))"] (* 3D Band Edge DOS in 1/(cm^3) *)
LogPlot[n[2,2,1,0.2,300,η]*(10^(-4)),{η,-10,5},Frame->True, FrameLabel->{"\!\(\*FractionBox[\(\*SubscriptBox[\(E\), \(F\)] - \*SubscriptBox[\(E\), \(c\)]\), \(kb\\ T\)]\)","n2d((cm^(-2)))"},BaseStyle->{FontSize->15},PlotStyle->Thick] (* 2D Carrier Density as a function of ((Subscript[E, F])-(Subscript[E, c]))/kb*T *)
Print["jd[1,(10^6),2,1,0.2,300,1]=",jd[1,(10^6),2,1,0.2,300,1],"A"] (* 1D Current in A *)
Print["jd[2,(10^13),2,1,0.2,300,1]=",jd[2,(10^13),2,1,0.2,300,1],"A/m"] (* 2D Current in A *)
Print["jd[3,(10^18),2,1,0.2,300,1]=",jd[3,(10^18),2,1,0.2,300,1],"A/(m^2)"] (* 3D Current in A *)

(* 1d Metal *)
(*ηs[d_,nd_,gs_,gv_,meff_,T_,v_]
EFs[d_,nd_,gs_,gv_,meff_,T_,v_]
jd[d_,nd_,gs_,gv_,meff_,T_,v_]*)
EFs[1,(10^7),2,1,1,300,1]
Plot[{EFs[1,0.1*(10^7),2,1,1,300,v],EFs[1,0.1*(10^7),2,1,1,300,v]-v,EFs[1,0.1*(10^7),2,1,1,300,0],EFs[1,5*(10^7),2,1,1,300,v],EFs[1,5*(10^7),2,1,1,300,v]-v,EFs[1,5*(10^7),2,1,1,300,0]},{v,0,15},Frame->True, FrameLabel->{"V (Volt)","(Subscript[E, Fs])-(Subscript[E, c]),(Subscript[E, Fd])-\!\(\*SubscriptBox[\(E\), \(\(c\)\(\\ \)\)]\)( in eV)"},BaseStyle->{FontSize->15},PlotRange->{-10,10},PlotStyle->{{Thick,Blue},{Thick,Blue},{Dashed,Blue},{Thick,Red},{Thick,Red},{Dashed,Red}}]
jd[1,(10^7),2,2,1,300,1]
Plot[{jd[1,0.1*(10^7),2,1,1,300,v],jd[1,(10^7),2,1,1,300,v],jd[1,2*(10^7),2,1,1,300,v],jd[1,3*(10^7),2,1,1,300,v],jd[1,4*(10^7),2,1,1,300,v],jd[1,5*(10^7),2,1,1,300,v]},{v,0,15},Frame->True, FrameLabel->{"V (Volt)","(Subscript[J, 1  d]) (Amp)"},BaseStyle->{FontSize->15},PlotStyle->{{Thick,RGBColor[0.085069, 0.157214, 0.922347]},{Thick,RGBColor[0, 0.699168, 0.901411]},{Thick,RGBColor[0, 0.92462, 0.228962]},{Thick,RGBColor[0.907683, 0.742321, 0.19472]},{Thick,RGBColor[0.956741, 0.282475, 0.188861]},{Thick,RGBColor[0.97229, 0, 0.0719463]}}]

(* Effective Classical Ensemble Velocity in d-dimensions *)
vd[d_,nd_,gs_,gv_,meff_,T_,v_]:=(Sqrt[((kb*T)/(2*π*me))])*(((FermiDiracIntegral[((d-1)/2),((EFs[d,nd,gs,gv,meff,T,v])/((kb*T)/q))]-FermiDiracIntegral[((d-1)/2),((EFs[d,nd,gs,gv,meff,T,v]-v)/((kb*T)/q))])/(FermiDiracIntegral[((d-2)/2),((EFs[d,nd,gs,gv,meff,T,0])/((kb*T)/q))])));
LogLinearPlot[{vd[1,n1*(10^7),2,1,1,300,5]*(10^2),vd[1,n1*(10^7),2,1,1,77,5]*(10^2),vd[1,n1*(10^7),2,1,1,4,5]*(10^2),vd[2,n1*(10^14),2,1,1,300,5]*(10^2),vd[2,n1*(10^14),2,1,1,77,5]*(10^2),vd[2,n1*(10^14),2,1,1,4,5]*(10^2),vd[3,n1*(10^23),2,1,1,300,5]*(10^2),vd[3,n1*(10^23),2,1,1,77,5]*(10^2),vd[3,n1*(10^23),2,1,1,4,5]*(10^2)},{n1,0.001,3},PlotRange->{0,5*(10^7)},Frame->True]

(* Effective Classical Ensemble Velocity in d-dimensions *)
vdEf[d_,x_,T_,v_]:=(Sqrt[((kb*T)/(2*π*me))])*((FermiDiracIntegral[((d-1)/2),(x/((kb*T)/q))]-FermiDiracIntegral[((d-1)/2),((x-v)/((kb*T)/q))])/FermiDiracIntegral[((d-2)/2),(x/((kb*T)/q))]); 
Quiet[Plot[{vdEf[1,x,300,5]*(10^2),vdEf[1,x,77,5]*(10^2),vdEf[1,x,4,5]*(10^2),vdEf[2,x,300,5]*(10^2),vdEf[2,x,77,5]*(10^2),vdEf[2,x,4,5]*(10^2),vdEf[3,x,300,5]*(10^2),vdEf[3,x,77,5]*(10^2),vdEf[3,x,4,5]*(10^2)},{x,-0.2,1.0},PlotRange->{0,2*(10^7)},Frame->True, FrameLabel->{"(Subscript[E, F]) (eV)","<v(Subscript[>, d]) (cm/s)"},BaseStyle->{FontSize->15},AspectRatio->1,PlotStyle->{{Thick,Red},{Thick,Green},{Thick,Blue},{Thick,Red},{Thick,Green},{Thick,Blue},{Thick,Red},{Thick,Green},{Thick,Blue}}]]

EFs[1,5*(10^7),2,1,1,300,0]
(2^(2/1))*EFs[1,5*(10^7),2,1,1,300,0]

(*EFs[d_,nd_,gs_,gv_,meff_,T_,v_]*)
Vsat[d_,nd_,gs_,gv_,meff_,T_]:=EFs[d,nd,gs,gv,meff,T,0]*(2^(2/d));  (* Saturation voltage for d dimensions *)
Plot[Vsat[1,n1d,2,1,1,300],{n1d,0,5*(10^7)}]
Plot[Vsat[2,n2d,2,1,1,300],{n2d,0,(10^15)}]
Plot[Vsat[3,n3d,2,1,1,300],{n3d,0,(10^23)}]

Gamma[(5/2)]

(* 1d Semiconductor *)
(*ηs[d_,nd_,gs_,gv_,meff_,T_,v_]
EFs[d_,nd_,gs_,gv_,meff_,T_,v_]
jd[d_,nd_,gs_,gv_,meff_,T_,v_]*)
EFs[1,(10^6),2,1,0.2,300,1]
Plot[{EFs[1,5*(10^6),2,1,0.2,300,v],EFs[1,5*(10^6),2,1,0.2,300,v]-v},{v,0,2},Frame->True, FrameLabel->{"V (Volt)","(Subscript[E, Fs])-(Subscript[E, c]),(Subscript[E, Fd])-\!\(\*SubscriptBox[\(E\), \(\(c\)\(\\ \)\)]\)( in eV)"},BaseStyle->{FontSize->15},PlotStyle->Thick]
jd[1,(10^6),2,2,0.2,300,1]
Plot[{jd[1,5*(10^6),2,1,0.2,300,v],jd[1,4*(10^6),2,1,0.2,300,v],jd[1,3*(10^6),2,1,0.2,300,v],jd[1,2*(10^6),2,1,0.2,300,v],jd[1,1*(10^6),2,1,0.2,300,v],jd[1,1*(10^5),2,1,0.2,300,v]},{v,0,2},Frame->True, FrameLabel->{"V (Volt)","(Subscript[J, 1  d]) (Amp)"},BaseStyle->{FontSize->15},PlotStyle->Thick]
Plot[{jd[1,5*(10^6),2,1,0.2,77,v],jd[1,4*(10^6),2,1,0.2,77,v],jd[1,3*(10^6),2,1,0.2,77,v],jd[1,2*(10^6),2,1,0.2,77,v],jd[1,1*(10^6),2,1,0.2,77,v],jd[1,1*(10^5),2,1,0.2,77,v]},{v,0,2},Frame->True, FrameLabel->{"V (Volt)","(Subscript[J, 1  d]) (Amp)"},BaseStyle->{FontSize->15},PlotStyle->Thick]

(* 2d Metal *)
(*ηs[d_,nd_,gs_,gv_,meff_,T_,v_]
jd[d_,nd_,gs_,gv_,meff_,T_,v_]*)
EFs[2,(10^14),2,1,1,300,1]
Plot[{EFs[2,(10^15),2,1,1,300,v],EFs[2,(10^15),2,1,1,300,v]-v},{v,0,10},Frame->True, FrameLabel->{"V (Volt)","(Subscript[E, Fs])-(Subscript[E, c]),(Subscript[E, Fd])-\!\(\*SubscriptBox[\(E\), \(\(c\)\(\\ \)\)]\)( in eV)"},BaseStyle->{FontSize->15},PlotStyle->Thick]
jd[2,(10^15),2,2,1,300,1]
Plot[{jd[2,1*(10^14),2,1,1,300,v],jd[2,2*(10^14),2,1,1,300,v],jd[2,4*(10^14),2,1,1,300,v],jd[2,8*(10^14),2,1,1,300,v],jd[2,(10^15),2,1,1,300,v]},{v,0,15},Frame->True, FrameLabel->{"V (Volt)","(Subscript[J, 2  d]) (Amp/m)"},BaseStyle->{FontSize->15},PlotStyle->Thick]

(* 2d Semiconductor *)
(*ηs[d_,nd_,gs_,gv_,meff_,T_,v_]
jd[d_,nd_,gs_,gv_,meff_,T_,v_]*)
EFs[2,5*(10^13),2,1,0.2,300,1]
Plot[{EFs[2,5*(10^13),2,1,0.2,300,v],EFs[2,5*(10^13),2,1,0.2,300,v]-v},{v,0,2},Frame->True, FrameLabel->{"V (Volt)","(Subscript[E, Fs])-(Subscript[E, c]),(Subscript[E, Fd])-\!\(\*SubscriptBox[\(E\), \(\(c\)\(\\ \)\)]\)( in eV)"},BaseStyle->{FontSize->15},PlotStyle->Thick]
jd[2,(10^15),2,2,1,300,1]
Plot[(10^(-3))*{jd[2,5*(10^13),2,1,0.2,300,v],jd[2,4*(10^13),2,1,0.2,300,v],jd[2,3*(10^13),2,1,0.2,300,v],jd[2,2*(10^13),2,1,0.2,300,v],jd[2,1*(10^13),2,1,0.2,300,v],jd[2,1*(10^12),2,1,0.2,300,v]},{v,0,2},Frame->True, FrameLabel->{"V (Volt)","(Subscript[J, 2  d]) (mA/μm)"},BaseStyle->{FontSize->15},PlotStyle->Thick]
Plot[(10^(-3))*{jd[2,5*(10^13),2,1,0.2,77,v],jd[2,4*(10^13),2,1,0.2,77,v],jd[2,3*(10^13),2,1,0.2,77,v],jd[2,2*(10^13),2,1,0.2,77,v],jd[2,1*(10^13),2,1,0.2,77,v],jd[2,1*(10^12),2,1,0.2,77,v]},{v,0,2},Frame->True, FrameLabel->{"V (Volt)","(Subscript[J, 2  d]) (mA/μm)"},BaseStyle->{FontSize->15},PlotStyle->Thick]

(* 3d Metal *)
(*ηs[d_,nd_,gs_,gv_,meff_,T_,v_]
jd[d_,nd_,gs_,gv_,meff_,T_,v_]*)
EFs[3,(10^23),2,1,1,300,1]
Plot[{EFs[3,(10^22),2,1,1,300,v],EFs[3,(10^22),2,1,1,300,v]-v},{v,0,15},Frame->True, FrameLabel->{"V (Volt)","(Subscript[E, Fs])-(Subscript[E, c]),(Subscript[E, Fd])-\!\(\*SubscriptBox[\(E\), \(\(c\)\(\\ \)\)]\)( in eV)"},BaseStyle->{FontSize->15},PlotStyle->Thick]
jd[3,(10^23),2,2,1,300,1]
Plot[{jd[3,1*(10^22),2,1,1,300,v],jd[3,2*(10^22),2,1,1,300,v],jd[3,3*(10^22),2,1,1,300,v],jd[3,4*(10^22),2,1,1,300,v],jd[3,5*(10^22),2,1,1,300,v]},{v,0,15},Frame->True, FrameLabel->{"V (Volt)","(Subscript[J, 3  d]) (A/(m^2))"},BaseStyle->{FontSize->15},PlotStyle->Thick]

(* 3d Semiconductor *)
(*ηs[d_,nd_,gs_,gv_,meff_,T_,v_]
jd[d_,nd_,gs_,gv_,meff_,T_,v_]*)
EFs[3,(10^20),2,1,0.2,300,1]
Plot[{EFs[3,(10^20),2,1,0.2,300,v],EFs[3,(10^20),2,1,0.2,300,v]-v},{v,0,2},Frame->True, FrameLabel->{"V (Volt)","(Subscript[E, Fs])-(Subscript[E, c]),(Subscript[E, Fd])-\!\(\*SubscriptBox[\(E\), \(\(c\)\(\\ \)\)]\)( in eV)"},BaseStyle->{FontSize->15},PlotStyle->Thick]
jd[3,(10^20),2,2,1,300,1]
Plot[(10^(-4))*{jd[3,1*(10^20),2,1,0.2,300,v],jd[3,2*(10^20),2,1,0.2,300,v],jd[3,3*(10^20),2,1,0.2,300,v],jd[3,4*(10^20),2,1,0.2,300,v],jd[3,5*(10^20),2,1,0.2,300,v]},{v,0,2},Frame->True, FrameLabel->{"V (Volt)","(Subscript[J, 3  d]) (A/(cm^2))"},BaseStyle->{FontSize->15},PlotStyle->Thick]
Plot[(10^(-4))*{jd[3,1*(10^20),2,1,0.2,77,v],jd[3,2*(10^20),2,1,0.2,77,v],jd[3,3*(10^20),2,1,0.2,77,v],jd[3,4*(10^20),2,1,0.2,77,v],jd[3,5*(10^20),2,1,0.2,77,v]},{v,0,2},Frame->True, FrameLabel->{"V (Volt)","(Subscript[J, 3  d]) (A/(cm^2))"},BaseStyle->{FontSize->15},PlotStyle->Thick]

2D Carrier Statistics, Currents, Energies, and Specific Heat

(* Carrier Statistics And Quantum Currents in 1, 2, and 3 Dimensions *)
ClearAll
h:=6.63*(10^(-34)); (* Planck's constant in J.s *)
q:=1.6*(10^(-19));               (* Electron charge, Coulomb *)
hbar:=(6.63/(2π))*(10^(-34));   (* Reduced Planck's constant, J.s *)
kb:=1.38*(10^(-23));           (* Boltzmann constant, J/K *)
me:=9.1*(10^(-31));             (* Electron rest mass, Kg *)
FermiDiracIntegral[d_,η_]:=Switch[d,-1,((E^η)/(1+(E^η))),-(1/2),-PolyLog[(1/2),-(E^η)],0,Log[1+(E^η)],(1/2),-PolyLog[(3/2),-(E^η)],1,-PolyLog[2,-(E^η)],(3/2),-PolyLog[(5/2),-(E^η)],2,-PolyLog[3,-(E^η)],(5/2),-PolyLog[(7/2),-(E^η)],3,-PolyLog[4,-(E^η)]]; (* Define Fermi-Dirac Integrals as PolyLog functions for fast numerical evaluation *)
Nc[d_,gs_,gv_,meff_,T_]:=gs*gv*(((((2*π*meff*me*kb*T)/(h^2))))^(d/2)); (* Band-edge DOS in 1/(m^d) for d dimensions *)
n[d_,gs_,gv_,meff_,T_,Ef_]:=Nc[d,gs,gv,meff,T]*FermiDiracIntegral[((d-2)/2),(Ef/((kb*T)/q))];  (* Carrier density in: 1d in 1/(m^1), 2d in 1/(m^2), and 3d in 1/(m^3) *)
LogPlot[{(10^(-4))*n[2,2,1,1,77,Ef],(10^(-4))*n[2,2,1,1,200,Ef],(10^(-4))*n[2,2,1,1,300,Ef]},{Ef,-1,1},Frame->True, FrameLabel->{"(Subscript[E, F]) (eV)","(Subscript[n, 2  d])(1/(cm^2))"},BaseStyle->{FontSize->15},PlotStyle->{{Thick,Blue},{Thick,Green},{Thick,Red}},AspectRatio->GoldenRatio,PlotRange->{(10^0),(10^15)}]
(* Carrier density in: 1d in 1/(cm^2), LogScale *)
Plot[{(10^(-4))*n[2,2,1,1,77,Ef],(10^(-4))*n[2,2,1,1,200,Ef],(10^(-4))*n[2,2,1,1,300,Ef]},{Ef,-1,1},Frame->True, FrameLabel->{"(Subscript[E, F]) (eV)","(Subscript[n, 2  d])(1/(cm^2))"},BaseStyle->{FontSize->15},PlotStyle->{{Thick,Blue},{Thick,Green},{Thick,Red}},AspectRatio->GoldenRatio,PlotRange->{-(10^13),4*(10^14)}]
(* Carrier density in: 1d in 1/(cm^1), LinearScale *)
EF[d_,nd_,gs_,gv_,meff_,T_]:=Chop[Quiet[FindRoot[FermiDiracIntegral[((d-2)/2),(Efs/((kb*T)/q))]+FermiDiracIntegral[((d-2)/2),(Efs/((kb*T)/q))]==((2*nd)/(Nc[d,gs,gv,meff,T]*(10^(-2*d)))),{Efs,0}]][[1,2]]]; (* nd is in 1/(cm^d) *)
EF[2,1*(10^14),2,1,1,300]
(*Plot[{EF[1,n1d*(10^7),2,1,1,77],EF[1,n1d*(10^7),2,1,1,200],EF[1,n1d*(10^7),2,1,1,300]},{n1d,0,10},Frame->True, FrameLabel->{"(Subscript[n, 1 d])(1/cm)","(Subscript[E, F]) (eV)"},BaseStyle->{FontSize->15},PlotStyle->{{Thick,Blue},{Thick,Green},{Thick,Red}},AspectRatio->GoldenRatio,PlotRange->All]
LogLinearPlot[{EF[1,n1d*(10^7),2,1,1,77],EF[1,n1d*(10^7),2,1,1,200],EF[1,n1d*(10^7),2,1,1,300]},{n1d,(10^(-5)),(10^1)},Frame->True, FrameLabel->{"(Subscript[n, 1 d])(1/cm)","(Subscript[E, F]) (eV)"},BaseStyle->{FontSize->15},PlotStyle->{{Thick,Blue},{Thick,Green},{Thick,Red}},AspectRatio->GoldenRatio,PlotRange->All]*)
ud[d_,T_,η_]:=(d/2)*((kb*T)/q)*((FermiDiracIntegral[(d/2),η])/(FermiDiracIntegral[((d-2)/2),η])); (* Quantum Equipartition of Energy in eV *)
LogPlot[{ud[2,4,((q*Ef)/(kb*4))],ud[2,300,((q*Ef)/(kb*300))],ud[2,1000,((q*Ef)/(kb*1000))]},{Ef,-0.5,1},Frame->True, FrameLabel->{"(Subscript[E, F])(eV)","(Subscript[u, 2  d])(eV)"},BaseStyle->{FontSize->15},PlotStyle->{{Thick,Blue},{Thick,Green},{Thick,Red}},AspectRatio->GoldenRatio,PlotRange->{(10^(-5)),1}]

ud[d_,T_,η_]:=(d/2)*kb*T*((FermiDiracIntegral[(d/2),η])/(FermiDiracIntegral[((d-2)/2),η])); (* Quantum Equipartition of Energy *)

(* Carrier Statistics And Quantum Currents in 1, 2, and 3 Dimensions *)
ClearAll
h:=6.63*(10^(-34)); (* Planck's constant in J.s *)
q:=1.6*(10^(-19));               (* Electron charge, Coulomb *)
hbar:=(6.63/(2π))*(10^(-34));   (* Reduced Planck's constant, J.s *)
kb:=1.38*(10^(-23));           (* Boltzmann constant, J/K *)
me:=9.1*(10^(-31));             (* Electron rest mass, Kg *)
n[d_,gs_,gv_,meff_,T_,Ef_]:=Nc[d,gs,gv,meff,T]*Switch[d,1,-PolyLog[(1/2),-(E^η)],2,Log[1+(E^η)],3,-PolyLog[(3/2),-(E^η)]];  (* Carrier density in 1d, 2d, and 3d in 1/(m^3) *)
Nc[d_,gs_,gv_,meff_,T_]:=gs*gv*(((((2*π*meff*me*kb*T)/(h^2))))^(d/2)); (* Band-edge DOS in 1/(m^d) for d dimensions *)
FermiDiracIntegral[d_,η_]:=Switch[d,-1,((E^η)/(1+(E^η))),-(1/2),-PolyLog[(1/2),-(E^η)],0,Log[1+(E^η)],(1/2),-PolyLog[(3/2),-(E^η)],1,-PolyLog[2,-(E^η)],(3/2),-PolyLog[(5/2),-(E^η)],2,-PolyLog[3,-(E^η)],(5/2),-PolyLog[(7/2),-(E^η)],3,-PolyLog[4,-(E^η)]]; (* Define Fermi-Dirac Integrals as PolyLog functions for fast numerical evaluation *)
EFs[d_,nd_,gs_,gv_,meff_,T_,v_]:=Chop[Quiet[FindRoot[FermiDiracIntegral[((d-2)/2),(Efs/((kb*T)/q))]+FermiDiracIntegral[((d-2)/2),((Efs-v)/((kb*T)/q))]==((2*nd)/(Nc[d,gs,gv,meff,T]*(10^(-2*d)))),{Efs,0}]][[1,2]]]; (* nd is in 1/(cm^d) *)
CurrentPreFactor[d_]:=Switch[d,1,2,2,(1/2),3,((3(Sqrt[π]))/4)]  (* Prefactor for the current density *)
jd[d_,nd_,gs_,gv_,meff_,T_,v_]:=Chop[((q^2)/h)*Nc[d-1,gs,gv,meff,T]*((kb*T)/q)*(FermiDiracIntegral[((d-1)/2),((EFs[d,nd,gs,gv,meff,T,v])/((kb*T)/q))]-FermiDiracIntegral[((d-1)/2),((EFs[d,nd,gs,gv,meff,T,v]-v)/((kb*T)/q))])];
(* Current Density in A/(m^(d-1)) *)

(* n[d_,gs_,gv_,meff_,T_,Ef_] *)
(* Nc[d_,gs_,gv_,meff_,T_] *)
Print["Nc[1,2,1,0.2,300]*(10^(-2))=",Nc[1,2,1,0.2,300]*(10^(-2)),"(cm^(-1))"] (* 1D Band Edge DOS in 1/cm *)
Print["Nc[2,2,1,0.2,300]*(10^(-2))=",Nc[2,2,1,0.2,300]*(10^(-4)),"(cm^(-2))"] (* 2D Band Edge DOS in 1/(cm^2) *)
Print["Nc[3,2,1,0.2,300]*(10^(-2))=",Nc[3,2,1,0.2,300]*(10^(-6)),"(cm^(-3))"] (* 3D Band Edge DOS in 1/(cm^3) *)
LogPlot[n[2,2,1,0.2,300,η]*(10^(-4)),{η,-10,5},Frame->True, FrameLabel->{"\!\(\*FractionBox[\(\*SubscriptBox[\(E\), \(F\)] - \*SubscriptBox[\(E\), \(c\)]\), \(kb\\ T\)]\)","n2d((cm^(-2)))"},BaseStyle->{FontSize->15},PlotStyle->Thick] (* 2D Carrier Density as a function of ((Subscript[E, F])-(Subscript[E, c]))/kb*T *)
Print["jd[1,(10^6),2,1,0.2,300,1]=",jd[1,(10^6),2,1,0.2,300,1],"A"] (* 1D Current in A *)
Print["jd[2,(10^13),2,1,0.2,300,1]=",jd[2,(10^13),2,1,0.2,300,1],"A/m"] (* 2D Current in A *)
Print["jd[3,(10^18),2,1,0.2,300,1]=",jd[3,(10^18),2,1,0.2,300,1],"A/(m^2)"] (* 3D Current in A *)

(* 1d Metal *)
(*ηs[d_,nd_,gs_,gv_,meff_,T_,v_]
jd[d_,nd_,gs_,gv_,meff_,T_,v_]*)
EFs[1,(10^7),2,1,1,300,1]
Plot[{EFs[1,5*(10^7),2,1,1,300,v],EFs[1,5*(10^7),2,1,1,300,v]-v},{v,0,15},Frame->True, FrameLabel->{"V (Volt)","(Subscript[E, Fs])-(Subscript[E, c]),(Subscript[E, Fd])-\!\(\*SubscriptBox[\(E\), \(\(c\)\(\\ \)\)]\)( in eV)"},BaseStyle->{FontSize->15},PlotStyle->Thick]
jd[1,(10^7),2,2,1,300,1]
Plot[{jd[1,(10^7),2,1,1,300,v],jd[1,2*(10^7),2,1,1,300,v],jd[1,3*(10^7),2,1,1,300,v],jd[1,4*(10^7),2,1,1,300,v],jd[1,5*(10^7),2,1,1,300,v]},{v,0,15},Frame->True, FrameLabel->{"V (Volt)","(Subscript[J, 1  d]) (Amp)"},BaseStyle->{FontSize->15},PlotStyle->Thick]

(* 1d Semiconductor *)
(*ηs[d_,nd_,gs_,gv_,meff_,T_,v_]
jd[d_,nd_,gs_,gv_,meff_,T_,v_]*)
EFs[1,(10^6),2,1,0.2,300,1]
Plot[{EFs[1,5*(10^6),2,1,0.2,300,v],EFs[1,5*(10^6),2,1,0.2,300,v]-v},{v,0,2},Frame->True, FrameLabel->{"V (Volt)","(Subscript[E, Fs])-(Subscript[E, c]),(Subscript[E, Fd])-\!\(\*SubscriptBox[\(E\), \(\(c\)\(\\ \)\)]\)( in eV)"},BaseStyle->{FontSize->15},PlotStyle->Thick]
jd[1,(10^6),2,2,0.2,300,1]
Plot[{jd[1,5*(10^6),2,1,0.2,300,v],jd[1,4*(10^6),2,1,0.2,300,v],jd[1,3*(10^6),2,1,0.2,300,v],jd[1,2*(10^6),2,1,0.2,300,v],jd[1,1*(10^6),2,1,0.2,300,v],jd[1,1*(10^5),2,1,0.2,300,v]},{v,0,2},Frame->True, FrameLabel->{"V (Volt)","(Subscript[J, 1  d]) (Amp)"},BaseStyle->{FontSize->15},PlotStyle->Thick]
Plot[{jd[1,5*(10^6),2,1,0.2,77,v],jd[1,4*(10^6),2,1,0.2,77,v],jd[1,3*(10^6),2,1,0.2,77,v],jd[1,2*(10^6),2,1,0.2,77,v],jd[1,1*(10^6),2,1,0.2,77,v],jd[1,1*(10^5),2,1,0.2,77,v]},{v,0,2},Frame->True, FrameLabel->{"V (Volt)","(Subscript[J, 1  d]) (Amp)"},BaseStyle->{FontSize->15},PlotStyle->Thick]

(* 2d Metal *)
(*ηs[d_,nd_,gs_,gv_,meff_,T_,v_]
jd[d_,nd_,gs_,gv_,meff_,T_,v_]*)
EFs[2,(10^14),2,1,1,300,1]
Plot[{EFs[2,5*(10^14),2,1,1,300,0],EFs[2,5*(10^14),2,1,1,300,v],EFs[2,5*(10^14),2,1,1,300,v]-v,EFs[2,(10^13),2,1,1,300,v],EFs[2,(10^13),2,1,1,300,v]-v},{v,0,10},Frame->True, FrameLabel->{"V (Volt)","(Subscript[E, Fs])-(Subscript[E, c]),(Subscript[E, Fd])-\!\(\*SubscriptBox[\(E\), \(\(c\)\(\\ \)\)]\)( in eV)"},BaseStyle->{FontSize->15},PlotStyle->{{Dashed,Red},{Thick,Red},{Thick,Red},{Thick,Blue},{Thick,Blue}},PlotRange->{-6,3}]
jd[2,(10^15),2,2,1,300,1]
Plot[{jd[2,0.1*(10^14),2,1,1,300,v],jd[2,1*(10^14),2,1,1,300,v],jd[2,2*(10^14),2,1,1,300,v],jd[2,3*(10^14),2,1,1,300,v],jd[2,4*(10^14),2,1,1,300,v],jd[2,5*(10^14),2,1,1,300,v]},{v,0,10},Frame->True, FrameLabel->{"V (Volt)","(Subscript[J, 2  d]) (Amp/m)"},BaseStyle->{FontSize->15},PlotStyle->{{Thick,RGBColor[0.085069, 0.157214, 0.922347]},{Thick,RGBColor[0, 0.699168, 0.901411]},{Thick,RGBColor[0, 0.92462, 0.228962]},{Thick,RGBColor[0.907683, 0.742321, 0.19472]},{Thick,RGBColor[0.956741, 0.282475, 0.188861]},{Thick,RGBColor[0.97229, 0, 0.0719463]}}]

(* 2d Semiconductor *)
(*ηs[d_,nd_,gs_,gv_,meff_,T_,v_]
jd[d_,nd_,gs_,gv_,meff_,T_,v_]*)
EFs[2,5*(10^13),2,1,0.2,300,1]
Plot[{EFs[2,5*(10^13),2,1,0.2,300,v],EFs[2,5*(10^13),2,1,0.2,300,v]-v},{v,0,2},Frame->True, FrameLabel->{"V (Volt)","(Subscript[E, Fs])-(Subscript[E, c]),(Subscript[E, Fd])-\!\(\*SubscriptBox[\(E\), \(\(c\)\(\\ \)\)]\)( in eV)"},BaseStyle->{FontSize->15},PlotStyle->Thick]
jd[2,(10^15),2,2,1,300,1]
Plot[(10^(-3))*{jd[2,5*(10^13),2,1,0.2,300,v],jd[2,4*(10^13),2,1,0.2,300,v],jd[2,3*(10^13),2,1,0.2,300,v],jd[2,2*(10^13),2,1,0.2,300,v],jd[2,1*(10^13),2,1,0.2,300,v],jd[2,1*(10^12),2,1,0.2,300,v]},{v,0,2},Frame->True, FrameLabel->{"V (Volt)","(Subscript[J, 2  d]) (mA/μm)"},BaseStyle->{FontSize->15},PlotStyle->Thick]
Plot[(10^(-3))*{jd[2,5*(10^13),2,1,0.2,77,v],jd[2,4*(10^13),2,1,0.2,77,v],jd[2,3*(10^13),2,1,0.2,77,v],jd[2,2*(10^13),2,1,0.2,77,v],jd[2,1*(10^13),2,1,0.2,77,v],jd[2,1*(10^12),2,1,0.2,77,v]},{v,0,2},Frame->True, FrameLabel->{"V (Volt)","(Subscript[J, 2  d]) (mA/μm)"},BaseStyle->{FontSize->15},PlotStyle->Thick]

(* 3d Metal *)
(*ηs[d_,nd_,gs_,gv_,meff_,T_,v_]
jd[d_,nd_,gs_,gv_,meff_,T_,v_]*)
EFs[3,(10^23),2,1,1,300,1]
Plot[{EFs[3,(10^22),2,1,1,300,v],EFs[3,(10^22),2,1,1,300,v]-v},{v,0,15},Frame->True, FrameLabel->{"V (Volt)","(Subscript[E, Fs])-(Subscript[E, c]),(Subscript[E, Fd])-\!\(\*SubscriptBox[\(E\), \(\(c\)\(\\ \)\)]\)( in eV)"},BaseStyle->{FontSize->15},PlotStyle->Thick]
jd[3,(10^23),2,2,1,300,1]
Plot[{jd[3,1*(10^22),2,1,1,300,v],jd[3,2*(10^22),2,1,1,300,v],jd[3,3*(10^22),2,1,1,300,v],jd[3,4*(10^22),2,1,1,300,v],jd[3,5*(10^22),2,1,1,300,v]},{v,0,15},Frame->True, FrameLabel->{"V (Volt)","(Subscript[J, 3  d]) (A/(m^2))"},BaseStyle->{FontSize->15},PlotStyle->Thick]

(* 3d Semiconductor *)
(*ηs[d_,nd_,gs_,gv_,meff_,T_,v_]
jd[d_,nd_,gs_,gv_,meff_,T_,v_]*)
EFs[3,(10^20),2,1,0.2,300,1]
Plot[{EFs[3,(10^20),2,1,0.2,300,v],EFs[3,(10^20),2,1,0.2,300,v]-v},{v,0,2},Frame->True, FrameLabel->{"V (Volt)","(Subscript[E, Fs])-(Subscript[E, c]),(Subscript[E, Fd])-\!\(\*SubscriptBox[\(E\), \(\(c\)\(\\ \)\)]\)( in eV)"},BaseStyle->{FontSize->15},PlotStyle->Thick]
jd[3,(10^20),2,2,1,300,1]
Plot[(10^(-4))*{jd[3,1*(10^20),2,1,0.2,300,v],jd[3,2*(10^20),2,1,0.2,300,v],jd[3,3*(10^20),2,1,0.2,300,v],jd[3,4*(10^20),2,1,0.2,300,v],jd[3,5*(10^20),2,1,0.2,300,v]},{v,0,2},Frame->True, FrameLabel->{"V (Volt)","(Subscript[J, 3  d]) (A/(cm^2))"},BaseStyle->{FontSize->15},PlotStyle->Thick]
Plot[(10^(-4))*{jd[3,1*(10^20),2,1,0.2,77,v],jd[3,2*(10^20),2,1,0.2,77,v],jd[3,3*(10^20),2,1,0.2,77,v],jd[3,4*(10^20),2,1,0.2,77,v],jd[3,5*(10^20),2,1,0.2,77,v]},{v,0,2},Frame->True, FrameLabel->{"V (Volt)","(Subscript[J, 3  d]) (A/(cm^2))"},BaseStyle->{FontSize->15},PlotStyle->Thick]

3D Carrier Statistics, Currents, Energies, and Specific Heat

(* Carrier Statistics And Quantum Currents in 1, 2, and 3 Dimensions *)
ClearAll
h:=6.63*(10^(-34)); (* Planck's constant in J.s *)
q:=1.6*(10^(-19));               (* Electron charge, Coulomb *)
hbar:=(6.63/(2π))*(10^(-34));   (* Reduced Planck's constant, J.s *)
kb:=1.38*(10^(-23));           (* Boltzmann constant, J/K *)
me:=9.1*(10^(-31));             (* Electron rest mass, Kg *)
FermiDiracIntegral[d_,η_]:=Switch[d,-1,((E^η)/(1+(E^η))),-(1/2),-PolyLog[(1/2),-(E^η)],0,Log[1+(E^η)],(1/2),-PolyLog[(3/2),-(E^η)],1,-PolyLog[2,-(E^η)],(3/2),-PolyLog[(5/2),-(E^η)],2,-PolyLog[3,-(E^η)],(5/2),-PolyLog[(7/2),-(E^η)],3,-PolyLog[4,-(E^η)]]; (* Define Fermi-Dirac Integrals as PolyLog functions for fast numerical evaluation *)
Nc[d_,gs_,gv_,meff_,T_]:=gs*gv*(((((2*π*meff*me*kb*T)/(h^2))))^(d/2)); (* Band-edge DOS in 1/(m^d) for d dimensions *)
n[d_,gs_,gv_,meff_,T_,Ef_]:=Nc[d,gs,gv,meff,T]*FermiDiracIntegral[((d-2)/2),(Ef/((kb*T)/q))];  (* Carrier density in: 1d in 1/(m^1), 2d in 1/(m^2), and 3d in 1/(m^3) *)
LogPlot[{(10^(-6))*n[3,2,1,1,77,Ef],(10^(-6))*n[3,2,1,1,200,Ef],(10^(-6))*n[3,2,1,1,300,Ef]},{Ef,-1,1},Frame->True, FrameLabel->{"(Subscript[E, F]) (eV)","(Subscript[n, 3  d])(1/(cm^3))"},BaseStyle->{FontSize->15},PlotStyle->{{Thick,Blue},{Thick,Green},{Thick,Red}},AspectRatio->GoldenRatio,PlotRange->{(10^5),(10^24)}]
(* Carrier density in: 3d in 1/(cm^3), LogScale *)
Plot[{(10^(-6))*n[3,2,1,1,77,Ef],(10^(-6))*n[3,2,1,1,200,Ef],(10^(-6))*n[3,2,1,1,300,Ef]},{Ef,-1,1},Frame->True, FrameLabel->{"(Subscript[E, F]) (eV)","(Subscript[n, 3  d])(1/(cm^3))"},BaseStyle->{FontSize->15},PlotStyle->{{Thick,Blue},{Thick,Green},{Thick,Red}},AspectRatio->GoldenRatio,PlotRange->{-(10^21),(10^22)}]
(* Carrier density in: 3d in 1/(cm^3), LinearScale *)
EF[d_,nd_,gs_,gv_,meff_,T_]:=Chop[Quiet[FindRoot[FermiDiracIntegral[((d-2)/2),(Efs/((kb*T)/q))]+FermiDiracIntegral[((d-2)/2),(Efs/((kb*T)/q))]==((2*nd)/(Nc[d,gs,gv,meff,T]*(10^(-2*d)))),{Efs,0}]][[1,2]]]; (* nd is in 1/(cm^d) *)
EF[3,1*(10^7),2,1,1,300]
(*Plot[{EF[1,n1d*(10^7),2,1,1,77],EF[1,n1d*(10^7),2,1,1,200],EF[1,n1d*(10^7),2,1,1,300]},{n1d,0,10},Frame->True, FrameLabel->{"(Subscript[n, 1 d])(1/cm)","(Subscript[E, F]) (eV)"},BaseStyle->{FontSize->15},PlotStyle->{{Thick,Blue},{Thick,Green},{Thick,Red}},AspectRatio->GoldenRatio,PlotRange->All]
LogLinearPlot[{EF[1,n1d*(10^7),2,1,1,77],EF[1,n1d*(10^7),2,1,1,200],EF[1,n1d*(10^7),2,1,1,300]},{n1d,(10^(-5)),(10^1)},Frame->True, FrameLabel->{"(Subscript[n, 1 d])(1/cm)","(Subscript[E, F]) (eV)"},BaseStyle->{FontSize->15},PlotStyle->{{Thick,Blue},{Thick,Green},{Thick,Red}},AspectRatio->GoldenRatio,PlotRange->All]*)
ud[d_,T_,η_]:=(d/2)*((kb*T)/q)*((FermiDiracIntegral[(d/2),η])/(FermiDiracIntegral[((d-2)/2),η])); (* Quantum Equipartition of Energy in eV *)
LogPlot[{ud[3,4,((q*Ef)/(kb*4))],ud[3,300,((q*Ef)/(kb*300))],ud[3,1000,((q*Ef)/(kb*1000))]},{Ef,-0.5,1},Frame->True, FrameLabel->{"(Subscript[E, F])(eV)","(Subscript[u, 3  d])(eV)"},BaseStyle->{FontSize->15},PlotStyle->{{Thick,Blue},{Thick,Green},{Thick,Red}},AspectRatio->GoldenRatio,PlotRange->{(10^(-5)),1}]

ud[d_,T_,η_]:=(d/2)*kb*T*((FermiDiracIntegral[(d/2),η])/(FermiDiracIntegral[((d-2)/2),η])); (* Quantum Equipartition of Energy *)

(* Carrier Statistics And Quantum Currents in 1, 2, and 3 Dimensions *)
ClearAll
h:=6.63*(10^(-34)); (* Planck's constant in J.s *)
q:=1.6*(10^(-19));               (* Electron charge, Coulomb *)
hbar:=(6.63/(2π))*(10^(-34));   (* Reduced Planck's constant, J.s *)
kb:=1.38*(10^(-23));           (* Boltzmann constant, J/K *)
me:=9.1*(10^(-31));             (* Electron rest mass, Kg *)
n[d_,gs_,gv_,meff_,T_,Ef_]:=Nc[d,gs,gv,meff,T]*Switch[d,1,-PolyLog[(1/2),-(E^η)],2,Log[1+(E^η)],3,-PolyLog[(3/2),-(E^η)]];  (* Carrier density in 1d, 2d, and 3d in 1/(m^3) *)
Nc[d_,gs_,gv_,meff_,T_]:=gs*gv*(((((2*π*meff*me*kb*T)/(h^2))))^(d/2)); (* Band-edge DOS in 1/(m^d) for d dimensions *)
FermiDiracIntegral[d_,η_]:=Switch[d,-1,((E^η)/(1+(E^η))),-(1/2),-PolyLog[(1/2),-(E^η)],0,Log[1+(E^η)],(1/2),-PolyLog[(3/2),-(E^η)],1,-PolyLog[2,-(E^η)],(3/2),-PolyLog[(5/2),-(E^η)],2,-PolyLog[3,-(E^η)],(5/2),-PolyLog[(7/2),-(E^η)],3,-PolyLog[4,-(E^η)]]; (* Define Fermi-Dirac Integrals as PolyLog functions for fast numerical evaluation *)
EFs[d_,nd_,gs_,gv_,meff_,T_,v_]:=Chop[Quiet[FindRoot[FermiDiracIntegral[((d-2)/2),(Efs/((kb*T)/q))]+FermiDiracIntegral[((d-2)/2),((Efs-v)/((kb*T)/q))]==((2*nd)/(Nc[d,gs,gv,meff,T]*(10^(-2*d)))),{Efs,0}]][[1,2]]]; (* nd is in 1/(cm^d) *)
CurrentPreFactor[d_]:=Switch[d,1,2,2,(1/2),3,((3(Sqrt[π]))/4)]  (* Prefactor for the current density *)
jd[d_,nd_,gs_,gv_,meff_,T_,v_]:=Chop[((q^2)/h)*Nc[d-1,gs,gv,meff,T]*((kb*T)/q)*(FermiDiracIntegral[((d-1)/2),((EFs[d,nd,gs,gv,meff,T,v])/((kb*T)/q))]-FermiDiracIntegral[((d-1)/2),((EFs[d,nd,gs,gv,meff,T,v]-v)/((kb*T)/q))])];
(* Current Density in A/(m^(d-1)) *)

(* n[d_,gs_,gv_,meff_,T_,Ef_] *)
(* Nc[d_,gs_,gv_,meff_,T_] *)
Print["Nc[1,2,1,0.2,300]*(10^(-2))=",Nc[1,2,1,0.2,300]*(10^(-2)),"(cm^(-1))"] (* 1D Band Edge DOS in 1/cm *)
Print["Nc[2,2,1,0.2,300]*(10^(-2))=",Nc[2,2,1,0.2,300]*(10^(-4)),"(cm^(-2))"] (* 2D Band Edge DOS in 1/(cm^2) *)
Print["Nc[3,2,1,0.2,300]*(10^(-2))=",Nc[3,2,1,0.2,300]*(10^(-6)),"(cm^(-3))"] (* 3D Band Edge DOS in 1/(cm^3) *)
LogPlot[n[2,2,1,0.2,300,η]*(10^(-4)),{η,-10,5},Frame->True, FrameLabel->{"\!\(\*FractionBox[\(\*SubscriptBox[\(E\), \(F\)] - \*SubscriptBox[\(E\), \(c\)]\), \(kb\\ T\)]\)","n2d((cm^(-2)))"},BaseStyle->{FontSize->15},PlotStyle->Thick] (* 2D Carrier Density as a function of ((Subscript[E, F])-(Subscript[E, c]))/kb*T *)
Print["jd[1,(10^6),2,1,0.2,300,1]=",jd[1,(10^6),2,1,0.2,300,1],"A"] (* 1D Current in A *)
Print["jd[2,(10^13),2,1,0.2,300,1]=",jd[2,(10^13),2,1,0.2,300,1],"A/m"] (* 2D Current in A *)
Print["jd[3,(10^18),2,1,0.2,300,1]=",jd[3,(10^18),2,1,0.2,300,1],"A/(m^2)"] (* 3D Current in A *)

(* 1d Metal *)
(*ηs[d_,nd_,gs_,gv_,meff_,T_,v_]
jd[d_,nd_,gs_,gv_,meff_,T_,v_]*)
EFs[1,(10^7),2,1,1,300,1]
Plot[{EFs[1,5*(10^7),2,1,1,300,v],EFs[1,5*(10^7),2,1,1,300,v]-v},{v,0,15},Frame->True, FrameLabel->{"V (Volt)","(Subscript[E, Fs])-(Subscript[E, c]),(Subscript[E, Fd])-\!\(\*SubscriptBox[\(E\), \(\(c\)\(\\ \)\)]\)( in eV)"},BaseStyle->{FontSize->15},PlotStyle->Thick]
jd[1,(10^7),2,2,1,300,1]
Plot[{jd[1,(10^7),2,1,1,300,v],jd[1,2*(10^7),2,1,1,300,v],jd[1,3*(10^7),2,1,1,300,v],jd[1,4*(10^7),2,1,1,300,v],jd[1,5*(10^7),2,1,1,300,v]},{v,0,15},Frame->True, FrameLabel->{"V (Volt)","(Subscript[J, 1  d]) (Amp)"},BaseStyle->{FontSize->15},PlotStyle->Thick]

(* 1d Semiconductor *)
(*ηs[d_,nd_,gs_,gv_,meff_,T_,v_]
jd[d_,nd_,gs_,gv_,meff_,T_,v_]*)
EFs[1,(10^6),2,1,0.2,300,1]
Plot[{EFs[1,5*(10^6),2,1,0.2,300,v],EFs[1,5*(10^6),2,1,0.2,300,v]-v},{v,0,2},Frame->True, FrameLabel->{"V (Volt)","(Subscript[E, Fs])-(Subscript[E, c]),(Subscript[E, Fd])-\!\(\*SubscriptBox[\(E\), \(\(c\)\(\\ \)\)]\)( in eV)"},BaseStyle->{FontSize->15},PlotStyle->Thick]
jd[1,(10^6),2,2,0.2,300,1]
Plot[{jd[1,5*(10^6),2,1,0.2,300,v],jd[1,4*(10^6),2,1,0.2,300,v],jd[1,3*(10^6),2,1,0.2,300,v],jd[1,2*(10^6),2,1,0.2,300,v],jd[1,1*(10^6),2,1,0.2,300,v],jd[1,1*(10^5),2,1,0.2,300,v]},{v,0,2},Frame->True, FrameLabel->{"V (Volt)","(Subscript[J, 1  d]) (Amp)"},BaseStyle->{FontSize->15},PlotStyle->Thick]
Plot[{jd[1,5*(10^6),2,1,0.2,77,v],jd[1,4*(10^6),2,1,0.2,77,v],jd[1,3*(10^6),2,1,0.2,77,v],jd[1,2*(10^6),2,1,0.2,77,v],jd[1,1*(10^6),2,1,0.2,77,v],jd[1,1*(10^5),2,1,0.2,77,v]},{v,0,2},Frame->True, FrameLabel->{"V (Volt)","(Subscript[J, 1  d]) (Amp)"},BaseStyle->{FontSize->15},PlotStyle->Thick]

(* 2d Metal *)
(*ηs[d_,nd_,gs_,gv_,meff_,T_,v_]
jd[d_,nd_,gs_,gv_,meff_,T_,v_]*)
EFs[2,(10^14),2,1,1,300,1]
Plot[{EFs[2,(10^15),2,1,1,300,v],EFs[2,(10^15),2,1,1,300,v]-v},{v,0,10},Frame->True, FrameLabel->{"V (Volt)","(Subscript[E, Fs])-(Subscript[E, c]),(Subscript[E, Fd])-\!\(\*SubscriptBox[\(E\), \(\(c\)\(\\ \)\)]\)( in eV)"},BaseStyle->{FontSize->15},PlotStyle->Thick]
jd[2,(10^15),2,2,1,300,1]
Plot[{jd[2,1*(10^14),2,1,1,300,v],jd[2,2*(10^14),2,1,1,300,v],jd[2,4*(10^14),2,1,1,300,v],jd[2,8*(10^14),2,1,1,300,v],jd[2,(10^15),2,1,1,300,v]},{v,0,15},Frame->True, FrameLabel->{"V (Volt)","(Subscript[J, 2  d]) (Amp/m)"},BaseStyle->{FontSize->15},PlotStyle->Thick]

(* 2d Semiconductor *)
(*ηs[d_,nd_,gs_,gv_,meff_,T_,v_]
jd[d_,nd_,gs_,gv_,meff_,T_,v_]*)
EFs[2,5*(10^13),2,1,0.2,300,1]
Plot[{EFs[2,5*(10^13),2,1,0.2,300,v],EFs[2,5*(10^13),2,1,0.2,300,v]-v},{v,0,2},Frame->True, FrameLabel->{"V (Volt)","(Subscript[E, Fs])-(Subscript[E, c]),(Subscript[E, Fd])-\!\(\*SubscriptBox[\(E\), \(\(c\)\(\\ \)\)]\)( in eV)"},BaseStyle->{FontSize->15},PlotStyle->Thick]
jd[2,(10^15),2,2,1,300,1]
Plot[(10^(-3))*{jd[2,5*(10^13),2,1,0.2,300,v],jd[2,4*(10^13),2,1,0.2,300,v],jd[2,3*(10^13),2,1,0.2,300,v],jd[2,2*(10^13),2,1,0.2,300,v],jd[2,1*(10^13),2,1,0.2,300,v],jd[2,1*(10^12),2,1,0.2,300,v]},{v,0,2},Frame->True, FrameLabel->{"V (Volt)","(Subscript[J, 2  d]) (mA/μm)"},BaseStyle->{FontSize->15},PlotStyle->Thick]
Plot[(10^(-3))*{jd[2,5*(10^13),2,1,0.2,77,v],jd[2,4*(10^13),2,1,0.2,77,v],jd[2,3*(10^13),2,1,0.2,77,v],jd[2,2*(10^13),2,1,0.2,77,v],jd[2,1*(10^13),2,1,0.2,77,v],jd[2,1*(10^12),2,1,0.2,77,v]},{v,0,2},Frame->True, FrameLabel->{"V (Volt)","(Subscript[J, 2  d]) (mA/μm)"},BaseStyle->{FontSize->15},PlotStyle->Thick]

(* 3d Metal *)
(*ηs[d_,nd_,gs_,gv_,meff_,T_,v_]
jd[d_,nd_,gs_,gv_,meff_,T_,v_]*)
EFs[3,(10^23),2,1,1,300,1]
Plot[{EFs[3,0.1*(10^22),2,1,1,300,0],EFs[3,0.1*(10^22),2,1,1,300,v],EFs[3,0.1*(10^22),2,1,1,300,v]-v,EFs[3,(10^22),2,1,1,300,0],EFs[3,(10^22),2,1,1,300,v],EFs[3,(10^22),2,1,1,300,v]-v},{v,0,10},Frame->True, FrameLabel->{"V (Volt)","(Subscript[E, Fs])-(Subscript[E, c]),(Subscript[E, Fd])-\!\(\*SubscriptBox[\(E\), \(\(c\)\(\\ \)\)]\)( in eV)"},BaseStyle->{FontSize->15},PlotStyle->{{Dashed,Blue},{Thick,Blue},{Thick,Blue},{Dashed,Red},{Thick,Red},{Thick,Red}}]
jd[3,(10^23),2,2,1,300,1]
Plot[{jd[3,0.1*(10^22),2,1,1,300,v],jd[3,1*(10^22),2,1,1,300,v],jd[3,2*(10^22),2,1,1,300,v],jd[3,3*(10^22),2,1,1,300,v],jd[3,4*(10^22),2,1,1,300,v],jd[3,5*(10^22),2,1,1,300,v]},{v,0,10},Frame->True, FrameLabel->{"V (Volt)","(Subscript[J, 3  d]) (A/(m^2))"},BaseStyle->{FontSize->15},PlotStyle->{{Thick,RGBColor[0.085069, 0.157214, 0.922347]},{Thick,RGBColor[0, 0.699168, 0.901411]},{Thick,RGBColor[0, 0.92462, 0.228962]},{Thick,RGBColor[0.907683, 0.742321, 0.19472]},{Thick,RGBColor[0.956741, 0.282475, 0.188861]},{Thick,RGBColor[0.97229, 0, 0.0719463]}}]

(* 3d Semiconductor *)
(*ηs[d_,nd_,gs_,gv_,meff_,T_,v_]
jd[d_,nd_,gs_,gv_,meff_,T_,v_]*)
EFs[3,(10^20),2,1,0.2,300,1]
Plot[{EFs[3,(10^20),2,1,0.2,300,v],EFs[3,(10^20),2,1,0.2,300,v]-v},{v,0,2},Frame->True, FrameLabel->{"V (Volt)","(Subscript[E, Fs])-(Subscript[E, c]),(Subscript[E, Fd])-\!\(\*SubscriptBox[\(E\), \(\(c\)\(\\ \)\)]\)( in eV)"},BaseStyle->{FontSize->15},PlotStyle->Thick]
jd[3,(10^20),2,2,1,300,1]
Plot[(10^(-4))*{jd[3,1*(10^20),2,1,0.2,300,v],jd[3,2*(10^20),2,1,0.2,300,v],jd[3,3*(10^20),2,1,0.2,300,v],jd[3,4*(10^20),2,1,0.2,300,v],jd[3,5*(10^20),2,1,0.2,300,v]},{v,0,2},Frame->True, FrameLabel->{"V (Volt)","(Subscript[J, 3  d]) (A/(cm^2))"},BaseStyle->{FontSize->15},PlotStyle->Thick]
Plot[(10^(-4))*{jd[3,1*(10^20),2,1,0.2,77,v],jd[3,2*(10^20),2,1,0.2,77,v],jd[3,3*(10^20),2,1,0.2,77,v],jd[3,4*(10^20),2,1,0.2,77,v],jd[3,5*(10^20),2,1,0.2,77,v]},{v,0,2},Frame->True, FrameLabel->{"V (Volt)","(Subscript[J, 3  d]) (A/(cm^2))"},BaseStyle->{FontSize->15},PlotStyle->Thick]

New

(* Fermi Level variation for a 2D Electron Gas, n2d in /(cm^2) *)
EF[gs_,gv_,n2d_,T_]:=(((1.38*(10^(-23))*T)/(1.6*(10^(-19)))))*Log[Exp[((n2d*(10^4))/((((gs*gv*9.1*(10^(-31)))/(2*π*(((6.63*(10^(-34))))^2)))*1.38*(10^(-23))*T)))]-1];
LogLogPlot[{EF[2,1,n2d,4],EF[2,1,n2d,77],EF[2,1,n2d,300]},{n2d,(10^10),(10^15)},PlotRange->{(10^(-4)),(10^1)},Frame->True,FrameLabel->{"(Subscript[n, 2  d])(/(cm^2))","(Subscript[E, F])(eV)"},PlotStyle->{Blue,Green,Red},BaseStyle->{FontSize->15},AspectRatio->1,GridLines->Automatic]

FermiDiracIntegral[d_,η_]:=Switch[d,-1,((E^η)/(1+(E^η))),-(1/2),-PolyLog[(1/2),-(E^η)],0,Log[1+(E^η)],(1/2),-PolyLog[(3/2),-(E^η)],1,-PolyLog[2,-(E^η)],(3/2),-PolyLog[(5/2),-(E^η)],2,-PolyLog[3,-(E^η)],(5/2),-PolyLog[(7/2),-(E^η)],3,-PolyLog[4,-(E^η)]]; (* Define Fermi-Dirac Integrals as PolyLog functions for fast numerical evaluation *)
FindRoot[FermiDiracIntegral[-(1/2),x]+FermiDiracIntegral[-(1/2),x-10]==20,{x,0}](* Finding the Root of Fermi-Dirac Sums or Differences *)
FindRoot[FermiDiracIntegral[-(1/2),x]+FermiDiracIntegral[-(1/2),x-10]==0.8,{x,0}] (* Finding the Root of Fermi-Dirac Sums or Differences *)
Chop[FindRoot[FermiDiracIntegral[-(1/2),x]+FermiDiracIntegral[-(1/2),x-10]==0.8,{x,0}]][[1,2]](* Chop out the imaginary Residues *)
Plot[Chop[FindRoot[FermiDiracIntegral[-(1/2),x]+FermiDiracIntegral[-(1/2),x-v]==0.8,{x,0}]][[1,2]],{v,0,100}](* Plot the Solutions *)

Plot[Chop[FindRoot[-PolyLog[0.5,-Exp[x]]-PolyLog[0.5,-Exp[x-v]]==0.76,{x,0}][[1,2]]],{v,0,20},PlotRange->All]

Plot[{Chop[FindRoot[-PolyLog[(1/2),-Exp[x]]-PolyLog[(1/2),-Exp[x-v]]==0.76,{x,0}]][[1,2]],Chop[FindRoot[-PolyLog[(1/2),-Exp[x]]-PolyLog[(1/2),-Exp[x-v]]==0.76,{x,0}]][[1,2]]-v},{v,0,20},PlotRange->All]

(* Carrier Statistics And Quantum Currents in 1, 2, and 3 Dimensions *)
ClearAll
h:=6.63*(10^(-34)); (* Planck's constant in J.s *)
q:=1.6*(10^(-19));               (* Electron charge, Coulomb *)
hbar:=(6.63/(2π))*(10^(-34));   (* Reduced Planck's constant, J.s *)
kb:=1.38*(10^(-23));           (* Boltzmann constant, J/K *)
me:=9.1*(10^(-31));             (* Electron rest mass, Kg *)
n[d_,gs_,gv_,meff_,T_,Ef_]:=Nc[d,gs,gv,meff,T]*Switch[d,1,-PolyLog[(1/2),-(E^η)],2,Log[1+(E^η)],3,-PolyLog[(3/2),-(E^η)]];  (* Carrier density in 1d, 2d, and 3d in 1/(m^3) *)
Nc[d_,gs_,gv_,meff_,T_]:=gs*gv*(((((2*π*meff*me*kb*T)/(h^2))))^(d/2)); (* Band-edge DOS in 1/(m^d) for d dimensions *)
FermiDiracIntegral[d_,η_]:=Switch[d,-1,((E^η)/(1+(E^η))),-(1/2),-PolyLog[(1/2),-(E^η)],0,Log[1+(E^η)],(1/2),-PolyLog[(3/2),-(E^η)],1,-PolyLog[2,-(E^η)],(3/2),-PolyLog[(5/2),-(E^η)],2,-PolyLog[3,-(E^η)],(5/2),-PolyLog[(7/2),-(E^η)],3,-PolyLog[4,-(E^η)]]; (* Define Fermi-Dirac Integrals as PolyLog functions for fast numerical evaluation *)
EFs[d_,nd_,gs_,gv_,meff_,T_,v_]:=Chop[Quiet[FindRoot[FermiDiracIntegral[((d-2)/2),(Efs/((kb*T)/q))]+FermiDiracIntegral[((d-2)/2),((Efs-v)/((kb*T)/q))]==((2*nd)/(Nc[d,gs,gv,meff,T]*(10^(-2*d)))),{Efs,0}]][[1,2]]]; (* nd is in 1/(cm^d) *)
CurrentPreFactor[d_]:=Switch[d,1,2,2,(1/2),3,((3(Sqrt[π]))/4)]  (* Prefactor for the current density *)
jd[d_,nd_,gs_,gv_,meff_,T_,v_]:=Chop[((q^2)/h)*Nc[d-1,gs,gv,meff,T]*((kb*T)/q)*(FermiDiracIntegral[((d-1)/2),((EFs[d,nd,gs,gv,meff,T,v])/((kb*T)/q))]-FermiDiracIntegral[((d-1)/2),((EFs[d,nd,gs,gv,meff,T,v]-v)/((kb*T)/q))])];
(* Current Density in A/(m^(d-1)) *)

(* n[d_,gs_,gv_,meff_,T_,Ef_] *)
(* Nc[d_,gs_,gv_,meff_,T_] *)
Print["Nc[1,2,1,0.2,300]*(10^(-2))=",Nc[1,2,1,0.2,300]*(10^(-2)),"(cm^(-1))"] (* 1D Band Edge DOS in 1/cm *)
Print["Nc[2,2,1,0.2,300]*(10^(-2))=",Nc[2,2,1,0.2,300]*(10^(-4)),"(cm^(-2))"] (* 2D Band Edge DOS in 1/(cm^2) *)
Print["Nc[3,2,1,0.2,300]*(10^(-2))=",Nc[3,2,1,0.2,300]*(10^(-6)),"(cm^(-3))"] (* 3D Band Edge DOS in 1/(cm^3) *)
LogPlot[n[2,2,1,0.2,300,η]*(10^(-4)),{η,-10,5},Frame->True, FrameLabel->{"\!\(\*FractionBox[\(\*SubscriptBox[\(E\), \(F\)] - \*SubscriptBox[\(E\), \(c\)]\), \(kb\\ T\)]\)","n2d((cm^(-2)))"},BaseStyle->{FontSize->15},PlotStyle->Thick] (* 2D Carrier Density as a function of ((Subscript[E, F])-(Subscript[E, c]))/kb*T *)
Print["jd[1,(10^6),2,1,0.2,300,1]=",jd[1,(10^6),2,1,0.2,300,1],"A"] (* 1D Current in A *)
Print["jd[2,(10^13),2,1,0.2,300,1]=",jd[2,(10^13),2,1,0.2,300,1],"A/m"] (* 2D Current in A *)
Print["jd[3,(10^18),2,1,0.2,300,1]=",jd[3,(10^18),2,1,0.2,300,1],"A/(m^2)"] (* 3D Current in A *)

(* 1d Metal *)
(*ηs[d_,nd_,gs_,gv_,meff_,T_,v_]
jd[d_,nd_,gs_,gv_,meff_,T_,v_]*)
EFs[1,(10^7),2,1,1,300,1]
Plot[{EFs[1,5*(10^7),2,1,1,300,v],EFs[1,5*(10^7),2,1,1,300,v]-v},{v,0,15},Frame->True, FrameLabel->{"V (Volt)","(Subscript[E, Fs])-(Subscript[E, c]),(Subscript[E, Fd])-\!\(\*SubscriptBox[\(E\), \(\(c\)\(\\ \)\)]\)( in eV)"},BaseStyle->{FontSize->15},PlotStyle->Thick]
jd[1,(10^7),2,2,1,300,1]
Plot[{jd[1,(10^7),2,1,1,300,v],jd[1,2*(10^7),2,1,1,300,v],jd[1,3*(10^7),2,1,1,300,v],jd[1,4*(10^7),2,1,1,300,v],jd[1,5*(10^7),2,1,1,300,v]},{v,0,15},Frame->True, FrameLabel->{"V (Volt)","(Subscript[J, 1  d]) (Amp)"},BaseStyle->{FontSize->15},PlotStyle->Thick]

(* 1d Semiconductor *)
(*ηs[d_,nd_,gs_,gv_,meff_,T_,v_]
jd[d_,nd_,gs_,gv_,meff_,T_,v_]*)
EFs[1,(10^6),2,1,0.2,300,1]
Plot[{EFs[1,5*(10^6),2,1,0.2,300,v],EFs[1,5*(10^6),2,1,0.2,300,v]-v},{v,0,2},Frame->True, FrameLabel->{"V (Volt)","(Subscript[E, Fs])-(Subscript[E, c]),(Subscript[E, Fd])-\!\(\*SubscriptBox[\(E\), \(\(c\)\(\\ \)\)]\)( in eV)"},BaseStyle->{FontSize->15},PlotStyle->Thick]
jd[1,(10^6),2,2,0.2,300,1]
Plot[{jd[1,5*(10^6),2,1,0.2,300,v],jd[1,4*(10^6),2,1,0.2,300,v],jd[1,3*(10^6),2,1,0.2,300,v],jd[1,2*(10^6),2,1,0.2,300,v],jd[1,1*(10^6),2,1,0.2,300,v],jd[1,1*(10^5),2,1,0.2,300,v]},{v,0,2},Frame->True, FrameLabel->{"V (Volt)","(Subscript[J, 1  d]) (Amp)"},BaseStyle->{FontSize->15},PlotStyle->Thick]
Plot[{jd[1,5*(10^6),2,1,0.2,77,v],jd[1,4*(10^6),2,1,0.2,77,v],jd[1,3*(10^6),2,1,0.2,77,v],jd[1,2*(10^6),2,1,0.2,77,v],jd[1,1*(10^6),2,1,0.2,77,v],jd[1,1*(10^5),2,1,0.2,77,v]},{v,0,2},Frame->True, FrameLabel->{"V (Volt)","(Subscript[J, 1  d]) (Amp)"},BaseStyle->{FontSize->15},PlotStyle->Thick]

(* 2d Metal *)
(*ηs[d_,nd_,gs_,gv_,meff_,T_,v_]
jd[d_,nd_,gs_,gv_,meff_,T_,v_]*)
EFs[2,(10^14),2,1,1,300,1]
Plot[{EFs[2,(10^15),2,1,1,300,v],EFs[2,(10^15),2,1,1,300,v]-v},{v,0,10},Frame->True, FrameLabel->{"V (Volt)","(Subscript[E, Fs])-(Subscript[E, c]),(Subscript[E, Fd])-\!\(\*SubscriptBox[\(E\), \(\(c\)\(\\ \)\)]\)( in eV)"},BaseStyle->{FontSize->15},PlotStyle->Thick]
jd[2,(10^15),2,2,1,300,1]
Plot[{jd[2,1*(10^14),2,1,1,300,v],jd[2,2*(10^14),2,1,1,300,v],jd[2,4*(10^14),2,1,1,300,v],jd[2,8*(10^14),2,1,1,300,v],jd[2,(10^15),2,1,1,300,v]},{v,0,15},Frame->True, FrameLabel->{"V (Volt)","(Subscript[J, 2  d]) (Amp/m)"},BaseStyle->{FontSize->15},PlotStyle->Thick]

(* 2d Semiconductor *)
(*ηs[d_,nd_,gs_,gv_,meff_,T_,v_]
jd[d_,nd_,gs_,gv_,meff_,T_,v_]*)
EFs[2,5*(10^13),2,1,0.2,300,1]
Plot[{EFs[2,5*(10^13),2,1,0.2,300,v],EFs[2,5*(10^13),2,1,0.2,300,v]-v},{v,0,2},Frame->True, FrameLabel->{"V (Volt)","(Subscript[E, Fs])-(Subscript[E, c]),(Subscript[E, Fd])-\!\(\*SubscriptBox[\(E\), \(\(c\)\(\\ \)\)]\)( in eV)"},BaseStyle->{FontSize->15},PlotStyle->Thick]
jd[2,(10^15),2,2,1,300,1]
Plot[(10^(-3))*{jd[2,5*(10^13),2,1,0.2,300,v],jd[2,4*(10^13),2,1,0.2,300,v],jd[2,3*(10^13),2,1,0.2,300,v],jd[2,2*(10^13),2,1,0.2,300,v],jd[2,1*(10^13),2,1,0.2,300,v],jd[2,1*(10^12),2,1,0.2,300,v]},{v,0,2},Frame->True, FrameLabel->{"V (Volt)","(Subscript[J, 2  d]) (mA/μm)"},BaseStyle->{FontSize->15},PlotStyle->Thick]
Plot[(10^(-3))*{jd[2,5*(10^13),2,1,0.2,77,v],jd[2,4*(10^13),2,1,0.2,77,v],jd[2,3*(10^13),2,1,0.2,77,v],jd[2,2*(10^13),2,1,0.2,77,v],jd[2,1*(10^13),2,1,0.2,77,v],jd[2,1*(10^12),2,1,0.2,77,v]},{v,0,2},Frame->True, FrameLabel->{"V (Volt)","(Subscript[J, 2  d]) (mA/μm)"},BaseStyle->{FontSize->15},PlotStyle->Thick]

(* 3d Metal *)
(*ηs[d_,nd_,gs_,gv_,meff_,T_,v_]
jd[d_,nd_,gs_,gv_,meff_,T_,v_]*)
EFs[3,(10^23),2,1,1,300,1]
Plot[{EFs[3,(10^22),2,1,1,300,v],EFs[3,(10^22),2,1,1,300,v]-v},{v,0,15},Frame->True, FrameLabel->{"V (Volt)","(Subscript[E, Fs])-(Subscript[E, c]),(Subscript[E, Fd])-\!\(\*SubscriptBox[\(E\), \(\(c\)\(\\ \)\)]\)( in eV)"},BaseStyle->{FontSize->15},PlotStyle->Thick]
jd[3,(10^23),2,2,1,300,1]
Plot[{jd[3,1*(10^22),2,1,1,300,v],jd[3,2*(10^22),2,1,1,300,v],jd[3,3*(10^22),2,1,1,300,v],jd[3,4*(10^22),2,1,1,300,v],jd[3,5*(10^22),2,1,1,300,v]},{v,0,15},Frame->True, FrameLabel->{"V (Volt)","(Subscript[J, 3  d]) (A/(m^2))"},BaseStyle->{FontSize->15},PlotStyle->Thick]

(* 3d Semiconductor *)
(*ηs[d_,nd_,gs_,gv_,meff_,T_,v_]
jd[d_,nd_,gs_,gv_,meff_,T_,v_]*)
EFs[3,(10^20),2,1,0.2,300,1]
Plot[{EFs[3,(10^20),2,1,0.2,300,v],EFs[3,(10^20),2,1,0.2,300,v]-v},{v,0,2},Frame->True, FrameLabel->{"V (Volt)","(Subscript[E, Fs])-(Subscript[E, c]),(Subscript[E, Fd])-\!\(\*SubscriptBox[\(E\), \(\(c\)\(\\ \)\)]\)( in eV)"},BaseStyle->{FontSize->15},PlotStyle->Thick]
jd[3,(10^20),2,2,1,300,1]
Plot[(10^(-4))*{jd[3,1*(10^20),2,1,0.2,300,v],jd[3,2*(10^20),2,1,0.2,300,v],jd[3,3*(10^20),2,1,0.2,300,v],jd[3,4*(10^20),2,1,0.2,300,v],jd[3,5*(10^20),2,1,0.2,300,v]},{v,0,2},Frame->True, FrameLabel->{"V (Volt)","(Subscript[J, 3  d]) (A/(cm^2))"},BaseStyle->{FontSize->15},PlotStyle->Thick]
Plot[(10^(-4))*{jd[3,1*(10^20),2,1,0.2,77,v],jd[3,2*(10^20),2,1,0.2,77,v],jd[3,3*(10^20),2,1,0.2,77,v],jd[3,4*(10^20),2,1,0.2,77,v],jd[3,5*(10^20),2,1,0.2,77,v]},{v,0,2},Frame->True, FrameLabel->{"V (Volt)","(Subscript[J, 3  d]) (A/(cm^2))"},BaseStyle->{FontSize->15},PlotStyle->Thick]

Old--->>>

First

(* Free Electron Energy Bandstructure *)
h:=6.63*(10^(-34)); (* J.s, Planck's constant *)
hbar:=((6.63*(10^(-34)))/(2π)); (* J.s, Reduced Planck's constant *)
q:=1.6*(10^(-19 ));  (* Coulomb, electron charge *)
m0:=9.1*(10^(-31)); (* kg, free electron rest mass *)
freeenergy[k_]:=(((hbar^2)*(((k*(10^9)))^2))/(2*m0*q));  (* particle in a box energies, in eV units. Note k in (nm^(-1)) *)
Plot[{freeenergy[k]},{k,-10,10},Frame->True,FrameLabel->{"k ((nm^(-1)))","E(k) (eV)"},PlotStyle->Thick,BaseStyle->{FontSize->15},AspectRatio->GoldenRatio]  (* Note the energy and wavevector scales *)

(* Compton wavelength, and the de Broglie relation for free electrons *)
c=3*(10^8);  (* Speed of light, m/s *)
λc=((6.63*(10^(-34)))/(9.1*(10^(-31))*c));  (* Compton wavelength of electron *)
λc
v[λ_]:=(c/(Sqrt[1+((((λ/(λc))))^2)]));  (* relativistic velocity *)
LogLogPlot[{v[0],v[λ]},{λ,(10^(-15)),(10^(-6))},Frame->True,FrameLabel->{"λ (m)","v (m/s)"},PlotStyle->{Dashed,Thick},BaseStyle->{FontSize->15},AspectRatio->GoldenRatio,GridLines->Automatic]

(* Fermi Level variation for a 2D Electron Gas, n2d in /(cm^2) *)
EF[gs_,gv_,n2d_,T_]:=(((1.38*(10^(-23))*T)/(1.6*(10^(-19)))))*Log[Exp[((n2d*(10^4))/((((gs*gv*9.1*(10^(-31)))/(2*π*(((6.63*(10^(-34))))^2)))*1.38*(10^(-23))*T)))]-1];
LogLogPlot[{EF[2,1,n2d,4],EF[2,1,n2d,77],EF[2,1,n2d,300]},{n2d,(10^10),(10^15)},PlotRange->{(10^(-4)),(10^1)},Frame->True,FrameLabel->{"(Subscript[n, 2  d])(/(cm^2))","(Subscript[E, F])(eV)"},PlotStyle->{Blue,Green,Red},BaseStyle->{FontSize->15},AspectRatio->1,GridLines->Automatic]

(* Group Velocities of 2D Free Electrons *)
VectorPlot[{kx,ky},{kx,-3,3},{ky,-3,3},FrameLabel->{"(Subscript[k, x])","(Subscript[k, y])"},BaseStyle->{FontSize->15}]
VectorPlot[{y,-x},{x,-3,3},{y,-3,3},FrameLabel->{"(Subscript[k, x])","(Subscript[k, y])"},BaseStyle->{FontSize->15}]

(* Defining the Fermi-Dirac Integral *)
FD[j_,η_]:=(1/(Gamma[j+1]))(∫_0^∞)((x^j)/(1+Exp[x-η]))\[DifferentialD]x;
Print["FD[(1/2),η] =",FD[(1/2),η]]  (* Fermi-Dirac integral for j=(1/2) *)
Print["FD[0,η] =",FD[0,η]]  (* Fermi-Dirac integral for j=0 *)
Print["FD[-(1/2),η] =",FD[-(1/2),η]]  (* Fermi-Dirac integral for j=-(1/2) *)
Print["FD[1,η] =",FD[1,η]]  (* Fermi-Dirac integral for j=1 *)
(* Plot of The Fermi-Dirac Integrals above *)
LogPlot[{-PolyLog[(3/2),-(E^η)],Log[1+(E^η)],-PolyLog[(1/2),-(E^η)],-PolyLog[2,-(E^η)]},{η,-8,8},Frame->True, AspectRatio->GoldenRatio,FrameLabel->{"η","Fermi-Dirac Integrals"},BaseStyle->{FontSize->15},PlotStyle->Thick,GridLines->Automatic]

(* 1D Electron Gas: Quantum Currents *)
ClearAll
q:=1.6*(10^(-19)); (* Coulomb, electron charge *)
me:=9.1*(10^(-31)); (* kg, electron rest mass *)
kb:=1.38*(10^(-23)); (* J.K, Boltzmann constant *)
hbar:=(6.63/(2π))*(10^(-34)); (* J.s, Reduced Planck's constant *)
FullSimplify[qC*((gs*gv)/(((2π))^1))*((∫_0^∞)((((h1/m)*k))/(1+Exp[(((((h1^2)*(k^2))/(2*m))-EFL)/(kB*T))]))\[DifferentialD]k) ] (* Evaluating the right-going 1D quantum current *)
J1dR[gs_,gv_,T_,EFL_]:=q* ((gs*gv*kb*T*Log[1+(E^((q*EFL)/(kb*T)))])/(2*π*hbar)); (* Right-going 1D Quantum Current, A *)
LogPlot[{J1dR[2,1,10,EFL],J1dR[2,1,77,EFL],J1dR[2,1,300,EFL]},{EFL,-1,1},Frame->True, AspectRatio->GoldenRatio,FrameLabel->{"(Subscript[E, FS])(eV)","(Subscript[J, 1  d]^R)(A)"},BaseStyle->{FontSize->15},PlotRange->{(10^(-19)),(10^(-3))},PlotStyle->{Blue,Green,Red},GridLines->Automatic](* Evaluating the right-going 1D quantum current *)
FullSimplify[((gs*gv)/(((2π))^1))*((∫_0^∞)(1/(1+Exp[(((((h1^2)*(k^2))/(2*m))-EFL)/(kB*T))]))\[DifferentialD]k) ] (* Evaluating the right-going 1D carrier density *)
n1dR[gs_,gv_,T_,EFL_]:=(10^(-2))* ((gs*gv *(-PolyLog[(1/2),-(E^((q*EFL)/(kb*T)))]))/(2*(Sqrt[2*π])*(Sqrt[((hbar^2)/(kb*me*T))]))) ;(* Right-going 1D Carrier Density in 1/cm *)
LogPlot[{n1dR[2,1,10,EFL],n1dR[2,1,77,EFL],n1dR[2,1,300,EFL]},{EFL,-1,1},Frame->True, AspectRatio->GoldenRatio,FrameLabel->{"(Subscript[E, FS])(eV)","(Subscript[n, 1  d]^R)(1/cm)"},BaseStyle->{FontSize->15},PlotRange->{(10^(-8)),(10^8)},PlotStyle->{Blue,Green,Red},GridLines->Automatic](* Evaluating the right-going 2D quantum current *)
v1dR[T_,EFL_]:=(10^2)*(((gs*gv*kb*T*Log[1+(E^((q*EFL)/(kb*T)))])/(2*π*hbar))/( ((gs*gv *(-PolyLog[(1/2),-(E^((q*EFL)/(kb*T)))]))/(2*(Sqrt[2*π])*(Sqrt[((hbar^2)/(kb*me*T))]))))) ;(* Right-going 1D Ensemble Velocity, cm/s *)
Quiet[LogPlot[{v1dR[4,EFL],v1dR[77,EFL],v1dR[300,EFL]},{EFL,-1,1},Frame->True, AspectRatio->GoldenRatio,FrameLabel->{"(Subscript[E, FS])(eV)","(Subscript[v, 1  d]^R)(cm/s)"},BaseStyle->{FontSize->15},PlotRange->{(10^5),(10^8)},PlotStyle->{Blue,Green,Red},GridLines->Automatic]](* Evaluating the Ensemble Average Right-Going Velocity *)

(* 2D Electron Gas: Quantum Currents *)
ClearAll
q:=1.6*(10^(-19)); (* Coulomb, electron charge *)
me:=9.1*(10^(-31)); (* kg, electron rest mass *)
kb:=1.38*(10^(-23)); (* J.K, Boltzmann constant *)
hbar:=(6.63/(2π))*(10^(-34)); (* J.s, Reduced Planck's constant *)
FullSimplify[qC*((gs*gv)/(((2π))^2))*(∫_0^∞)(∫_(-(π/2))^(π/2))((((h1/m)*k*Cos[θ])*k)/(1+Exp[(((((h1^2)*(k^2))/(2*m))-EFL)/(kB*T))]))\[DifferentialD]θ\[DifferentialD]k ] (* Evaluating the right-going 2D quantum current *)
J2dR[gs_,gv_,T_,EFL_]:=(10^(-3))* q*((gs*gv *hbar)/(2 (Sqrt[2])*me*(π^(3/2))*(((((hbar^2)/(kb*me*T))))^(3/2))))*(-PolyLog[(3/2),-(E^((q*EFL)/(kb*T)))]); (* Right-going 2D Quantum Current, mA/um *)
LogPlot[{J2dR[2,1,10,EFL],J2dR[2,1,77,EFL],J2dR[2,1,300,EFL]},{EFL,-1,1},Frame->True, AspectRatio->GoldenRatio,FrameLabel->{"(Subscript[E, FS])(eV)","(Subscript[J, 2  d]^R)(mA/μm)"},BaseStyle->{FontSize->15},PlotRange->{(10^(-13)),(10^3)},PlotStyle->{Blue,Green,Red},GridLines->Automatic](* Evaluating the right-going 2D quantum current *)
FullSimplify[((gs*gv)/(((2π))^2))*(∫_0^∞)(∫_(-(π/2))^(π/2))(k/(1+Exp[(((((h1^2)*(k^2))/(2*m))-EFL)/(kB*T))]))\[DifferentialD]θ\[DifferentialD]k ] (* Evaluating the right-going 2D carrier density *)
n2dR[gs_,gv_,T_,EFL_]:=(10^(-4))*((gs*gv*kb*me*T*Log[1+(E^((q*EFL)/(kb*T)))])/(4*π*(hbar^2))) ;(* Right-going 2D carrier density, 1/(cm^2) *)
LogPlot[{n2dR[2,1,10,EFL],n2dR[2,1,77,EFL],n2dR[2,1,300,EFL]},{EFL,-1,1},Frame->True, AspectRatio->GoldenRatio,FrameLabel->{"(Subscript[E, FS])(eV)","(Subscript[n, 2  d]^R)(1/(cm^2))"},BaseStyle->{FontSize->15},PlotRange->{(10^0),(10^16)},PlotStyle->{Blue,Green,Red},GridLines->Automatic](* Evaluating the right-going 2D quantum current *)
v2dR[T_,EFL_]:=(10^2)*((((gs*gv *hbar)/(2 (Sqrt[2])*me*(π^(3/2))*(((((hbar^2)/(kb*me*T))))^(3/2))))*(-PolyLog[(3/2),-(E^((q*EFL)/(kb*T)))]))/( ((gs*gv*kb*me*T*Log[1+(E^((q*EFL)/(kb*T)))])/(4*π*(hbar^2))))) ;(* Right-going 2D Ensemble Velocity, cm/s *)
Quiet[LogPlot[{v2dR[4,EFL],v2dR[77,EFL],v2dR[300,EFL]},{EFL,-1,1},Frame->True, AspectRatio->GoldenRatio,FrameLabel->{"(Subscript[E, FS])(eV)","(Subscript[v, 2  d]^R)(cm/s)"},BaseStyle->{FontSize->15},PlotRange->{(10^5),(10^8)},PlotStyle->{Blue,Green,Red},GridLines->Automatic]](* Evaluating the Ensemble Average Right-Going Velocity *)

(* 3D Electron Gas: Quantum Currents *)
ClearAll
q:=1.6*(10^(-19)); (* Coulomb, electron charge *)
me:=9.1*(10^(-31)); (* kg, electron rest mass *)
kb:=1.38*(10^(-23)); (* J.K, Boltzmann constant *)
hbar:=(6.63/(2π))*(10^(-34)); (* J.s, Reduced Planck's constant *)
FullSimplify[qC*((gs*gv)/(((2π))^2))*(∫_0^∞)(∫_0^(2π))(∫_0^(π/2))((((h1/m)*k*Cos[θ])*(k^2)*Sin[θ])/(1+Exp[(((((h1^2)*(k^2))/(2*m))-EFL)/(kB*T))]))\[DifferentialD]θ\[DifferentialD]ϕ\[DifferentialD]k] (* Evaluating the right-going 2D quantum current *)
J3dR[gs_,gv_,T_,EFL_]:= (10^(-4))*q*((gs*gv*me*(kb^2)*(T^2) *(-PolyLog[2,-(E^((q*EFL)/(kb*T)))]))/(4*(π^2) *(hbar^3))); (* Right-going 3D Quantum Current, A/(cm^2) *)
LogPlot[{J3dR[2,1,4,EFL],J3dR[2,1,77,EFL],J3dR[2,1,300,EFL]},{EFL,-1,1},Frame->True, AspectRatio->GoldenRatio,FrameLabel->{"(Subscript[E, FS])(eV)","(Subscript[J, 3  d]^R)(A/(cm^2))"},BaseStyle->{FontSize->15},PlotRange->{(10^(-4)),(10^12)},PlotStyle->{Blue,Green,Red},GridLines->Automatic](* Evaluating the right-going 2D quantum current *)
FullSimplify[((gs*gv)/(((2π))^3))*(∫_0^∞)(∫_0^(2π))(∫_0^(π/2))(((k^2)*Sin[θ])/(1+Exp[(((((h1^2)*(k^2))/(2*m))-EFL)/(kB*T))]))\[DifferentialD]θ\[DifferentialD]ϕ\[DifferentialD]k ] (* Evaluating the right-going 3D carrier density *)
n3dR[gs_,gv_,T_,EFL_]:=(10^(-6))*((gs*gv*(-PolyLog[(3/2),-(E^((q*EFL)/(kb*T)))]))/(4*(Sqrt[2])*(π^(3/2))*(((((hbar^2)/(kb*me*T))))^(3/2)))) ;(* Right-going 3D carrier density, 1/(cm^3) *)
LogPlot[{n3dR[2,1,4,EFL],n3dR[2,1,77,EFL],n3dR[2,1,300,EFL]},{EFL,-1,1},Frame->True, AspectRatio->GoldenRatio,FrameLabel->{"(Subscript[E, FS])(eV)","(Subscript[n, 3  d]^R)(1/(cm^3))"},BaseStyle->{FontSize->15},PlotRange->{(10^8),(10^24)},PlotStyle->{Blue,Green,Red},GridLines->Automatic](* Evaluating the right-going 2D quantum current *)
v3dR[T_,EFL_]:=(10^2)*(((gs*gv*me*(kb^2)*(T^2) *(-PolyLog[2,-(E^((q*EFL)/(kb*T)))]))/(4*(π^2) *(hbar^3)))/((gs*gv*(-PolyLog[(3/2),-(E^((q*EFL)/(kb*T)))]))/(4*(Sqrt[2])*(π^(3/2))*(((((hbar^2)/(kb*me*T))))^(3/2))))) ;(* Right-going 3D Ensemble Velocity, cm/s *)
Quiet[LogPlot[{v3dR[4,EFL],v3dR[77,EFL],v3dR[300,EFL]},{EFL,-1,1},Frame->True, AspectRatio->GoldenRatio,FrameLabel->{"(Subscript[E, FS])(eV)","(Subscript[v, 3  d]^R)(cm/s)"},BaseStyle->{FontSize->15},PlotRange->{(10^5),(10^8)},PlotStyle->{Blue,Green,Red},GridLines->Automatic]](* Evaluating the Ensemble Average Right-Going Velocity *)

(* Ensemble velocities for 1D, 2D, and 3D Electron Systems *)
Quiet[LogPlot[{v1dR[4,EFL],v1dR[77,EFL],v1dR[300,EFL],v2dR[4,EFL],v2dR[77,EFL],v2dR[300,EFL],v3dR[4,EFL],v3dR[77,EFL],v3dR[300,EFL]},{EFL,-1,1},Frame->True, AspectRatio->GoldenRatio,FrameLabel->{"(Subscript[E, FS])(eV)","(v^R)(cm/s)"},BaseStyle->{FontSize->15},PlotRange->{(10^5),(10^8)},PlotStyle->Thick,GridLines->Automatic]]
Quiet[Plot[{v1dR[4,EFL],v1dR[77,EFL],v1dR[300,EFL],v2dR[4,EFL],v2dR[77,EFL],v2dR[300,EFL],v3dR[4,EFL],v3dR[77,EFL],v3dR[300,EFL]},{EFL,-1,1},Frame->True, AspectRatio->GoldenRatio,FrameLabel->{"(Subscript[E, FS])(eV)","(v^R)(cm/s)"},BaseStyle->{FontSize->15},PlotRange->{0,4*(10^7)},GridLines->Automatic]]

(* Non-Degenerate Limit of Ensemble Velocities for (EFL/(kb*T))<<-1 *)
ClearAll
Print["v1dND=",FullSimplify[(((gs*gv*kB*T)/(2*π*h))/( ((gs*gv )/(2*(Sqrt[2*π])*(Sqrt[((h^2)/(kB*m*T))])))))]]
Print["v2dND=",FullSimplify[(((gs*gv *h)/(2 (Sqrt[2])*m*(π^(3/2))*(((((h^2)/(kB*m*T))))^(3/2))))/( ((gs*gv*kB*m*T)/(4*π*(h^2)))))]]
Print["v3dND=",FullSimplify[(((gs*gv*m*(kB^2)*(T^2))/(4*(π^2) *(h^3)))/((gs*gv)/(4*(Sqrt[2])*(π^(3/2))*(((((h^2)/(kB*m*T))))^(3/2)))))]]
N[(Sqrt[(2/π)])]

v1dPlot=Quiet[Plot[{v1dR[4,EFL],v1dR[77,EFL],v1dR[300,EFL]},{EFL,-0.2,1},Frame->True, AspectRatio->1,FrameLabel->{"(Subscript[E, FS])(eV)","(v^R)(cm/s)"},BaseStyle->{FontSize->15},PlotStyle->{Blue,Green,Red},PlotRange->{0,4*(10^7)},GridLines->Automatic]];
v2dPlot=Quiet[Plot[{v2dR[4,EFL],v2dR[77,EFL],v2dR[300,EFL]},{EFL,-0.2,1},Frame->True, AspectRatio->1,FrameLabel->{"(Subscript[E, FS])(eV)","(v^R)(cm/s)"},BaseStyle->{FontSize->15},PlotStyle->{Blue,Green,Red},PlotRange->{0,4*(10^7)},GridLines->Automatic]];
v3dPlot=Quiet[Plot[{v3dR[4,EFL],v3dR[77,EFL],v3dR[300,EFL]},{EFL,-1,1},Frame->True, AspectRatio->1,FrameLabel->{"(Subscript[E, FS])(eV)","(v^R)(cm/s)"},BaseStyle->{FontSize->15},PlotStyle->{Blue,Green,Red},PlotRange->{0,4*(10^7)},GridLines->Automatic]];
Show[v1dPlot,v2dPlot,v3dPlot]

(* Fermi Energies and Specific Heat for 1D, 2D, and 3D free Electrons *)
Series[x*Log[Exp[(1/x)]-1]+((π^2)/12)*x*(1/(Log[Exp[(1/x)]-1])),{x,0,5}]
Series[(Sqrt[1-((2*(π^2))/3)*(x^2)]),{x,0,5}]
Series[Log[1-x],{x,0,5}]
Series[Log[1+x],{x,0,5}]
LogPlot[{x*Log[Exp[(1/x)]-1],1-((π^2)/(10*6))*(x^2)},{x,0,0.3},Frame->True]
FullSimplify[Solve[(x^2)-x*ef0+((π^2)/6)*(((kb*T))^2)==0,x]]
LogLogPlot[{Log[1+Exp[x]],x+((π^2)/6)*(1/x)},{x,1,1000},Frame->True]
f[x_]:=x*Log[Exp[(1/x)]-1]+((π^2)/12)*(x/(Log[Exp[(1/x)]-1]));
FullSimplify[D[x*Log[Exp[(1/x)]-1]+((π^2)/12)*(x/(Log[Exp[(1/x)]-1])),x]]
FullSimplify[Series[D[x*Log[Exp[(1/x)]-1]+((π^2)/12)*(x/(Log[Exp[(1/x)]-1])),x],{x,0,5}]]
Plot[{(1/(-x+(E^(-1/x)) x))+(((E^(1/x)) (π^2))/(12 (-1+(E^(1/x))) x ((Log[-1+(E^(1/x))])^2)))+((π^2)/(12 Log[-1+(E^(1/x))]))+Log[-1+(E^(1/x))]},{x,0,1}]
Plot[{1-((π^2)/12)*(x^2),x*Log[Exp[(1/x)]-1],1+((π^2)/12)*(x^(5/2))},{x,0,1}]

Quantum Currents for the Free Electron Gas

ClearAll
q:=1.6*(10^(-19)); (* Electron charge in Coulomb *)
me:=9.1*(10^(-31)); (* Free electron rest mass in kg *)
h:=(6.63/(2*π))*(10^(-34));  (* Reduced Planck's constant in J.s *)
kb:=1.38*(10^(-23)); (* Boltzmann constant in J/K *)
FD[j_,η_]:=(1/(Gamma[j+1]))(∫_0^∞)((x^j)/(1+Exp[x-η]))\[DifferentialD]x;  (* Defining the Fermi-Dirac integrals *)
Nc[d_,gs_,gv_,meff_,T_]:=((gs*gv)/(((2*π))^d))*(((((2*π*meff*me*kb*T)/(h^2))))^(d/2));  (* Effective band-edge DOS in d dimensions *)
n[d_,gs_,gv_,meff_,T_,ηF_]:=Nc[d,gs,gv,meff,T]*FD[((d-2)/2),ηF]; (* Effective SI carrier density in d dimensions *)
n1d[gs_,gv_,meff_,T_,ηF_]:=Nc[1,gs,gv,meff,T]*(-PolyLog[(1/2),-(E^(ηF))]); (* Effective SI carrier density in 1D *)
n2d[gs_,gv_,meff_,T_,ηF_]:=Nc[2,gs,gv,meff,T]*(Log[1+(E^(ηF))]); (* Effective SI carrier density in 2D *)
n3d[gs_,gv_,meff_,T_,ηF_]:=Nc[3,gs,gv,meff,T]*(-PolyLog[(3/2),-(E^(ηF))]); (* Effective SI carrier density in 3D *)
ηS[d_,gs_,gv_,meff_,T_,v_,nd_]:=Quiet[Chop[FindRoot[(1/2)*(n[d,gs,gv,meff,T,η]+n[d,gs,gv,meff,T,η-v])==nd,{η,1}]][[1]][[2]]]; (* Source Fermi Level at voltage v in d dimensions *)
ηS1d[gs_,gv_,meff_,T_,v_,nd_]:=Quiet[Chop[FindRoot[(1/2)*(n1d[gs,gv,meff,T,η]+n1d[gs,gv,meff,T,η-v])==nd,{η,1}]][[1]][[2]]]; (* Source Fermi Level at voltage v in 1D *)
ηS2d[gs_,gv_,meff_,T_,v_,nd_]:=Quiet[Chop[FindRoot[(1/2)*(n2d[gs,gv,meff,T,η]+n2d[gs,gv,meff,T,η-v])==nd,{η,1}]][[1]][[2]]]; (* Source Fermi Level at voltage v in 2D *)
ηS3d[gs_,gv_,meff_,T_,v_,nd_]:=Quiet[Chop[FindRoot[(1/2)*(n3d[gs,gv,meff,T,η]+n3d[gs,gv,meff,T,η-v])==nd,{η,1}]][[1]][[2]]];(* Source Fermi Level at voltage v in 3D *)
(* Dimensionless Source Fermi level in in the Ballistic Limit at a drain voltage v *)
J1dQuantum[gs_,gv_,meff_,T_,v_,nd_]:=Chop[(q/(2*π*h))*(((Nc[0,gs,gv,meff,T])/1))*(kb*T)*Log[((1+Exp[ηS1d[gs,gv,meff,T,v,nd]])/(1+Exp[ηS1d[gs,gv,meff,T,v,nd]-v]))]]; (* Quantum Current at voltage v in 1D *)
J2dQuantum[gs_,gv_,meff_,T_,v_,nd_]:=Chop[(q/(2*π*h))*(((Nc[1,gs,gv,meff,T])/2))*(kb*T)*(-PolyLog[(3/2),-(E^(ηS2d[gs,gv,meff,T,v,nd]))]-(-PolyLog[(3/2),-(E^((ηS2d[gs,gv,meff,T,v,nd]-v)))]))]; (* Quantum Current at voltage v in 2D *)
J3dQuantum[gs_,gv_,meff_,T_,v_,nd_]:=Chop[(q/(2*π*h))*((((2*π))^3)*((Nc[2,gs,gv,meff,T])/1))*(kb*T)*(-PolyLog[2,-(E^(ηS3d[gs,gv,meff,T,v,nd]))]-(-PolyLog[2,-(E^((ηS3d[gs,gv,meff,T,v,nd]-v)))]))];(* Quantum Current at voltage v in 3D *)

1D

(* ηS1d[gs_,gv_,meff_,T_,v_,nd_] *)
ηS1d[2,1,1,300,0.0,(10^4)*(10^2)]
ηS1d[2,1,1,300,0.0,(10^8)*(10^2)]
LogLinearPlot[{ηS1d[2,1,1,300,0.0,(10^4)*(10^2)],ηS1d[2,1,1,300,0.0,n1d*(10^2)]},{n1d,(10^4),(10^8)}]

(* ηS1d[gs_,gv_,meff_,T_,v_,nd_] *)LogLinearPlot[{ηS1d[2,1,1,4,0.0,n1d*(10^2)],ηS1d[2,1,1,77,0.0,n1d*(10^2)],ηS1d[2,1,1,300,0.0,n1d*(10^2)]},{n1d,(10^4),(10^8)},Frame->True,PlotStyle->{{Blue,Thick},{Green,Thick},{Red,Thick}},FrameLabel->{"(Subscript[n, 1  d])(1/cm)","(Subscript[η, F])=(Subscript[E, F])/(Subscript[k, B])T"},PlotRange->{-20,100},BaseStyle->{FontSize->15},AspectRatio->GoldenRatio,GridLines->Automatic]
Plot[{ηS1d[2,1,1,300,0,(10^10)],ηS1d[2,1,1,300,v,(10^10)],ηS1d[2,1,1,300,v,(10^10)]-v},{v,0,40},Frame->True,PlotStyle->{{Black,Dashed},{Black,Thick},{Black,Thick}},FrameLabel->{"(Subscript[v, d])((Subscript[k, B])T/q)","η"},PlotRange->Automatic,BaseStyle->{FontSize->15},AspectRatio->GoldenRatio,GridLines->Automatic]
Plot[{ηS1d[2,1,1,300,0,(10^7)],ηS1d[2,1,1,300,v,(10^7)],ηS1d[2,1,1,300,v,(10^7)]-v},{v,0,40},Frame->True,PlotStyle->{{Black,Dashed},{Black,Thick},{Black,Thick}},FrameLabel->{"(Subscript[v, d])((Subscript[k, B])T/q)","η"},PlotRange->Automatic,BaseStyle->{FontSize->15},AspectRatio->GoldenRatio,GridLines->Automatic]

(* 1D Fermi Level Split with Drain Voltage *)
(* ηS1d[gs_,gv_,meff_,T_,v_,nd_] *)
Plot[{ηS1d[2,1,1,300,0,(10^8)],ηS1d[2,1,1,300,v,(10^8)],ηS1d[2,1,1,300,v,(10^8)]-v,ηS1d[2,1,1,300,0,5*(10^9)],ηS1d[2,1,1,300,v,5*(10^9)],ηS1d[2,1,1,300,v,5*(10^9)]-v},{v,0,500},Frame->True,PlotStyle->{{Blue,Dashed},{Blue,Thick},{Blue,Thick},{Red,Dashed},{Red,Thick},{Red,Thick}},FrameLabel->{"(Subscript[v, d])((Subscript[k, B])T/q)","η"},PlotRange->{-200,400},BaseStyle->{FontSize->15},AspectRatio->GoldenRatio,GridLines->Automatic]

(* 1D Quantum Current *)
(*J1dQuantum[gs_,gv_,meff_,T_,v_,nd_]:=Chop[(q/(2*π*h))*(((Nc[0,gs,gv,meff,T])/1))*(kb*T)*Log[((1+Exp[ηS1d[gs,gv,meff,T,v,nd]])/(1+Exp[ηS1d[gs,gv,meff,T,v,nd]-v]))]];*)
Plot[{(10^6)*J1dQuantum[2,1,1,300,v,1*(10^8)],(10^6)*J1dQuantum[2,1,1,300,v,(10^9)],(10^6)*J1dQuantum[2,1,1,300,v,2*(10^9)],(10^6)*J1dQuantum[2,1,1,300,v,3*(10^9)],(10^6)*J1dQuantum[2,1,1,300,v,4*(10^9)],(10^6)*J1dQuantum[2,1,1,300,v,5*(10^9)]},{v,0,500},Frame->True,PlotStyle->{{Thick,Blue},{Thick,Black},{Thick,Green},{Thick,Orange},{Thick,Cyan},{Thick,Red}},FrameLabel->{"(Subscript[v, d])((Subscript[k, B])T/q)","(Subscript[J, 1  d])(μA)"},PlotRange->{0,800},BaseStyle->{FontSize->15},AspectRatio->GoldenRatio,GridLines->Automatic]
Plot[{(10^6)*J1dQuantum[2,1,1,77,v,1*(10^8)],(10^6)*J1dQuantum[2,1,1,77,v,(10^9)],(10^6)*J1dQuantum[2,1,1,77,v,2*(10^9)],(10^6)*J1dQuantum[2,1,1,77,v,3*(10^9)],(10^6)*J1dQuantum[2,1,1,77,v,4*(10^9)],(10^6)*J1dQuantum[2,1,1,77,v,5*(10^9)]},{v,0,500},Frame->True,PlotStyle->{{Thick,Blue},{Thick,Black},{Thick,Green},{Thick,Orange},{Thick,Cyan},{Thick,Red}},FrameLabel->{"(Subscript[v, d])((Subscript[k, B])T/q)","(Subscript[J, 1  d])(μA)"},PlotRange->{0,800},BaseStyle->{FontSize->15},AspectRatio->GoldenRatio,GridLines->Automatic]

J1dQuantum[2,1,1,300,500,(10^4)*(10^2)]
J1dQuantum[2,1,1,300,500,(10^8)*(10^2)]
LogLogPlot[{J1dQuantum[2,1,1,300,1,n1d]},{n1d,(10^6),(10^10)}]

(*J1dQuantum[gs_,gv_,meff_,T_,v_,nd_]:=((q*gs*gv)/(2*π*h))*(kb*T)*Log[((1+Exp[ηS1d[gs,gv,meff,T,v,nd]])/(1+Exp[ηS1d[gs,gv,meff,T,v,nd]-v]))];*)
LogLogPlot[{J1dQuantum[2,1,1,300,0.01,n1d*(10^2)],J1dQuantum[2,1,1,300,0.1,n1d*(10^2)],J1dQuantum[2,1,1,300,1,n1d*(10^2)]},{n1d,(10^4),(10^8)},Frame->True,PlotStyle->Thick,FrameLabel->{"(Subscript[n, 1  d])(/cm)","(Subscript[J, 1  d])(μA/μm)"},PlotRange->All,BaseStyle->{FontSize->15},AspectRatio->GoldenRatio,GridLines->Automatic]

2D

(* 2D Fermi Level Split with Drain Voltage *)
(* ηS2d[gs_,gv_,meff_,T_,v_,nd_] *)
Plot[{ηS2d[2,1,1,300,0,(10^17)],ηS2d[2,1,1,300,v,(10^17)],ηS2d[2,1,1,300,v,(10^17)]-v,ηS2d[2,1,1,300,0,(10^19)],ηS2d[2,1,1,300,v,(10^19)],ηS2d[2,1,1,300,v,(10^19)]-v},{v,0,500},Frame->True,PlotStyle->{{Blue,Dashed},{Blue,Thick},{Blue,Thick},{Red,Dashed},{Red,Thick},{Red,Thick}},FrameLabel->{"(Subscript[v, d])((Subscript[k, B])T/q)","η"},PlotRange->{-200,400},BaseStyle->{FontSize->15},AspectRatio->GoldenRatio,GridLines->Automatic]

(*J2dQuantum[gs_,gv_,meff_,T_,v_,nd_]:=(q/(2*π*h))*(((Nc[1,gs,gv,meff,T])/2))*(kb*T)*(-PolyLog[(3/2),-(E^(ηS2d[gs,gv,meff,T,v,nd]))]-(-PolyLog[(3/2),-(E^((ηS2d[gs,gv,meff,T,v,nd]-v)))]));*)
LogLogPlot[{J2dQuantum[2,1,1,300,0.01,nd*(10^16)],J2dQuantum[2,1,1,300,0.1,nd*(10^16)],J2dQuantum[2,1,1,300,1,nd*(10^16)]},{nd,(10^(-4)),(10^4)},Frame->True,PlotStyle->Thick,FrameLabel->{"(Subscript[n, 2  d])((10^12)/(cm^2))","(Subscript[J, 2  d])(μA/μm)"},PlotRange->All,BaseStyle->{FontSize->15},AspectRatio->GoldenRatio,GridLines->Automatic]

(*J2dQuantum[gs_,gv_,meff_,T_,v_,nd_]:=(q/(2*π*h))*(((Nc[1,gs,gv,meff,T])/2))*(kb*T)*(-PolyLog[(3/2),-(E^(ηS2d[gs,gv,meff,T,v,nd]))]-(-PolyLog[(3/2),-(E^((ηS2d[gs,gv,meff,T,v,nd]-v)))]));*)
J2dQuantum[2,1,1,300,100,(10^16)]
Plot[{J2dQuantum[2,1,1,300,v,1*(10^14)],J2dQuantum[2,1,1,300,v,(10^19)],J2dQuantum[2,1,1,300,v,2*(10^19)],J2dQuantum[2,1,1,300,v,3*(10^19)],J2dQuantum[2,1,1,300,v,4*(10^19)],J2dQuantum[2,1,1,300,v,5*(10^19)]},{v,0,1000},Frame->True,PlotStyle->{{Thick,Blue},{Thick,Black},{Thick,Green},{Thick,Orange},{Thick,Cyan},{Thick,Red}},FrameLabel->{"(Subscript[v, d])((Subscript[k, B])T/q)","(Subscript[J, 2  d])(μA/μm)"},PlotRange->{0,5*(10^6)},BaseStyle->{FontSize->15},AspectRatio->GoldenRatio,GridLines->Automatic]
Plot[{J2dQuantum[2,1,1,77,v,1*(10^14)],J2dQuantum[2,1,1,77,v,(10^19)],J2dQuantum[2,1,1,77,v,2*(10^19)],J2dQuantum[2,1,1,77,v,3*(10^19)],J2dQuantum[2,1,1,77,v,4*(10^19)],J2dQuantum[2,1,1,77,v,5*(10^19)]},{v,0,1000},Frame->True,PlotStyle->{{Thick,Blue},{Thick,Black},{Thick,Green},{Thick,Orange},{Thick,Cyan},{Thick,Red}},FrameLabel->{"(Subscript[v, d])((Subscript[k, B])T/q)","(Subscript[J, 2  d])(μA/μm)"},PlotRange->{0,5*(10^6)},BaseStyle->{FontSize->15},AspectRatio->GoldenRatio,GridLines->Automatic]

(*J2dQuantum[gs_,gv_,meff_,T_,v_,nd_]:=(q/(2*π*h))*(((Nc[1,gs,gv,meff,T])/2))*(kb*T)*(-PolyLog[(3/2),-(E^(ηS2d[gs,gv,meff,T,v,nd]))]-(-PolyLog[(3/2),-(E^((ηS2d[gs,gv,meff,T,v,nd]-v)))]));*)
J2dQuantum[2,1,1,300,100,(10^16)]
Plot[{J2dQuantum[2,1,1,300,v/(((kb*300)/q)),1*(10^14)],J2dQuantum[2,1,1,300,v/(((kb*300)/q)),(10^19)],J2dQuantum[2,1,1,300,v/(((kb*300)/q)),2*(10^19)],J2dQuantum[2,1,1,300,v/(((kb*300)/q)),3*(10^19)],J2dQuantum[2,1,1,300,v/(((kb*300)/q)),4*(10^19)],J2dQuantum[2,1,1,300,v/(((kb*300)/q)),5*(10^19)]},{v,0,26},Frame->True,PlotStyle->{{Thick,Blue},{Thick,Black},{Thick,Green},{Thick,Orange},{Thick,Cyan},{Thick,Red}},FrameLabel->{"(Subscript[V, d])(Volt)","(Subscript[J, 2  d])(μA/μm)"},PlotRange->{-1*(10^5),5*(10^6)},BaseStyle->{FontSize->15},AspectRatio->GoldenRatio,GridLines->Automatic]
Plot[{J2dQuantum[2,1,1,77,v/(((kb*77)/q)),1*(10^14)],J2dQuantum[2,1,1,77,v/(((kb*77)/q)),(10^19)],J2dQuantum[2,1,1,77,v/(((kb*77)/q)),2*(10^19)],J2dQuantum[2,1,1,77,v/(((kb*77)/q)),3*(10^19)],J2dQuantum[2,1,1,77,v/(((kb*77)/q)),4*(10^19)],J2dQuantum[2,1,1,77,v/(((kb*77)/q)),5*(10^19)]},{v,0,26},Frame->True,PlotStyle->{{Thick,Blue},{Thick,Black},{Thick,Green},{Thick,Orange},{Thick,Cyan},{Thick,Red}},FrameLabel->{"(Subscript[V, d])(Volt)","(Subscript[J, 2  d])(μA/μm)"},PlotRange->{-1*(10^5),5*(10^6)},BaseStyle->{FontSize->15},AspectRatio->GoldenRatio,GridLines->Automatic]

J2dQuantum[2,1,1,300,0.2,(10^16)]
J2dQuantum[2,1,1,300,1000,(10^16)]
J2dQuantum[2,1,1,300,1000,(10^19)]
LogLogPlot[J2dQuantum[2,1,1,300,200,n2d],{n2d,(10^15),(10^19)},Frame->True,PlotStyle->Thick,FrameLabel->{"(Subscript[n, 1  d])(/(cm^2))","(Subscript[J, 1  d])(A)"},PlotRange->All,BaseStyle->{FontSize->15},AspectRatio->GoldenRatio,GridLines->Automatic]

(* ηS[d_,gs_,gv_,meff_,T_,v_,nd_] *)
LogLinearPlot[{ηS2d[2,2,1,1,4,0.0,n2d*(10^4)],ηS[2,2,1,1,77,0.0,n2d*(10^4)],ηS[2,2,1,1,300,0.0,n2d*(10^4)]},{n2d,(10^10),(10^16)},Frame->True,PlotStyle->{{Blue,Thick},{Green,Thick},{Red,Thick}},FrameLabel->{"(Subscript[n, 2  d])(1/(cm^2))","(Subscript[η, F])=(Subscript[E, F])/(Subscript[k, B])T"},PlotRange->{-20,100},BaseStyle->{FontSize->15},AspectRatio->GoldenRatio,GridLines->Automatic]
Plot[{ηS[2,2,1,1,300,0,(10^19)],ηS[2,2,1,1,300,v,(10^19)],ηS[2,2,1,1,300,v,(10^19)]-v},{v,0,40},Frame->True,PlotStyle->{{Black,Dashed},{Black,Thick},{Black,Thick}},FrameLabel->{"v((Subscript[k, B])T/q)","η"},PlotRange->Automatic,BaseStyle->{FontSize->15},AspectRatio->GoldenRatio,GridLines->Automatic]
Plot[{ηS[2,2,1,1,300,0,(10^16)],ηS[2,2,1,1,300,v,(10^16)],ηS[2,2,1,1,300,v,(10^16)]-v},{v,0,40},Frame->True,PlotStyle->{{Black,Dashed},{Black,Thick},{Black,Thick}},FrameLabel->{"v((Subscript[k, B])T/q)","η"},PlotRange->Automatic,BaseStyle->{FontSize->15},AspectRatio->GoldenRatio,GridLines->Automatic]

3D

(* 3D Fermi Level Split with Drain Voltage *)
(* ηS3d[gs_,gv_,meff_,T_,v_,nd_] *)
Plot[{ηS3d[2,1,1,300,0,(10^21)],ηS3d[2,1,1,300,v,(10^21)],ηS3d[2,1,1,300,v,(10^21)]-v,ηS3d[2,1,1,300,0,(10^29)],ηS3d[2,1,1,300,v,(10^29)],ηS3d[2,1,1,300,v,(10^29)]-v},{v,0,1000},Frame->True,PlotStyle->{{Blue,Dashed},{Blue,Thick},{Blue,Thick},{Red,Dashed},{Red,Thick},{Red,Thick}},FrameLabel->{"(Subscript[v, d])((Subscript[k, B])T/q)","η"},PlotRange->All,BaseStyle->{FontSize->15},AspectRatio->GoldenRatio,GridLines->Automatic]

(*J3dQuantum[gs_,gv_,meff_,T_,v_,nd_]:=Chop[(q/(2*π*h))*((((2*π))^3)*((Nc[2,gs,gv,meff,T])/1))*(kb*T)*(-PolyLog[2,-(E^(ηS3d[gs,gv,meff,T,v,nd]))]-(-PolyLog[2,-(E^((ηS3d[gs,gv,meff,T,v,nd]-v)))]))];*)
Plot[{J3dQuantum[2,1,1,300,v,1*(10^24)],J3dQuantum[2,1,1,300,v,(10^29)],J3dQuantum[2,1,1,300,v,2*(10^29)],J3dQuantum[2,1,1,300,v,3*(10^29)],J3dQuantum[2,1,1,300,v,4*(10^29)],J3dQuantum[2,1,1,300,v,5*(10^29)]},{v,0,1000},Frame->True,PlotStyle->{{Thick,Blue},{Thick,Black},{Thick,Green},{Thick,Orange},{Thick,Cyan},{Thick,Red}},FrameLabel->{"(Subscript[v, d])((Subscript[k, B])T/q)","(Subscript[J, 3  d])(A/(m^2))"},PlotRange->{0,2.5*(10^19)},BaseStyle->{FontSize->15},AspectRatio->GoldenRatio,GridLines->Automatic]
Plot[{J3dQuantum[2,1,1,77,v,1*(10^24)],J3dQuantum[2,1,1,77,v,(10^29)],J3dQuantum[2,1,1,77,v,2*(10^29)],J3dQuantum[2,1,1,77,v,3*(10^29)],J3dQuantum[2,1,1,77,v,4*(10^29)],J3dQuantum[2,1,1,77,v,5*(10^29)]},{v,0,1000},Frame->True,PlotStyle->{{Thick,Blue},{Thick,Black},{Thick,Green},{Thick,Orange},{Thick,Cyan},{Thick,Red}},FrameLabel->{"(Subscript[v, d])((Subscript[k, B])T/q)","(Subscript[J, 3  d])(A/(m^2))"},PlotRange->{0,2.5*(10^19)},BaseStyle->{FontSize->15},AspectRatio->GoldenRatio,GridLines->Automatic]

Plot[{J3dQuantum[2,1,1,300,v/(((kb*300)/q)),1*(10^24)],J3dQuantum[2,1,1,300,v/(((kb*300)/q)),(10^29)],J3dQuantum[2,1,1,300,v/(((kb*300)/q)),2*(10^29)],J3dQuantum[2,1,1,300,v/(((kb*300)/q)),3*(10^29)],J3dQuantum[2,1,1,300,v/(((kb*300)/q)),4*(10^29)],J3dQuantum[2,1,1,300,v/(((kb*300)/q)),5*(10^29)]},{v,0,26},Frame->True,PlotStyle->{{Thick,Blue},{Thick,Black},{Thick,Green},{Thick,Orange},{Thick,Cyan},{Thick,Red}},FrameLabel->{"(Subscript[V, d])(Volt)","(Subscript[J, 3  d])(A/(m^2))"},PlotRange->All,BaseStyle->{FontSize->15},AspectRatio->GoldenRatio,GridLines->Automatic]
Plot[{J3dQuantum[2,1,1,77,v/(((kb*77)/q)),1*(10^24)],J3dQuantum[2,1,1,77,v/(((kb*77)/q)),(10^29)],J3dQuantum[2,1,1,77,v/(((kb*77)/q)),2*(10^29)],J3dQuantum[2,1,1,77,v/(((kb*77)/q)),3*(10^29)],J3dQuantum[2,1,1,77,v/(((kb*77)/q)),4*(10^29)],J3dQuantum[2,1,1,77,v/(((kb*77)/q)),5*(10^29)]},{v,0,26},Frame->True,PlotStyle->{{Thick,Blue},{Thick,Black},{Thick,Green},{Thick,Orange},{Thick,Cyan},{Thick,Red}},FrameLabel->{"(Subscript[V, d])(Volt)","(Subscript[J, 3  d])(A/(m^2))"},PlotRange->All,BaseStyle->{FontSize->15},AspectRatio->GoldenRatio,GridLines->Automatic]

6. Red or Blue Pill: Befriending the Matrix

M={{1,1,0},{0,2,0},{0,1,1}};  (* See Basics 2 above, time to use matrices! *)
MatrixForm[M]
MatrixForm[Inverse[M]]
MatrixForm[M.Inverse[M]]
Eigenvalues[M]
Eigenvectors[M]
MatrixForm[Eigenvectors[M]]
MatrixForm[Eigensystem[M]]

(* The Eigenvalues and Eigenfunctions of a 2x2 matrix *)
H2={{h11,h12},{h12,h22}};
Print["H2=",MatrixForm[H2]]
Print["EigenSystem[H2]=",MatrixForm[Eigensystem[H2]]]

M={{1,1,0},{0,2,0},{0,1,1}};
Print["M=",MatrixForm[M]]
Print["EigenSystem[M]=",MatrixForm[Eigensystem[M]]]

H5={{E0,-t,0,0,-t},{-t,E0,-t,0,0},{0,-t,E0,-t,0},{0,0,-t,E0,-t},{-t,0,0,-t,E0}};
Print["H5=",MatrixForm[H5]]
Print["EigenSystem[H5]=",MatrixForm[Eigensystem[H5]]]
Print["Trace[H5]=",Tr[H5]]

7. Perturbations to the Electron’s Freedom

The variational method

Variational Method: ⟨E⟩=((⟨ψ|H|ψ⟩)/(⟨ψ|ψ⟩))≥(Subscript[E, 0]) holds for ANY trial wavefunction ψ. Here (Subscript[E, 0]) is the ground state energy

(((∫_(-∞)^∞)Exp[-α*(x^2)]*(-((ℏ^2)/(2*me))*D[Exp[-α*(x^2)],{x,2}]+(1/2)*me*(ω^2)*(x^2)*Exp[-α*(x^2)])\[DifferentialD]x)/((∫_(-∞)^∞)Exp[-α*(x^2)]*Exp[-α*(x^2)]\[DifferentialD]x))

Solve[D[(((me^2) (ω^2)+4 (α^2) (ℏ^2))/(8 me α)),{α,1}]==0,α]

(((me^2) (ω^2)+4 (α^2) (ℏ^2))/(8 me α))/.α->((me ω)/(2 ℏ))

The mechanics of perturbation theory

Simple example of the mechanics of perturbation theory: Solve the equation (x^3)-4.001x+0.002=0

eq1[x_]:=(x^3)-4.001*x+0.002;
Solve[eq1[x]==0,x]  (* This is the 'exact solution' *)

eq2[x_,eps_]:=(x^3)-(4+eps)*x+2*eps;
xp[a0_,a1_,a2_,a3_,eps_]:=a0+a1*eps+a2*(eps^2)+a3*(eps^3);
Collect[eq2[xp[a0,a1,a2,a3,eps],eps],eps]

Clear[a0,a1,a2,a3]
Solve[-4* a0+(a0^3)==0,a0]
a0=Solve[-4* a0+(a0^3)==0,a0][[1,1,2]]
a1=Solve[2-a0-4 a1+3 (a0^2) a1==0,a1][[1,1,2]]
a2=Solve[-a1+3 a0 (a1^2)-4 a2+3 (a0^2) a2==0,a2][[1,1,2]]
a3=Solve[(a1^3)-a2+6 a0 a1 a2-4 a3+3 (a0^2) a3==0,a3][[1,1,2]]
xp[a0,a1,a2,a3,0.001]

8 lines of corrupt data deleted.

(* note how Mathematica selects out numbers from a list expression *)
Clear[a0,a1,a2,a3]
Solve[-4* a0+(a0^3)==0,a0]
Solve[-4* a0+(a0^3)==0,a0][[1]]  (* Use 'esc [[ esc' to get [[ *)
Solve[-4* a0+(a0^3)==0,a0][[1,1]]
Solve[-4* a0+(a0^3)==0,a0][[1,1,2]]
Solve[-4* a0+(a0^3)==0,a0][[2,1,2]]
Solve[-4* a0+(a0^3)==0,a0][[3,1,2]]

The 2-state perturbation problem

Here we investigate the perturbation of a 2-state system to illustrate a few features of perturbation theory.

ClearAll[h2]
h2[a_,b_]:={{a,0},{0,b}};
Eigensystem[h2[a,b]]//MatrixForm  (* The command Eigensystem[...] gives: First row: eigenvalues, Second row: corresponding eigenvectors *)

ClearAll[w]
w[ϵ_]:=ϵ*{{0,1},{1,0}};  (* perturbation to the 2-state system *)
Eigensystem[h2[a,b]+w[ϵ]]//MatrixForm
Eigenvalues[h2[a,b]+w[ϵ]][[1]]
Eigenvalues[h2[a,b]+w[ϵ]][[2]]

Plot[Eigenvalues[{h2[1,2]+w[ϵ]}],{ϵ,-1,1},Frame->True,FrameLabel->{"ϵ","Eigenvalues"},PlotStyle->Thick,BaseStyle->{FontSize->15}]

Plot3D[{Re[Eigenvalues[h2[1,2]+w[ϵx+I*ϵy]]][[1]],Re[Eigenvalues[h2[1,2]+w[ϵx+I*ϵy]]][[2]]},{ϵx,-2,2},{ϵy,-2,2},PlotStyle->Opacity[0.5],Mesh->None]

Plot3D[{Im[Eigenvalues[h2[1,2]+w[ϵx+I*ϵy]]][[1]],Im[Eigenvalues[h2[1,2]+w[ϵx+I*ϵy]]][[2]]},{ϵx,-2,2},{ϵy,-2,2},PlotStyle->Opacity[0.5],Mesh->None]

Degenerate perturbation theory

When a quantum system with degenerate eigenvalues is perturbed, the new eigenvalues and eigenstates can be calculated by matrix methods.  Matrix methods always work - also for non-degenerate perturbation problems! The general formulation is to find the eigenvalues and eigenvectors of the Hamiltonian (Subscript[H, 0])+W, where (Subscript[H, 0]) is the exactly solved Hamiltonian, and W is the perturbation.

ClearAll[h0,w]
h0={{1,0,0},{0,1,0},{0,0,1.5}};
MatrixForm[h0]
Eigensystem[h0]//MatrixForm
w[ϵ_]:=ϵ*{{0,0.1,0.2},{0.1,0,0.1},{0.2,0.1,0}};
MatrixForm[w[1]]
Plot[Eigenvalues[h0+w[ϵ]],{ϵ,0,1},Frame->True,FrameLabel->{"ϵ","Eigenvalues"},PlotStyle->Thick,BaseStyle->{FontSize->15}]
Eigensystem[h0+w[0]]//MatrixForm

Example of a perturbation to a particle in a Box

(* Particle in a Box & Perturbation *)
Manipulate[Plot[Sin[n*(π/1)*x],{x,0,1},PlotRange->{-1.2,1.2},Frame->True,FrameLabel->{"x/L","(Subscript[ψ, n])(x)"},Filling->Axis,GridLines->{{0.5},{}},PlotStyle->Thick,BaseStyle->{FontSize->15}],{n,1,30,1}]
Manipulate[Plot[((Sin[n*(π/1)*x])^2),{x,0,1},PlotRange->{-1.2,1.2},Frame->True,FrameLabel->{"x/L","(Subscript[ψ, n]^2)(x)"},Filling->Axis,GridLines->{{0.5},{}},PlotStyle->Thick,BaseStyle->{FontSize->15}],{n,1,30,1}]
Manipulate[Plot[W0*(UnitStep[x-((1/2)-(0.1/2))]-UnitStep[x-((1/2)+(0.1/2))]),{x,0,1},PlotRange->{-1.2,1.2},Frame->True,Filling->Axis,GridLines->{{0.5},{}},FrameLabel->{"x/L","(Subscript[W, 0])(x)"},PlotStyle->Thick,BaseStyle->{FontSize->15}],{W0,0.2,1}]

(* Matrix elements, perturbation terms in energy *)
Eu[n_,L_]=((((n/L)))^2) (* Unperturbed Eigenvalues *)
Wdiag[a_,L_,W0_,m_]=FullSimplify[((2*W0)/L)*(∫_((L-a)/2)^((L+a)/2))Sin[m*(π/L)*x]Sin[m*(π/L)*x]\[DifferentialD]x,{m ∈ Integers}]
Woffdiag[a_,L_,W0_,m_,n_]=FullSimplify[((2*W0)/L)*(∫_((L-a)/2)^((L+a)/2))Sin[m*(π/L)*x]Sin[n*(π/L)*x]\[DifferentialD]x,{m ∈ Integers, n ∈ Integers}]
(* Wdiag[a_,L_,W0_,m_], Woffdiag[a_,L_,W0_,m_,n_] *)
DiscretePlot[Eu[m,1],{m,1,5},Frame->True,FrameLabel->{"m","(Subscript[E, m])/(Subscript[E, 1])"},PlotStyle->Thick,BaseStyle->{FontSize->15}]
DiscretePlot[Wdiag[0.1,1,1,m],{m,1,50},Frame->True,PlotRange->{0,0.3},FrameLabel->{"m","(Subscript[W, mm])"},PlotStyle->Thick,BaseStyle->{FontSize->15}]
DiscretePlot[Woffdiag[0.1,1,1,1,m],{m,2,50},Frame->True,PlotRange->{-0.2,0.2},FrameLabel->{"m","(Subscript[W, m1])"},PlotStyle->Thick,BaseStyle->{FontSize->15}]
DiscretePlot[-(((Woffdiag[0.1,1,1,1,m])^2)/(1-(m^2))),{m,2,20},Frame->True,PlotRange->All,FrameLabel->{"m","-(((Subscript[W, m1])^2)/((Subscript[E, 1]) - (Subscript[E, m])))"},PlotStyle->Thick,BaseStyle->{FontSize->15}]

FullSimplify[((2*W0)/L)*(∫_((L-a)/2)^((L+a)/2))Sin[1*(π/L)*x]Sin[n*(π/L)*x]\[DifferentialD]x,{ n ∈ Integers}]
Series[((Sin[x])/x),{x,0,6}]

(* Brillouin-Wigner Analytical Solution *)
Solve[(x-e1-e0)*(x-e3-e0)==(e0^2),x]

(* Perturbed Eignefunctions from Rayleigh-Schrodinger Theory *)
Eu[n_,L_]:=((((n/L)))^2);
psi1P[x_,a_,W0_,L_,Nmin_,Nmax_]:=(Sqrt[(2/L)])*(Sin[1*(π/L)*x]+(∑_(m=2*Nmin+1)^(2*Nmax+1))Sin[((m*π)/2)] *((2*a*W0)/(L*Eu[1,L]))*((Sin[m*(π/L)*x])/((1^2)-(m^2)-((2*a*W0)/(L*Eu[1,L])))));
N1[L_]:=NIntegrate[((((Sqrt[(2/L)])*(Sin[1*(π/L)*x])))^2),{x,0,1}]
Np[a_,W0_,L_,Nmin_,Nmax_]:=NIntegrate[((psi1P[x,a,W0,L,Nmin,Nmax])^2),{x,0,1}]
Plot[{(((Sqrt[(2/1)])*(Sin[1*(π/1)*x]))/(Sqrt[N1[1]])),((psi1P[x,0.1,0.5,1,1,5])/(Sqrt[Np[0.1,0.5,1,1,5]])),((psi1P[x,0.1,1,1,1,5])/(Sqrt[Np[0.1,1,1,1,5]]))},{x,0,1},GridLines->{{0.5},{}},Frame->True,PlotStyle->{{Black},{Gray},{Black,Dashed}},BaseStyle->{FontSize->15},FrameLabel->{"x/L","(Subscript[ψ, 1])(x),ψ(x)"}]
Plot[{((((((Sqrt[(2/1)])*(Sin[1*(π/1)*x]))/(Sqrt[N1[1]]))))^2),(((((psi1P[x,0.1,0.5,1,1,5])/(Sqrt[Np[0.1,0.5,1,1,5]]))))^2),(((((psi1P[x,0.1,1,1,1,5])/(Sqrt[Np[0.1,1,1,1,5]]))))^2)},{x,0,1},GridLines->{{0.5},{}},Frame->True,PlotStyle->{{Black},{Gray},{Black,Dashed}},BaseStyle->{FontSize->15},FrameLabel->{"x/L","(ψ^2)(x)"}]

(* Matrix perturbation theory *)
Eu[n_,L_]:=((((n/L)))^2);
HpW[a_,L_,W0_]={{Eu[1,1]+Wdiag[a,L,W0,1],Woffdiag[a,L,W0,1,2],Woffdiag[a,L,W0,1,3]},{Woffdiag[a,L,W0,2,1],Eu[2,1]+Wdiag[a,L,W0,2],Woffdiag[a,L,W0,1,3]},{Woffdiag[a,L,W0,3,1],Woffdiag[a,L,W0,3,2],Eu[3,1]+Wdiag[a,L,W0,3]}};
MatrixForm[HpW[0.1,1,0]]
MatrixForm[HpW[0.1,1,0.5]]
Eigensystem[HpW[0.1,1,0.1]]
Plot[{Eigensystem[HpW[0.1,1,W0]][[1]][[1]],Eigensystem[HpW[0.1,1,W0]][[1]][[2]],Eigensystem[HpW[0.1,1,W0]][[1]][[3]]},{W0,0,30},Frame->True,FrameLabel->{"(Subscript[W, 0])(eV)","Eigenvalues ((Subscript[E, 1]))"},BaseStyle->{FontSize->15},AspectRatio->GoldenRatio]

(* How to get numbers from Mathematica Lists *)
Eigensystem[HpW[0.01,1,1]]
Eigensystem[HpW[0.01,1,1]][[2]]
Eigensystem[HpW[0.01,1,1]][[2]][[1]][[1]]
Eigensystem[HpW[0.01,1,1]][[2]][[3]][[1]]

(* Perturbed Eigenfunctions from the Matrix method *)
(* Ground State m=1 *)
psi1Pmatrix[x_,a_,L_,W0_]:=Eigensystem[HpW[a,L,W0]][[2]][[3]][[1]]*(Sqrt[(2/L)])*(Sin[1*(π/L)*x])+Eigensystem[HpW[a,L,W0]][[2]][[3]][[2]]*(Sqrt[(2/L)])*(Sin[2*(π/L)*x])+Eigensystem[HpW[a,L,W0]][[2]][[3]][[3]]*(Sqrt[(2/L)])*(Sin[3*(π/L)*x])
Plot[{(Sqrt[(2/1)])*(Sin[1*(π/1)*x]),psi1Pmatrix[x,0.1,1,2]},{x,0,1},GridLines->{{0.5},{}},Frame->True,PlotStyle->{{Black},{Black,Dashed}},BaseStyle->{FontSize->15},FrameLabel->{"x/L","(Subscript[ψ, 1])(x),ψ(x)"}]
Plot[{((((Sqrt[(2/1)])*(Sin[1*(π/1)*x])))^2),((psi1Pmatrix[x,0.1,1,2])^2)},{x,0,1},GridLines->{{0.5},{}},Frame->True,PlotStyle->{{Black},{Black,Dashed}},BaseStyle->{FontSize->15},FrameLabel->{"x/L","(ψ^2)(x)"}]

(* Perturbed Eigenfunctions from the Matrix method *)
(* State m=2 *)
psi2Pmatrix[x_,a_,L_,W0_]:=Eigensystem[HpW[a,L,W0]][[2]][[2]][[1]]*(Sqrt[(2/L)])*(Sin[1*(π/L)*x])+Eigensystem[HpW[a,L,W0]][[2]][[2]][[2]]*(Sqrt[(2/L)])*(Sin[2*(π/L)*x])+Eigensystem[HpW[a,L,W0]][[2]][[2]][[3]]*(Sqrt[(2/L)])*(Sin[3*(π/L)*x])
Plot[{(Sqrt[(2/1)])*(Sin[2*(π/1)*x]),psi2Pmatrix[x,0.1,1,2]},{x,0,1},GridLines->{{0.5},{}},Frame->True,PlotStyle->{{Black},{Black,Dashed}},BaseStyle->{FontSize->15},FrameLabel->{"x/L","(Subscript[ψ, 1])(x),ψ(x)"}]
Plot[{((((Sqrt[(2/1)])*(Sin[2*(π/1)*x])))^2),((psi2Pmatrix[x,0.1,1,2])^2)},{x,0,1},GridLines->{{0.5},{}},Frame->True,PlotStyle->{{Black},{Black,Dashed}},BaseStyle->{FontSize->15},FrameLabel->{"x/L","(ψ^2)(x)"}]

(* Perturbed Eigenfunctions from the Matrix method *)
(* State m=3 *)
psi3Pmatrix[x_,a_,L_,W0_]:=Eigensystem[HpW[a,L,W0]][[2]][[1]][[1]]*(Sqrt[(2/L)])*(Sin[1*(π/L)*x])+Eigensystem[HpW[a,L,W0]][[2]][[1]][[2]]*(Sqrt[(2/L)])*(Sin[2*(π/L)*x])+Eigensystem[HpW[a,L,W0]][[2]][[1]][[3]]*(Sqrt[(2/L)])*(Sin[3*(π/L)*x])
Plot[{(Sqrt[(2/1)])*(Sin[3*(π/1)*x]),psi3Pmatrix[x,0.1,1,2]},{x,0,1},GridLines->{{0.5},{}},Frame->True,PlotStyle->{{Black},{Black,Dashed}},BaseStyle->{FontSize->15},FrameLabel->{"x/L","(Subscript[ψ, 1])(x),ψ(x)"}]
Plot[{((((Sqrt[(2/1)])*(Sin[3*(π/1)*x])))^2),((psi3Pmatrix[x,0.1,1,2])^2)},{x,0,1},GridLines->{{0.5},{}},Frame->True,PlotStyle->{{Black},{Black,Dashed}},BaseStyle->{FontSize->15},FrameLabel->{"x/L","(ψ^2)(x)"}]

WKB method

The Schrodinger equation reads ψ’’=Qψ, where Q=((2m(V-E))/(ℏ^2)).  We look for solutions of the form ψ=(e^S), where we assume S=(1/ϵ)∑(Subscript[S, n]) (ϵ^n)

S[eps_,x_]:=Exp[(1/eps)*(S0[x]+S1[x]*eps+S2[x]*(eps^2)+S3[x]*(eps^3))];
Collect[(((eps^2)*D[S[eps,x],{x,2}])/(S[eps,x]))//Simplify,eps]

Module II: Bands, Doping, and Heterostructures

8. Electrons in a Crystal get their Bands, Gaps and Masses

Free Electron Bandstructure in 1D & Bandgap Opening

(* Free Electron Bandstructure in 1D *)
h=(6.63/(2*π))*(10^(-34)) (* Reduced Planck's constant, J.s *);
q=1.6*(10^(-19)) (* Electron Charge, Coulomb *);
me=9.1*(10^(-31)) (* Free electron mass, kg *);
E0[k_]=(((h^2)*(((k*(10^10)))^2))/(2*me*q))  (* in eV, k is in 1/m, of the order of 2π/a, a is in ~angstroms *);
Plot[E0[k],{k,-2,2},Frame->True,FrameLabel->{"k ((10^10)/m)","Energy (eV)"},GridLines->Automatic,PlotStyle->Thick,BaseStyle->{FontSize->15},AspectRatio->GoldenRatio]

(* Matrix elements *)
Hkk[n_,m_,L_]=(1/L)*(∫_0^L)Exp[I*((2π)/L)*(n-m)*x]\[DifferentialD]x
Hkk[1,2,1]

(* Periodic Potential in 1D, opening of bandgaps *)
G=2; (* (10^10) 1/m, reciprocal lattice vector *)
F=E0[G/2];
Eprime[k_,G_]=E0[k-(G/2)];
Epos[k_,G_,UG_]=F+Eprime[k,G]+(Sqrt[4*F*Eprime[k,G]+(((UG))^2)]);
Eneg[k_,G_,UG_]=F+Eprime[k,G]-(Sqrt[4*F*Eprime[k,G]+(((UG))^2)]);
Manipulate[Plot[{UnitStep[k-1]*Epos[k,2,UG]+UnitStep[-k+1]*Eneg[k,2,UG],E0[k]},{k,0,2},Frame->True,FrameLabel->{"k ((10^10)/m)","Energy (eV)"},GridLines->Automatic,PlotStyle->Thick,PlotRange->{-1,10},BaseStyle->{FontSize->15},AspectRatio->GoldenRatio],{UG,0,2}]

(* Section 14.4: Higher Order Interactions *)
(* 2x2 matrix, -G/2 <-> +G/2 interaction *)
m22[F1_,UG_]={{F1,-UG},{-UG,F1}}
MatrixForm[m22[F1,UG]]
Eigenvalues[m22[F1,UG]]
Manipulate[Eigenvalues[m22[4.0,UG]],{UG,0,1}]

(* 4x4 matrix, -3G/2 <-> -G/2 <-> +G/2 <-> +3G/2 interaction *)
m44[F1_,UG_]={{9*F1,-UG,0,0},{-UG,F1,-UG,0},{0,-UG,F1,-UG},{0,0,-UG,9*F1}}
MatrixForm[m44[F1,UG]]
Eigenvalues[m44[F1,UG]]
Manipulate[Eigenvalues[m44[4.0,UG]],{UG,0,1}]

1D Bandstructure from the Brillouin-Wigner theory

(* Brilloiun-Wigner perturbation model for 1D electrons in W(x)=-2*UG*Cos[Gx] *)
(* Solve[x-E0[k]-2*(UG^2)*((x-E0[k]-4*F)/((((x-E0[k]-4*F))^2)-16*E0[k]*F))==0,x]; *)
(* Solve[x-A-2*(UG^2)*((x-A-4*B)/((((x-A-4*B))^2)-16*A*B))==0,x]; *)
h=(6.63/(2*π))*(10^(-34)) (* Reduced Planck's constant, J.s *);
q=1.6*(10^(-19)) (* Electron Charge, Coulomb *);
me=9.1*(10^(-31)) (* Free electron mass, kg *);
E0[k_]=(((h^2)*(((k*(10^10)))^2))/(2*me*q))  (* in eV, k is in 1/m, of the order of 2π/a, a is in ~angstroms *);
G=2; (* (10^10) 1/m, reciprocal lattice vector *)
F=E0[G/2];
Esol1[A_,B_,UG_]=Solve[x-A-2*(UG^2)*((x-A-4*B)/((((x-A-4*B))^2)-16*A*B))==0,x][[1]][[1]][[2]];
Esol2[A_,B_,UG_]=Solve[x-A-2*(UG^2)*((x-A-4*B)/((((x-A-4*B))^2)-16*A*B))==0,x][[2]][[1]][[2]];
Esol3[A_,B_,UG_]=Solve[x-A-2*(UG^2)*((x-A-4*B)/((((x-A-4*B))^2)-16*A*B))==0,x][[3]][[1]][[2]];
Ebw1[k_,UG_]=Esol1[E0[k],F,UG];
Ebw2[k_,UG_]=Esol2[E0[k],F,UG];
Ebw3[k_,UG_]=Esol3[E0[k],F,UG];
Manipulate[Plot[{E0[k],E0[k+G],E0[k-G],Ebw1[k,UG],Ebw2[k,UG],Ebw3[k,UG]},{k,-1.5,1.5},Frame->True,FrameLabel->{"k ((10^10)/m)","Energy (eV)"},GridLines->Automatic,PlotRange->{-2,20},PlotStyle->{{Gray,Dashed},{Gray,Dashed},{Gray,Dashed},{Black,Thick},{Black,Thick},{Black,Thick}},BaseStyle->{FontSize->20},AspectRatio->1.5],{UG,0,3}]

(* Compare energy of state |0> from Brillouin-Wigner Theory *)(* with results from Degenerate perturbation theory for states |G> <-> |0> <-> |-G> *) 
m33[F1_,UG_]={{4*F1,-UG,0},{-UG,0,-UG},{0,-UG,4*F1}}
MatrixForm[m33[F1,UG]]
Eigenvalues[m33[F1,UG]]
Manipulate[Eigenvalues[m33[F1,UG]],{F1,3.9,4.1},{UG,0,1}]

Periodic Potentials in Real and k-Space

(* Example of a Fourier Transform in Mathematica *)
Needs["FourierSeries`"]
f[x_]:=UnitStep[1-x] UnitStep[x+1]
Plot[f[x],{x,-3,3}]
Plot[NFourierTransform[f[x],x,k],{k,-20,20},PlotRange->All]

v0[x_,n_]:=-(1/(n^2))*Cos[n*2*π*x];  (* Example sinusoid periodic potential *)
vper[x_,N_]:=(∑_(n=1)^N)v0[x,n]; (* Example non-sinusoid periodic potential *)
(*vksp[k_,N_]:=NIntegrate[Exp[I*2*π*k*x]*vper[x,N],{x,-40,40}]; (* Fourier Transform of non-sinusoid periodic potential *)*)
Plot[vper[x,10],{x,-3,3},PlotRange->All,Frame->True,Filling->Bottom,BaseStyle->{FontSize->15},FrameLabel->{"x","V(x)"}]
(*Needs["FourierSeries`"]*)
(*Plot[NFourierTransform[vper[x,10],x,k],{k,-15,15}]*)(* Fourier Transform of non-sinusoid periodic potential using predef FT *)

Plot[-Convolve[FourierTransform[vper[x,5],x,k],(1/(Sqrt[π*(0.1^2)]))*Exp[-((((k/0.1)))^2)],k,kp],{kp,-20*π,20*π},PlotRange->All,FrameLabel->{"k","V(k)"},Frame->True,Filling->Bottom,BaseStyle->{FontSize->15}]

v2d[x_,y_,n_]:=-(1/(n^2))*Cos[n*2*π*x]*Cos[n*2*π*y];  (* Example 2D sinusoid periodic potential *)
vper2d[x_,y_,N_]:=(∑_(n=1)^N)v2d[x,y,n]; (* Example non-sinusoid periodic potential *)
(*vksp[k_,N_]:=NIntegrate[Exp[I*2*π*k*x]*vper[x,N],{x,-40,40}]; (* Fourier Transform of non-sinusoid periodic potential *)*)
Plot3D[vper2d[x,y,10],{x,-3,3},{y,-3,3},PlotRange->All,BaseStyle->{FontSize->15}]
(*Needs["FourierSeries`"]*)
(*Plot[NFourierTransform[vper[x,10],x,k],{k,-15,15}]*)(* Fourier Transform of non-sinusoid periodic potential using predef FT *)

Module III: Quantum Electronic Devices

Module IV: Quantum Photonic Devices

Appendix I: The Mathematica Toolbox

0. Basics of Mathematica I

q=1.6*(10^(-19))  (* Electron charge in Coulombs, how to Define a constant and comment *)

f[x_]:=(x^2)-x+1;  (* How to define a function *)

Solve[f[x]==0,x]  (* How to solve a polynomial equation *)

Solve[f[x]==0,x][[1]][[1]][[2]] (* How to get the 1st solution of the equation *)

Solve[f[x]==0,x][[2]][[1]][[2]] (* How to get the 2nd solution of the equation *)

Plot[{x,Cos[x]},{x,0,2*π},Frame->True,FrameLabel->{"x","x and cos[x]"}, BaseStyle->{FontSize->15}]  (* How to plot and label *)

lis[x_]:={x,(x^2),Sin[x]};(* How to plot lists, and Sort them *)
N[lis[2]]
Plot[N[lis[x]],{x,-2,2},Frame->True,FrameLabel->{"x","Functions"}, BaseStyle->{FontSize->15}]
Sort[N[lis[2]]]
Plot[{Sort[N[lis[x]]][[1]],Sort[N[lis[x]]][[2]],Sort[N[lis[x]]][[3]]},{x,-2,2},PlotStyle->{Red,Green,Blue},Frame->True,FrameLabel->{"x","Sorted Functions"}, BaseStyle->{FontSize->15}]

FindRoot[Cos[x]==x,{x,0}]  (* How to find numerical roots of difficult equations *)
FindRoot[Cos[x]==x,{x,0}][[1]][[2]]  
FindRoot[Cos[x]==x,{x,0}][[1,2]]  (* How to get the numerical value of the root *)

(* How to Use FindRoot to find Multiple And Variable Number of Roots and Sort them *)
f[x_]:=Sin[10*x];
Plot[f[x],{x,0,1}]
FindRoot[f[x],{x,0,1}] (* Ordinary FindRoot *)
Sort@DeleteDuplicates[Select[FindRoot[f[x],{x,Range[0,1,0.1]}][[1,2]],(0≤#<1)&],Abs[#2-#1]<(10^(-8))&]

∫(Sqrt[Tan[x]])\[DifferentialD]x  (* How to evaluate an indefinite integral analytically *)
(∫_0^∞)(x^2)*Exp[-((((x/a)))^2)]\[DifferentialD]x  (* How to evaluate a definite integral analytically *)
NIntegrate[(x^2)*Exp[-(x^2)],{x,0,1000}] (* How to evaluate any integral numerically *)
(∫_0^(π/2))Sin[θ]*Cos[θ]\[DifferentialD]θ (* A Simple Definite Integal *)

(* How to define conditional functions *)
f[d_,x_]:=Switch[d,1,Sin[x],2,Cos[x],3,Tan[x]];
Plot[{f[1,x],f[2,x],f[3,x]},{x,-5,5},PlotRange->{-2,2},PlotStyle->{Red,Blue,Green}]

0. Basics of Mathematica 2

A={{1,2},{0,1}}  (* How to define a simple matrix *)
MatrixForm[A]  (* How to display a matrix in the normal format *)

B={{2,1},{1,2}}  (* How to define a simple matrix *)
MatrixForm[A.B] (* Product of two matrices *)
MatrixForm[B.A] (* Product of two matrices; note A.B≠B.A *)

MatrixForm[MatrixPower[A,20]]  (* How to take powers of matrices *)
MatrixForm[N[MatrixFunction[Sin,A]]]  (* How to take numerical functions of matrices *)

Eigenvalues[A]  (* How to find eigenvalues of matrices *)
Eigenvectors[A]  (* How to find eigenvectors of matrices *)

MatrixForm[IdentityMatrix[10]] (* How to display a 10x10 identity matrix *)

tightbinding[center_,left_,right_]=center*IdentityMatrix[10]+SparseArray[{Band[{1,2}]->left,Band[{2,1}]->right},{10,10}] (* How to create a 10x10 tight-binding matrix *)
MatrixForm[tightbinding[1.0,-0.1,-0.1]] (* How to display the 10x10 tight-binding matrix *)
MatrixPlot[tightbinding[1.0,-0.1,-0.1]] (* How to "Plot" the 10x10 tight-binding matrix *)
ListPlot[Eigenvalues[tightbinding[1.0,-0.1,-0.1]]]  (* How to Plot the eigenvalues of the 10x10 tight-binding matrix *)

tightbindingBig[center_,left_,right_]=center*IdentityMatrix[100]+SparseArray[{Band[{1,2}]->left,Band[{2,1}]->right},{100,100}] ;(* How to create a 100x100 tight-binding matrix *)
Manipulate[Quiet[ListPlot[Eigenvalues[tightbindingBig[1.0,left,left]],PlotRange->{0,2}]]  ,{left,0,-0.3}](* How to Plot the eigenvalues of the 100x100 tight-binding matrix with a tunable hopping parameter.  This is exactly how electrons develop a bandstructure in solids! *)

Figures and plots of shapes in Mathematica

Spheres

ℛ=Sphere[{0,0,0},1];
ContourPlot3D[Evaluate@RegionDistance[ℛ,{x,y,z}],{x,-2,2},{y,0,2},{z,-2,2},Mesh->None,Contours->{0.25,0.5,1},BoxRatios->Automatic]

A=ContourPlot3D[Evaluate@RegionDistance[ℛ,{x,y,z}],{x,0,2},{y,-2,2},{z,-2,2},Mesh->None,Contours->{1},BoxRatios->Automatic,Axes->None]

B=ContourPlot3D[Evaluate@RegionDistance[ℛ,{x,y,z}],{x,-2,0},{y,-2,2},{z,-2,2},Mesh->None,Contours->{0.07},BoxRatios->Automatic]

Show[A,B]

Graphics3D[{Opacity[0.7],Sphere[]}]

Graphics3D[{Opacity[0.7],Sphere[{0,0,0},1]}]

Graphics3D[{Opacity[0.5],FaceForm[Yellow,Blue],Sphere[{0,0,0},1]},PlotRange->{{-1,1},{-1,1},{-1,1}}]

Manipulate[Graphics3D[{Opacity[0.7],FaceForm[Yellow,Blue],Sphere[{0,0,0},0.5]},PlotRange->{{-1,0},{-1,1},{-1,1}}]+Graphics3D[{Opacity[0.7],FaceForm[Yellow,Blue],Sphere[{0,0,0},0.5]},PlotRange->{{0,1},{-1,1},{-1,1}}],{r,0,1}]