# 3-Level Rate Equation modelling using Matlab

• xzyan
In summary, the problem is to use Matlab to model the 3-level rate equation of a solid-state laser. The simplified energy levels and transitions of the laser are shown in the picture. The problem is that the MATLAB keeps running for a long time and no output is obtained. The Attempt at a Solution is to write down the rate equations according to the energy levels and transitions. The main idea of this modelling is to observe the change of the population of the state in each level and the change of the intensity of both the pump and signal lasers.
xzyan
Hello everyone, I'm new here and new to Matlab. I hope I do place the post in the right place.

First of all, thanks for viewing my post. Please bear with my English. I tried my best to explain everything clear. If you have any question about anything I wrote or spot any fault, please tell me. Many thanks!

1. The problem statement
The problem is to use Matlab to modelling the 3-level rate equation of a solid-state laser.

The picture shows the simplified energy levels and transitions of the laser.

S0.2 is the absorption from ground level (L0) to upper lasing level (L2).

A2.0 is the spontaneous emission from upper lasing level (L2) to ground level (L0).

A2.1 is the spontaneous emission from upper lasing level (L2) to lower lasing level (L1).

SE2.1 is the stimulated emission from upper lasing level (L2) to lower lasing level (L1).

W1.0 is the fast non-radiative decay from lower lasing level (L1) to ground level (L0).

In this problem, the absorption S0.2 is triggered by pump laser source with intensity of Ip, stimulated emission SE2.1 is triggered by signal laser source with intensity of Is.

## Homework Equations

I did write a MATLAB function with a main script to do the modelling. However, the MATLAB keeps running for a long time and no output is obtained. (I've tried run it on the server for 3 days) I did try a simplified 2 level system with only one beam (i.e. the pump only), and I got the results in about half an hour.

I was thinking, is it because of the solving function I was using is not suitable for this problem? (I'm using ODE45)

## The Attempt at a Solution

I've written down the rate equations according to the energy levels and transitions like follows:

where n2, n1 and n0 are the populations of the state in upper lasing level (L2), lower lasing level (L1) and ground level (L0), respectively. σ02 and σ21 are the absorption and stimulation emission cross section of S0.2 and SE2.1 transitions. A20 and A21 is the Einstein coefficient of transition A2.1 and A2.0. W10 is the rate of the fast transition W1.0. h is the Planck's constant, vp and vs are the frequencies of the pump laser and signal laser. Ip and Is are the intensities of the pump and signal laser.

The main idea of this modelling is to observe the change of the population of the state in each level and the change of the intensity of both the pump and signal lasers.

I have written a Matlab code as below:

Function

function dy = rate_eq( t,y )
%Simplfied rate equation for Yb ions. Detailed equation is in Mini-thesis
dy = zeros(5,1);
h = 6.62606957E-34;
% Planck constant (m^2kg/s)
c = 299792458;
% Speed of light (m/s)
lambdaS = 1030;
% Wavelength of signal light (nm)
lambdaP = 975;
% Wavelength of pump light (nm)
vS = c/lambdaS;
% Frequency of signal (Hz)
vP = c/lambdaP;
% Frequency of pump (Hz)
A21 = 1/(600E-6);
% Spontaneous emission transition time from L2 -> L1 of 600 us
A20 = 1/(600E-6);
% Spontaneous emission transition time from L2 -> L1 of 600 us
W10 = 1/(500E-12);
% Fast decay transition time from L1 -> L0 of 500 ps
sigma21 = 0.5E-20;
% Stimulation emission cross section from L2 -> L1 of 0.5E-20 cm^2
sigma02 = 2.5E-20;
% Absorption cross section from L0 -> L2 of 2.5E-20 cm^2
dy(1) = sigma02*y(4)*y(3)/(h*vP)-A20*y(1)-A21*y(1)-sigma21*y(5)*y(1)/(h*vS)
% Change of population of upper lasing level (L2)
dy(2) = A21*y(1)+sigma21*y(5)*y(1)/(h*vS)-W10*y(2)
% Change of population of lower lasing level (L2)
dy(3) = -sigma02*y(4)*y(3)/(h*vP)+A20*y(1)+W10*y(2)
% Change of population of ground lasing level (L2)
dy(4) = (-sigma02*y(4)*y(3)/(h*vP)+A20*y(1))*(h*vP)
% Change of intensity of pump beam
dy(5) = (sigma21*y(5)*y(1)/(h*vS)+A21*y(1))*(h*vS)
% Change of intensity of signal beam
end

Main Script

TSPAN = [0 1];
Y0 =[0 0 8.75E20 1 0.2];
%y(1) = population of Upper Lasing level
%y(2) = population of Lower Lasing level
%y(3) = population of Ground Level
%y(4) = Optical intensity of Pump W/cm^2
%y(5) = Optical intensity of Signal W/cm^2

[T,Y] = ode45(@rate_eq,TSPAN,Y0);

subplot(3,1,1)
plot(T,Y(:,3),T,Y(:,2),T,Y(:,1))
legend('Ground level','Lower lasing level','Upper lasing level')
title('Population different levels')

subplot(3,1,2)
plot(T ,Y(:,4))
title('Pump Intensiry') % Pump intensity

subplot(3,1,3)
plot(T ,Y(:,5))
title('Singal Intensity') % Signal intensity

Thanks again for your patience and kindly help!

If I calculate the dy at initialization, I get

dy(1) 1.7123e+026
dy(2) 4.44032e+023
dy(3) -1.71674e+026
dy(4) -0.918435
dy(5) 0.836369

which is somewhat stiff (big ratios between derivatives).
Matlab will probably work in double precision (someone?) which means 16 digits. Any step it takes will be too big, and it can't reduce the step size enough to satisfy error criteria. Is my guess (don't know matlab).

Perhaps you want to do some scaling to get a numericallymore reasonable problem

BvU said:
If I calculate the dy at initialization, I get

dy(1) 1.7123e+026
dy(2) 4.44032e+023
dy(3) -1.71674e+026
dy(4) -0.918435
dy(5) 0.836369

which is somewhat stiff (big ratios between derivatives).
Matlab will probably work in double precision (someone?) which means 16 digits. Any step it takes will be too big, and it can't reduce the step size enough to satisfy error criteria. Is my guess (don't know matlab).

Perhaps you want to do some scaling to get a numericallymore reasonable problem

Thanks a lot for your reply. That's a good idea. I'll try to realize it. Many thanks!

Some scaling already occurs if you fix the units. v is now in GHz, not in Hz. In your program (and initally also in mine...)

Last edited:
Oh, and welcome to PF !

Thanks! I realized that after you mention about the scaling. Tried to modified it in the post to avoid misleading but it seems I cannot modify it anymore.

BvU said:
Some scaling already occurs if you fix the units. v is now in GHz, not in Hz. In your program (and initally also in mine...)

I've changed the 4th and 5th equations from intensity of the beam (W/cm2) to number of photons, which scaling the dy(4) and dy(5) to around E19 rather than 0.2 and 1. However, the MATLAB still keeps calculating and calculating.

Later, I tried another ode function called 'ode15s', which used for solving stiff differential equation, and I got the calculation done within a sec. Hence the problem is more likely to be an unsuitable choice of solver. Thanks for the clue you gave in previous reply 'somewhat stiff'!

Problem solved. Although there are still some problem of the model to be investigated as the current one returns an unusual result.

Thanks a lot for your help BvU! And of course everyone who spent time viewing my post and tried to help!

Unusual result sounds ominous... Your intial condition is pretty far from the working point. I played with the D.E. too and found little going on except filling up the upper level (and then running into integrator problems too).

Is the thing solvable for steady state? (all d/dt = 0)

Yes it is or should be solvable at steady state. The rate equations are derived by myself, that's why I thought the problem now is not how to solve them but how to correct them (if I was wrong). According to the unexpected result, I think I was wrong at some where in the equation deriving (probably a factor or missing a part).

From physical point of view, the process is: the dopants in the sample absorb the energy from pumping Ip, promoting themselves from ground level to upper lasing level, and then drop to different levels through different kinds of transitions.

The aim of the modeling is to see the growing/decaying in the population in different levels and beam intensities over time.

Last edited:
With a reasonable pumping level, I expect to see something similar like this:

where the green curve is the population of upper lasing level, blue curve is the population of ground state, orange curve is the population of lower lasing level.

At the moment, the model is not fully completed yet. I need to find out a way to continuously pumping the system for a while rather than just giving an initial pumping energy.

If I want to get the population inversion of upper lasing level and ground level (i.e. n2 > n0) like showed in the figure, I have to set the initial value of y(4) to 30000 (in the unit of W/cm2), which equals 30,000 J in energy in this modeling.

Hi xzyan,
Did you manage to get what you want?
If I am not wrong, there is 'area' missing in equations 5 and 6 if one expect change in intensity. As it is now, the unit will be energy per unit time- basically the power. But in other equations, Is means intensity, of course not the power.
It will be nice if you update this thread so that people like me may benefit or let's say save some time.
Thanks.

## 1. What is 3-Level Rate Equation modelling and how is it used in science?

3-Level Rate Equation modelling is a mathematical approach used to describe the kinetics of a system with three distinct energy levels. It is commonly used in various areas of science such as chemistry, physics, and biology to understand the behavior and interactions of molecules, particles, and organisms.

## 2. What is Matlab and why is it used for 3-Level Rate Equation modelling?

Matlab is a programming language and software environment commonly used in scientific research and data analysis. It is particularly useful for 3-Level Rate Equation modelling because it offers a wide range of tools and functions for solving complex mathematical equations and simulating systems with multiple variables.

## 3. How does 3-Level Rate Equation modelling using Matlab differ from other modelling methods?

3-Level Rate Equation modelling using Matlab is a computer-based approach that allows for the analysis of complex systems with multiple variables and interactions. It differs from other modelling methods, such as analytical or numerical approaches, by providing a more efficient and accurate way to solve complex equations and simulate real-world systems.

## 4. What are some applications of 3-Level Rate Equation modelling using Matlab?

3-Level Rate Equation modelling using Matlab has a wide range of applications in various fields of science, including chemical kinetics, quantum mechanics, population dynamics, and ecological modeling. It is also commonly used in industries such as pharmaceuticals, materials science, and environmental science to study the behavior of complex systems and predict their future outcomes.

## 5. Can beginners use Matlab for 3-Level Rate Equation modelling?

Yes, beginners can use Matlab for 3-Level Rate Equation modelling with some basic knowledge and understanding of programming and mathematics. There are also many online tutorials and resources available to help beginners learn how to use Matlab for scientific modelling and analysis.

Replies
4
Views
1K
Replies
5
Views
3K
Replies
1
Views
2K
• Earth Sciences
Replies
10
Views
4K
• MATLAB, Maple, Mathematica, LaTeX
Replies
1
Views
5K
• Differential Equations
Replies
1
Views
2K
Replies
1
Views
5K
• Engineering and Comp Sci Homework Help
Replies
1
Views
2K
• Quantum Physics
Replies
1
Views
1K
• MATLAB, Maple, Mathematica, LaTeX
Replies
0
Views
2K