LASER RATE EQUATIONS by P N V Anil Kumar,IIT Bombay
Version 1.0.0 (4.26 KB) by
Nagavenkata Anil Peddiraju
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
Cite As
Nagavenkata Anil Peddiraju (2024). LASER RATE EQUATIONS by P N V Anil Kumar,IIT Bombay (https://www.mathworks.com/matlabcentral/fileexchange/109540-laser-rate-equations-by-p-n-v-anil-kumar-iit-bombay), MATLAB Central File Exchange. Retrieved .
MATLAB Release Compatibility
Created with
R2022a
Compatible with any release
Platform Compatibility
Windows macOS LinuxTags
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!Discover Live Editor
Create scripts with code, output, and formatted text in a single executable document.
Version | Published | Release Notes | |
---|---|---|---|
1.0.0 |