RP Photonics logo
RP Photonics
Technical consulting services on lasers, nonlinear optics, fiber optics etc.
Profit from the knowledge and experience of a top expert!
Powerful simulation and design software.
Make computer models in order to get a comprehensive understanding of your devices!
Success comes from understanding – be it in science or in industrial development.
The famous Encyclopedia of Laser Physics and Technology – available online for free!
The ideal place for finding suppliers for many photonics products.
Advertisers: Make sure to have your products displayed here!
… combined with a great Buyer's Guide!
VLib part of the

The Photonics Spotlight

Rate Equations – An Example for Stiff Sets of Differential Equations

Dr. Rüdiger Paschotta

Ref.: encyclopedia articles on rate equation modeling

erbium energy levels and transitions

Figure 1: Level system of Er3+ ions.

The temporal evolution of the populations of electronic levels in some laser gain medium (containing e.g. rare earth ions) can be calculated with rate equations. These are systems of ordinary differential equations, with coupling terms resulting from processes like spontaneous emission, non-radiative transitions, stimulated emission and energy transfers. In principle, it is simple to numerically calculate the temporal evolution even for relatively complicated rate equation systems. This can also be used for approximately calculating the steady-state population values for given optical intensities e.g. of a pump wave and a signal wave.

However, there is a frequently encountered problem: in many situations of interest, the rate equations represent a so-called “stiff” set of differential equations. This happens when the transitions between different levels are caused by physical processes which operate on very different time scales. For example, consider the level system of erbium (as used in erbium-doped fiber amplifiers), which is shown in Figure 1. As explained in the encyclopedia article on rate equation modeling, the rate equations for a situation with some given pump and signal intensities read as follows:

transition rates in erbium ions

The simulated temporal evolution for some set of parameters is shown in Figure 2.

temporal evolution of erbium populations

It is an essential detail that the coefficient A32, describing non-radiative multi-phonon transitions from the pump level to the upper laser level, is usually much larger than the other coefficients, including the terms with the pumping rates. For that reason, ions pumped into level 3 will very quickly be transferred into level 2, so that the population of level 3 itself remains very small.

The curves shown look very smooth, so that one may expect it to be very easy and efficient to integrate the rate equations numerically with a step size of say 50 μs. If you try that, however, you will find that numerical results are absolutely wrong, wildly oscillating with quickly increasing amplitude. Good results are obtained only once you reduce the temporal step size to the order of 1/A32, which is 1 μs in this example. With the Euler method, where one step simply means to add the temporal derivatives times the temporal step size to the current populations, the instability starts for step sizes just above 2 μs. With the more refined Runge–Kutta algorithm, you can go only slightly further to a step size of about 3.7 μs. Of course, that behavior would be even far worse if A32 were substantially larger.

The observed behavior can be understood by first considering a simple exponential decay (with a single time constant). It is easy to see that the Euler method will deliver oscillating results when the temporal step size is larger than the time constant of the system, and it will even be unstable when the step size is at least twice the time constant. Figure 2 illustrates this: the lowest curve is for a step size of half the time constant (resulting in moderate errors), and the higher curves are for integer multiples of half the time constant.

Euler propagation

Figure 2: Behavior of the Euler method, applied to a simple exponential decay with different numerical step sizes. The correct solution is always shown as a thin blue curve. When the numerical step size is larger than twice the time constant (which is the case for the highest curve), the propagation becomes totally unstable, with exponentially increasing errors.

Now this problem doesn't go away if the system is more complex and has some slower processes in addition. Even when the exact value of the fast time constant of the system has nearly no effect on the longer-term evolution of the system, the corresponding process makes the numerical solution unstable when the step size is too large.

Examining the concrete rate equation system (see above) can also be instructive. When we start with all erbium ions being in the ground state (level 1), a first Euler step will pump some population from level 1 to 3. The calculated population of level 3 is then too high, as the decay from level 3 to 2 is not yet considered. In the next step, there will be a strong transfer from level 3 to 2, and if the step size is larger than twice the corresponding time constant, the population of level 3 will “overshoot” and get substantially negative (which is physically unreasonable, of course). In the third step, we obtain a transfer in the opposite direction, making the deviation from the correct solution even larger.

We can learn several things:

I actually got inspired to write about this by my work on the substantially extended version 2 of the software RP Fiber Power, which calculates the optical powers in fiber amplifiers and fiber lasers. While the previous version could be used only for simple laser-active ions with only a single metastable level, the new version allows the user to define complicated level systems process by process, including energy transfers even between different types of ions (such as Er3+ and Yb3+). This software (which is basically finished and currently undergoes intensive testing) must efficiently calculate the steady-state populations for arbitrary cases, including those where the rate equations are extremely stiff, and also possibly nonlinear. Therefore, I had to implement the multidimensional Newton–Raphson method, as mentioned above.

This article is a posting of the Photonics Spotlight, authored by Dr. Rüdiger Paschotta. You may link to this page, because its location is permanent. See also the Encyclopedia of Laser Physics and Technology.

Note that you can also receive the articles in the form of a newsletter or with an RSS feed.

If you like this article, share it with your friends and colleagues, e.g. via social media: