LASER RATE EQUATIONS by P N V Anil Kumar,IIT Bombay

Solves the LASER RATE EQUATIONS

%Transient Response of Coupled Rate Equations

% Copyright .This code is developed by P N V Anil Kumar,IIT Bombay & %SAC/ISRO,Ahmedabad

clc;

clear all;

close all;

tic

Time_initial = 0;

Time_final = 10^-9;

Time_steps = 10^-16;

Time = Time_initial:Time_steps:Time_final;

Frequency = 1.2*10^10;

N = zeros(1,length(Time));

I = 2.0*10^-2;

S = zeros(1,length(Time));

q = 1.6*10^-19;

V = 2*10^-11;

tau_n = 10^-9;

g0 = 1.5*10^-5;

Nth = 10^18;

epsilon = 1.5*10^-17;

tau_p = 10^-12;

beta = 10^-4;

gamma = 0.2;

alpha = 3.1;

phi = zeros(1,length(Time));

Optical_Output_Power_of_LASER = zeros(1,length(Time));

Quantum_Efficiency = [0.1 0.5 0.7 1];

h = 6.6261*10^-27;%in CGS unit

I_bias = linspace(10,50,10^4).*10^-3;%Bias Current in mA

Gain_Saturation_Parameter = 3.4*10^-17;

Relaxation_Oscillation_Frequency = zeros(1,length(I_bias));

Laser_Transfer_Function_Magnitude = zeros(1,length(Relaxation_Oscillation_Frequency));

Laser_Transfer_Function_Phase = zeros(1,length(Relaxation_Oscillation_Frequency));

for i =1:length(Time)

x1(1,i) = I/(q*V);

x2(1,i) = N(1,i)/tau_n;

x3(1,i) = (g0*(N(1,i) - Nth));

x4(1,i) = 1+(epsilon*S(1,i));

x5(1,i) = (x3(1,i)/x4(1,i))*S(1,i);

x6(1,i) = x1(1,i)-x2(1,i)- x5(1,i);

N(1,i+1) = N(1,i) + (Time_steps*x6(1,i));

y1(1,i) = (gamma*x3(1,i))-(1/tau_p);

y2(1,i) = (1/2)*alpha* y1(1,i);

phi(1,i+1) = phi(1,i) + (Time_steps*y2(1,i));

x7(1,i) = (((gamma*g0)*(N(1,i)-Nth))/x4(1,i))*S(1,i);

x8(1,i) = S(1,i)/tau_p;

x9(1,i) = (gamma*beta*N(1,i))/tau_n;

x10(1,i) = x7(1,i)- x8(1,i)+x9(1,i);

S(1,i+1) = S(1,i) + (Time_steps*x10(1,i));

z1(1,i) = S(1,i)*V*h*Frequency;

z2(1,i) = 2*gamma*tau_p;

Optical_Output_Power_of_LASER(1,i) = z1(1,i)/z2(1,i);

end

Pump_Rate = linspace(0,2,10^4);

tau_L = 1*10^-9;

for j = 1:length(Pump_Rate)

if Pump_Rate(j) < 1

Carrier_Concentration(j) = Pump_Rate(j)*Nth;

Photon_Concentration(j) = Pump_Rate(j)/(1-Pump_Rate(j));

else

Carrier_Concentration(j)=Nth;

Photon_Concentration(j) = (tau_p/tau_L)*Nth*(Pump_Rate(j)-1);

end

end

for k =1:length(I_bias)

z3 = (gamma*g0)/(q*V);

z4(1,k) = z3*(I_bias(1,k)-I);

z5 = (gamma*Gain_Saturation_Parameter*tau_p)/(q*V);

z6(1,k) = z5*(I_bias(1,k)-I);

z7(1,k) = 1-z6(1,k);

Capital_K(1,k) = z4(1,k)*z7(1,k);

z8 = z3*(tau_p+(Gain_Saturation_Parameter/g0));

z9(1,k) = z8*(I_bias(1,k)-I);

gamma_d(1,k) = (1/tau_n)+z9(1,k)*z7(1,k);

z10(1,k) = (1/2)*(gamma_d(1,k)^2);

z11(1,k) = Capital_K(1,k)-z10(1,k);

z12(1,k) = real(sqrt(z11(1,k)));

Relaxation_Oscillation_Frequency(1,k) = (1/(2*pi))*z12(1,k);

end

I_bias1 = [30 40 50]*10^-3;

clear j;

for ii = 1:length(I_bias1)

for jj = 1:length(Relaxation_Oscillation_Frequency)

pp1 = (gamma*g0)/(q*V);

pp2(ii,jj) = pp1*(I_bias1(ii)-I);

pp3 = (gamma*Gain_Saturation_Parameter*tau_p)/(q*V);

pp4(ii,jj) = pp3*(I_bias1(ii)-I);

pp5(ii,jj) = 1-pp4(ii,jj);

Capital_K1(ii,jj) = pp2(ii,jj)* pp5(ii,jj);

pp6 = pp1*(tau_p+(Gain_Saturation_Parameter/g0));

pp7(ii,jj) = pp6*(I_bias1(ii)-I);

gamma_d1(ii,jj) = (1/tau_n)+pp7(ii,jj)*pp5(ii,jj);

pp8(ii,jj) = j*(2*pi*Relaxation_Oscillation_Frequency(1,jj));

pp9(ii,jj) = pp8(ii,jj)*(pp8(ii,jj)+gamma_d1(ii,jj));

pp10(ii,jj) = pp9(ii,jj)+Capital_K1(ii,jj);

Laser_Transfer_Function_Magnitude(ii,jj) =abs(Capital_K1(ii,jj)/pp10(ii,jj));

Laser_Transfer_Function_Phase(ii,jj) = angle(Capital_K1(ii,jj)/pp10(ii,jj));

end

end

figure(1)

ax = axes;

plot(Time,N(1:length(N)-1),'m','Linewidth',2)

grid minor

xlabel('Time in Seconds','Interpreter','latex','FontSize',13);

ylabel('Carrier Concentration','Interpreter','latex','FontSize',13);

title('LASER Coupled Rate Equation','Interpreter','latex','FontSize',13);

legend('Transient Carrier Concentration','Location','southeast');

ax.YAxis.Color = [0 0 0];

ax.XAxis.Color = [0 0 0];

figure(2)

ax = axes;

plot(Time,S(1:length(N)-1),'r','Linewidth',2)

grid minor

xlabel('Time in Seconds','Interpreter','latex','FontSize',13);

ylabel('Photon Concentration','Interpreter','latex','FontSize',13);

title('LASER Coupled Rate Equation','Interpreter','latex','FontSize',13);

legend('Transient Photon Concentration');

ax.YAxis.Color = [0 0 0];

ax.XAxis.Color = [0 0 0];

figure(3)

ax = axes;

plot(N,S,'b','Linewidth',2)

grid minor

xlabel('Carrier Concentration','Interpreter','latex','FontSize',13);

ylabel('Photon Concentration','Interpreter','latex','FontSize',13);

title('Photon Concentration Vs Carrier Concentration','Interpreter','latex','FontSize',13);

ax.YAxis.Color = [0 0 0];

ax.XAxis.Color = [0 0 0];

figure(4)

ax = axes;

yyaxis left

plot(Time,S(1:length(N)-1),'r','Linewidth',2)

xlabel('Time in Seconds','Interpreter','latex','FontSize',13);

ylabel('Photon Concentration','Interpreter','latex','FontSize',13);

yyaxis right

plot(Time,N(1:length(N)-1),'m','Linewidth',2)

xlabel('Time in Seconds','Interpreter','latex','FontSize',13);

ylabel('Carrier Concentration','Interpreter','latex','FontSize',13);

title('LASER Rate Equations - Transient Response','Interpreter','latex','FontSize',12);

grid minor

ax.YAxis(1).Color = [1 0 0];

ax.YAxis(2).Color = [1 0 1];

figure(5)

ax = axes;

yyaxis left

plot(Pump_Rate,Carrier_Concentration,'r','Linewidth',2)

xlabel('Pump Rate','Interpreter','latex','FontSize',13);

ylabel('Carrier Concentration','Interpreter','latex','FontSize',13);

ylim([0,11*10^17])

xline(1,'--','Color','[0.4 0.6 0.2]','Linewidth',2);

y1 = yline(Nth,'--','Nth','Color','[1 0 0]','Linewidth',2);

y1.LabelHorizontalAlignment = 'left';

yl.Color = [1 0 1];

yyaxis right

plot(Pump_Rate,Photon_Concentration,'m','Linewidth',2)

xlabel('Pump Rate','Interpreter','latex','FontSize',13);

ylabel('Photon Concentration','Interpreter','latex','FontSize',13);

xline(1,'--','Color','[0.4 0.6 0.2]','Linewidth',2);

title('Steady State Solution Vs Pump Rate','Interpreter','latex','FontSize',12);

ax.YAxis(1).Color = [1 0 0];

ax.YAxis(2).Color = [1 0 1];

grid minor

figure(6)

ax = axes;

plot(Time,phi(1:length(phi)-1),'r','Linewidth',2)

xlabel('Time in Seconds','Interpreter','latex','FontSize',13);

ylabel('Optical Phase','Interpreter','latex','FontSize',13);

title('LASER Coupled Rate Equation','Interpreter','latex','FontSize',13);

ax.YAxis.Color = [0 0 0];

ax.XAxis.Color = [0 0 0];

grid minor

figure(7)

ax = axes;

yyaxis left

plot(Time,N(1:length(N)-1),'r','Linewidth',2)

xlabel('Time in Seconds','Interpreter','latex','FontSize',13);

ylabel('Carrier Concentration','Interpreter','latex','FontSize',13);

yyaxis right

plot(Time,phi(1:length(phi)-1),'m','Linewidth',2)

xlabel('Time in Seconds','Interpreter','latex','FontSize',13);

ylabel('Optical Phase','Interpreter','latex','FontSize',13);

title('LASER Coupled Rate Equation','Interpreter','latex','FontSize',13);

ax.YAxis(1).Color = [1 0 0];

ax.YAxis(2).Color = [1 0 1];

grid minor

figure(8)

ax = axes;

yyaxis left

plot(Time,S(1:length(N)-1),'r','Linewidth',2)

xlabel('Time in Seconds','Interpreter','latex','FontSize',13);

ylabel('Photon Concentration','Interpreter','latex','FontSize',13);

yyaxis right

plot(Time,phi(1:length(phi)-1),'m','Linewidth',2)

xlabel('Time in Seconds','Interpreter','latex','FontSize',13);

ylabel('Optical Phase','Interpreter','latex','FontSize',13);

title('LASER Coupled Rate Equation','Interpreter','latex','FontSize',13);

ax.YAxis(1).Color = [1 0 0];

ax.YAxis(2).Color = [1 0 1];

grid minor

figure(9)

ax = axes;

plot(Time,Optical_Output_Power_of_LASER.*Quantum_Efficiency(1),'r','Linewidth',2)

hold on

plot(Time,Optical_Output_Power_of_LASER.*Quantum_Efficiency(2),'g','Linewidth',2)

hold on

plot(Time,Optical_Output_Power_of_LASER.*Quantum_Efficiency(3),'b','Linewidth',2)

hold on

plot(Time,Optical_Output_Power_of_LASER.*Quantum_Efficiency(4),'m','Linewidth',2)

hold off

xlabel('Time in Seconds','Interpreter','latex','FontSize',13);

ylabel('Optical Output Power of LASER','Interpreter','latex','FontSize',13);

title('Optical Output Power of LASER Vs Time','Interpreter','latex','FontSize',13);

legend('Quantum Efficiency = 0.1','Quantum Efficiency = 0.5','Quantum Efficiency = 0.7','Quantum Efficiency = 1','Interpreter','latex','FontSize',11,'Location','northeast');

ax.YAxis.Color = [0 0 0];

ax.XAxis.Color = [0 0 0];

ylim([0,22.5])

grid minor

figure(10)

ax = axes;

plot(I_bias.*10^2,Relaxation_Oscillation_Frequency,'m','Linewidth',2)

xlabel('$I_{bias}$ (in mA)','Interpreter','latex','FontSize',13);

ylabel('Relaxation Oscillation Frequency (Hz)','Interpreter','latex','FontSize',13);

title('Relaxation Oscillation Frequency Vs $I_{bias}$','Interpreter','latex','FontSize',13);

grid minor

ax.YAxis.Color = [0 0 0];

ax.XAxis.Color = [0 0 0];

x1 = xline(I*10^2,'--','Threshold Current','Color','[0 0 1]','Linewidth',2);

x1.Interpreter = 'latex';

x1.Color = [0 0 1];

x1.FontSize =11;

M_0_e = 9.1*10^-31;

%Mn = 1.062*M_0_e;%4K Effective Mass

Mn = 1.182*M_0_e;%300K Effective Mass

T =300;

Kb = 1.38*10^-23;

x1 = (Kb*T);

h = 6.634*10^-34;

hbar = h/(2*pi);

deltaE_e = (linspace(0,6*x1,10^5));

gc = (Mn.*sqrt(2.*Mn.*deltaE_e))/((pi^2)*(hbar^3));

x = 1+exp((deltaE_e+(1.6*10^-20))./(Kb*T));

f = 1./x;

y = (gc.*f)*(10^-6)*(1.6*10^-19);

figure(11)

plot(y,(deltaE_e/(1.6*10^-19)),'b','linewidth',2)

xlabel('g_{C}(E) * f(E)');

ylabel('E-E_{C} in (eV)');

title('g_{C}(E) * f(E) Vs E-E_{C}');

yline((0.0259/2),'--r');

set(gca,'YTick',[(0.0259/2) 0.0259,0.0776,0.1294,0.1811,0.2329]);

set(gca,'YTickLabel',{'{KT}/_{2}','KT','3KT','5KT','7KT','9KT'});

grid minor

M_0_p = 9.1*10^-31;

%Mp =0.590*M_0_p;%4K Effective Mass

Mp =0.81*M_0_p;%300K Effective Mass

deltaE_p = (linspace(0,6*x1,10^5));

gv = (Mp.*sqrt(2.*Mp.*deltaE_p))/((pi^2)*(hbar^3));

x2 = 1+exp((deltaE_p+(1.6*10^-20))./(Kb*T));

f_p = 1./x2;

y2 = (gv.*f_p)*(10^-6)*(1.6*10^-19);

figure(12)

plot(y2,(deltaE_p/(1.6*10^-19)),'b','linewidth',2)

xlabel('g_{V}(E) * [1-f(E)]');

ylabel('E_{V} - E in (eV)');

title('g_{V}(E) * [1-f(E)] Vs E_{V} - E');

yline((0.0259/2),'--r');

set(gca,'YTick',[(0.0259/2) 0.0259,0.0776,0.1294,0.1811,0.2329]);

set(gca,'YTickLabel',{'{KT}/_{2}','KT','3KT','5KT','7KT','9KT'});

grid minor

figure(13)

ax = axes;

plot(y,(deltaE_e/(1.6*10^-19)),'b','linewidth',2)

hold on

plot(y2,(deltaE_p/(1.6*10^-19)),'m','linewidth',2)

hold off

set(gca,'YTick',[(0.0259/2) 0.0259,0.0776,0.1294,0.1811,0.2329]);

set(gca,'YTickLabel',{'{KT}/_{2}','KT','3KT','5KT','7KT','9KT'});

yline((0.0259/2),'--r');

xlabel('g_{C}(E) * f(E) & g_{V}(E) * [1-f(E)]');

ylabel('E-E_{C} in (eV) & E_{V} - E in (eV)');

title('Density of States (DOS) of CB & VB');

legend('g_{C}(E) * f(E) Vs E-E_{C} in (eV)','g_{V}(E) * [1-f(E)] Vs E_{V} - E in (eV)');

grid minor

ax.YAxis.Color = [0 0 0];

ax.XAxis.Color = [0 0 0];

Kb = 1.38*10^-23;

T = linspace(1,650,10^5);

M0 = 9.1*10^-31;

Mn = 1.182 *M0;

Mp = 0.81*M0;

h = 6.62607015 *10^-34;

ND = 10^15;

NA = 0;

EG_0K = 1.166;%eV For Silicon

alpha = 4.73*10^-4;%eV per Kelvin

beta = 636;%Kelvin

gD = 2;

EC_ED = 50*10^-3;%in eV

for i=1:length(T)

Kb_T_in_eV(i) = (Kb*T(i))/(1.6*10^-19);

x1(i) = (2*pi*Mn*Kb*T(i))/(h^2);

x2(i) = (x1(i))^(3/2);

Nc(i) = (2*x2(i))/(1*10^6);%in cm3

x3(i) = (2*pi*Mp*Kb*T(i))/(h^2);

x4(i) = (x3(i))^(3/2);

Nv(i) = (2*x4(i))/(1*10^6);%in cm3

x6(i) = alpha/(T(i)+beta);

x7(i) = (x6(i))* (T(i)^2);%in eV

Eg(i) = (EG_0K - x7(i));%in eV

x8(i) = sqrt(Nc(i)* Nv(i));%in cm3

x9(i) = -(Eg(i)/(2*Kb_T_in_eV(i)));%dimensionless

x10(i) = (x8(i))*exp(x9(i)); %in cm3

x11(i) = x10(i)/ND;%dimensionless

x12(i) = ((EC_ED)/Kb_T_in_eV(i));

N_zita(i) = ((Nc(i))/(gD))*(exp(-x12(i)));

x13(i) = 1+((4*ND)/(N_zita(i)));

x14(i) = (sqrt(x13(i)))-1;

n1(i) = ((N_zita(i))/2)*x14(i);

x15(i) = n1(i)/ND;

x16(i) = (((ND-NA)/2)^2) + ((x10(i))^2);

x17(i) = sqrt(x16(i));

n2(i) = ((ND-NA)/2)+ x17(i);

x18(i) = n2(i)/ND;

end

figure(14)

plot(T,Eg,'b','linewidth',2)

xlabel('Temperature in Kelvin');

ylabel('Band Gap in eV for Silicon');

title('Band Gap Vs Temperature For Silicon');

grid minor

Mn = 0.553 *M0;

Mp = 0.357*M0;

h = 6.62607015 *10^-34;

ND = 10^15;

NA = 0;

EG_0K = 0.7437;%eV

alpha = 4.77*10^-4;%eV per Kelvin

beta = 235;%Kelvin

gD = 2;

EC_ED = 12*10^-3;%in eV

for i=1:length(T)

Kb_T_in_eV(i) = (Kb*T(i))/(1.6*10^-19);

x1(i) = (2*pi*Mn*Kb*T(i))/(h^2);

x2(i) = (x1(i))^(3/2);

Nc(i) = (2*x2(i))/(1*10^6);%in cm3

x3(i) = (2*pi*Mp*Kb*T(i))/(h^2);

x4(i) = (x3(i))^(3/2);

Nv(i) = (2*x4(i))/(1*10^6);%in cm3

x6(i) = alpha/(T(i)+beta);

x7(i) = (x6(i))* (T(i)^2);%in eV

Eg(i) = (EG_0K - x7(i));%in eV

x8(i) = sqrt(Nc(i)* Nv(i));%in cm3

x9(i) = -(Eg(i)/(2*Kb_T_in_eV(i)));%dimensionless

x10(i) = (x8(i))*exp(x9(i)); %in cm3

x11(i) = x10(i)/ND;%dimensionless

x12(i) = ((EC_ED)/Kb_T_in_eV(i));

N_zita(i) = ((Nc(i))/(gD))*(exp(-x12(i)));

x13(i) = 1+((4*ND)/(N_zita(i)));

x14(i) = (sqrt(x13(i)))-1;

n1(i) = ((N_zita(i))/2)*x14(i);

x15(i) = n1(i)/ND;

x16(i) = (((ND-NA)/2)^2) + ((x10(i))^2);

x17(i) = sqrt(x16(i));

n2(i) = ((ND-NA)/2)+ x17(i);

x18(i) = n2(i)/ND;

end

figure(15)

plot(T,Eg,'b','linewidth',2)

xlabel('Temperature in Kelvin');

ylabel('Band Gap in eV for Germanium');

title('Band Gap Vs Temperature For Germanium');

grid minor

Eph = linspace(1.55,1.8,10^4);

T1 = 293;

T2 = 323;

Eg1 = 1.55;

Eg2 = 0.8;

Kb = 1.38*10^-23;

rr1 = (Eph-Eg1);

rr2 = -rr1/(0.0253);

rr3 = -rr1/(0.0279);

rr4 = (Eph-Eg2);

rr5 = -rr4/(0.0253);

rr6 = -rr4/(0.0279);

Power1 = rr1.*exp(rr2);

Power2 = rr1.*exp(rr3);

power3 = rr4.*exp(rr5);

power4 = rr4.*exp(rr6);

figure(16)

plot(Eph,Power1,'r','Linewidth',2)

hold on

plot(Eph,Power2,'m','Linewidth',2)

hold off

grid minor

xlabel('Photon Energy ($E_{ph})$','Interpreter','latex','FontSize',13);

ylabel('Normalized Power','Interpreter','latex','FontSize',13);

title('Emission Spectrum for $\lambda$ = 800 nm','Interpreter','latex','FontSize',13);

legend('T = 20 C','T = 50 C','Interpreter','latex','FontSize',11,'Location','northeast');

figure(17)

ax = axes;

plot(Relaxation_Oscillation_Frequency,Laser_Transfer_Function_Phase(1,:),'r','Linewidth',2)

hold on

plot(Relaxation_Oscillation_Frequency,Laser_Transfer_Function_Phase(2,:),'g','Linewidth',2)

hold on

plot(Relaxation_Oscillation_Frequency,Laser_Transfer_Function_Phase(3,:),'b','Linewidth',2)

hold off

xlabel('Relaxation Oscillation Frequency (Hz)','Interpreter','latex','FontSize',12);

ylabel('Laser Transfer Function Phase','Interpreter','latex','FontSize',12);

legend('$I_{bias}$ = 30 mA','$I_{bias}$ = 40 mA','$I_{bias}$ = 50 mA','Interpreter','latex','FontSize',12)

grid minor

ax.YAxis.Color = [0 0 0];

ax.XAxis.Color = [0 0 0];

figure(18)

ax = axes;

plot(Relaxation_Oscillation_Frequency,Laser_Transfer_Function_Magnitude(1,:),'r','Linewidth',2)

hold on

plot(Relaxation_Oscillation_Frequency,Laser_Transfer_Function_Magnitude(2,:),'g','Linewidth',2)

hold on

plot(Relaxation_Oscillation_Frequency,Laser_Transfer_Function_Magnitude(3,:),'b','Linewidth',2)

hold off

xlabel('Relaxation Oscillation Frequency (Hz)','Interpreter','latex','FontSize',12);

ylabel('Laser Transfer Function Magnitude','Interpreter','latex','FontSize',12);

legend('$I_{bias}$ = 30 mA','$I_{bias}$ = 40 mA','$I_{bias}$ = 50 mA','Interpreter','latex','FontSize',12)

grid minor

ax.YAxis.Color = [0 0 0];

ax.XAxis.Color = [0 0 1];

clear all;

toc