Abstract - A method for using the popular math packages MATLAB

and Simulink to simulate the behaviour of distributed feedback quan-

tum-well semiconductor laser diodes using the rate equations that de-
scribe them is presented. A large-signal model for threshold determi-
nation and operating point selection is discussed. Small-signal proper-
ties of interest can be investigated using a state-space description of a
linearized version of the rate equations, based on their Jacobian.
The techniques presented here are suitable for investigation of
both analog and digital modulation performance, for example frequency
and pulse response. Large signal gain-saturation and frequency modu-
lation effects are not considered.
I. INTRODUCTION
Modern high-speed data transmission is accomplished
with fiber optic systems, which generally consist of a laser
light source, the fiber itself, and a detector. The most com-
mon light source is the semiconductor laser diode because
of its small size, ruggedness, and manufacturing cost com-
pared to other laser emitters.
A cost-effective method of data transmission is direct
modulation of the laser diode drive current combined with
a photodetector at the other end. This scheme is known as
Intensity Modulation/Direct Detection, or IM/DD. One of
the factors limiting the data rate is the switching speed of
the laser diode. Improving the characterisitics of the diode
is thus an area of extreme interest.
II. A BRIEF REVIEW OF THE RATE EQUATIONS
The behaviour of a semiconductor laser diode is modeled
with a set of three differential equations which describe the
mechanism by which an electrical current causes stimu-
lated photon emission. There are many forms of these equa-
tions, but in this paper we use a modified form of Tuckers
equations [1] employed by Bjerkan et al. [3]:
A Dynamic Simulation Model for Semiconduc-
tor Laser Diodes
Toby Schaer, Robert Rusnov, Stephen Eagle, Jay Jastrebski, Steven Albanese and Xavier Fernando
) 1 (
) (
) ( 1
) ( ] ) ( [ ) ( ) (
n
o
o
a
t N
t S
t S N t N
g
V q
t I
dt
t dN


+

[ ]
) 2 ( ) (
) (
) ( 1
) ( ) ( ) (
t N
t S
t S
t S N t N
g
dt
t dS
n p
o
o


+
+

[ ] ) 3 (
1
) (
2
1 ) (
0
1
1
]
1

p
o
N t N g
dt
t d

) 4 (
2
) (
) (
p
o a
h V t S
t p

Coefficients are defined in Appendix A.

III. THEORY
The large signal single-mode rate equations as given in
Eq.1 and Eq.2 were used to build the initial model. This
models purpose was to identify the various regions of op-
eration of the laser diode based on the parameters found in
the article by Javro and Kang [2]. Building the large signal
model was essential in determining operating points within
the semi linear lasing region. It is these operating points,
once the lasing region was defined above threshold, that
were used to characterize the small signal response. As
depicted in Figure 2, this is in the linear stimulated emis-
sion region.
Figure 2. Full operating region of laser diode (large-signal).
The problem with finding a closed form solution to
the rate equations lies in the fact that the equations are
cross-coupled by the product term containing both the pho-
Equation (1) is the rate equation for carrier density, photon
density is given by Equation (2), and optical phase by Equa-
tion (3). This is a large-signal model that describes laser
diode behaviour in both the spontaneous emission and
stimulated emission region. The nonlinear saturation re-
gion (see Figure 2) due to junction heating is not included
in this model. Additionally, the optical power generated
(also described in [3]) by the diode is:
) 7 ( ) (
1
) (
) (
t S g
t N
dt
t dN
o
n

,
_

[ ] ) 8 ( ) (
) (
) (
o o
N t N g
t S
dt
t dN

,
_

) 9 ( ) (
) (
) (
n
o
t S g
t N
dt
t dS

,
_

[ ] ) 10 (
1
) (
) (
) (
p
o o
N t N g
t S
dt
t dS

,
_

In order to obtain values to simulate the small signal re-

sponse, the non-linear rate equations were linearized by
taking the partial derivative of each time-dependent term
with respect to photon and electron density. It is this new
model that allows the results of the large signal response to
be used in simulating the small signal response.
The process was as follows: Choose an arbitrary operat-
ing point close to the threshold, such that the gain satura-
tion term 1/(1 + S(t)) will approach 1 for values of S(t)
<< 1. This term will have negligible effect on Eq.1 and
Eq.2 and partial derivatives with respect to N(t) and S(t)
can safely be taken.
ton and electron density. For example, the standard Laplace
Transform approach would fail because it is suitable only
for linear, constant coefficient differential equations. Set-
ting each equation to zero and solving the resulting system
merely linearizes the rate equations about the threshold (N
th
).
Simply, this leaves the state-space matrix with linearly in-
separable terms. The solution to this was to linearize the
equations by creating a Jacobian matrix. This new matrix
contains constant coefficient terms when evaluated at spe-
cific operating points. The ultimate goal was to linearize
the rate equations and use them in a state-space model for
small signal simulation.
Now using the Jacobian matrix in a state-space model,
) 11 ( ) (
0
1
t I
V q
S
N
d
dn
n
dt
d
d
dn
n
dt
dn
S
N

1
1
1
1
]
1

+
1
1
1
1
]
1

1
1
1
1
1
1
1
1
]
1

,
_

,
_

,
_

,
_

1
1
1
1
]
1

[ ]
1
]
1

S
N
y 1 0
The final result was the following state-space model,
[ ]
[ ]
) 12 ( ) (
0
1
1
1
t I
V q
S
N
N N g S g
N N g S g
S
N
p
o op o
n
op o
o op o op o
n

1
1
1
1
]
1

+
1
1
1
1
]
1

1
1
1
1
]
1

1
1
1
1
]
1

[ ]
1
]
1

S
N
y 1 0
where S
op
and N
op
indicate the respective steady-state pho-
ton and electron densities for an arbitrary bias current. Fig-
ure 3 illustrates the selection of a particular operating point
within the lasing region of the diode.
Figure 3. Example of operating point.
IV. LARGE SIGNAL SIMULATION MODEL
The full equation model given by Eqs. (1) and (2) was
modeled in Simulink. The method used to transform the
equations into a signal-flow diagram is borrowed from the
field of analog computation the equations were mecha-
nized into a model that gives numerical solutions for arbi-
trary inputs and time scales that are of interest. The model
is shown in Figure 4.
V. LASER OPERATING REGIONS
A ramp input was applied to the model to obtain a graph
similar to Figure 2. The simulation was run over a time
period that was several orders of magnitude longer than
the turn-on transient effects that happen near threshold, to
avoid the oscillatory behaviour in both N(t) and S(t) in this
region. All simulations used a standard 4
th
-order Runge-
Kutta algorithm with a fixed step size of 1ps as an equation
solver. The model parameters are those used in [2], and are
given for reference in Appendix A. A maximum of 20mA
of pumping current was used, and simulations were run
including the e terms, and without. The resulting I(t) vs.
P(t) graphs are shown superimposed in Figure 5.
It is worth noting that the previous assumption of ne-
glecting the effect of (the gain saturation parameter) near
threshold is valid.
Figure 4: Large-Signal Simulink Model
Pumping Current (mA) Electron Density (m
-3
) Photon Density (m
-3
)
10mA 1.9531 x 10
18
1.9221 x 10
13
15mA 1.9571 x 10
18
1.7141 x 10
14
Table I. Operating Points.
VI. SMALL SIGNAL SIMULATIONS
With the operating points from Table I, the Jacobian ma-
trix of Eqs. (7) through (10) could be evaluated and a state-
space model generated for each. A MATLAB script used
for this purpose is shown in Appendix B.
This now permits an analysis of frequency response, for
analog modulation purposes, and pulse response for digital
modulation. Frequency modulation effects such as chirp-
ing were not simulated. Note that the output of the model is
photon density, not optical power, which explains some of
the large gains in the Bode plots!
VII. BODE PLOTS OF FREQUENCY RESPONSE
Both models were analyzed close to their corner frequen-
Figure 5: Laser Response Curves
From Figure 5 it is possible to determine the value of
the threshold current graphically, by extrapolation, as shown
in Figure 6. The threshold current I
th
= 9.32mA. This sets a
lower limit on the selection range of the operating points.
The Simulink model was again employed to find nu-
merical solutions for the operating points, given a constant
pumping current. Two bias points above threshold were
computed:
Figure 6: Laser Threshold Current
cies. From the state-space model it is clear that the laser
diode is a second-order system with no zeros, a fact veri-
fied by the Bode phase plots. Figures 7 & 8 show the re-
sponse of the 10mA system, while 9 & 10 show the 15mA
system.
VIII. PULSE MODULATION
Both models were modulated with a 2ns pulse of 1mA
above the operating point. Figures 11 and 12 show the 10mA
and 15mA models, respectively. The response worsens with
increased bias current, therefore only the 10mA model was
tried with a 1ns pulse of 1mA (see Figure 13).
Figure 10: Phase Response of 15mA System
Figure 11: 10mA Model Response to 2ns, 1mA Pulse
Figure 8: Phase Response of 10mA System
Figure 9: Magnitude Response of 15mA System
Figure 7: Magnitude Response of 10mA Model.
Figure 12: 15mA Model Response to 2ns, 1mA Pulse
IX. SUMMARY
This paper presented a method for simulating a quan-
tum-well distributed feedback (QW-DFB) semiconductor
laser diode based on the rate equations describing it. A large-
signal model was constructed to find the threshold current
of the diode, and to identify a relatively linear region above
threshold for operating point selection. The large signal
model also proved useful for finding numerical solutions
for the steady-state values of electron and photon densities
for currents above threshold.
The small-signal analysis was conducted using a
linearized state-space description, based on taking the
Jacobian of Eqs. (1) and (2), and evaluating at the operat-
ing points previously computed by the large-signal model.
This model was shown to be an aid in determining the fre-
quency response (and thus the bandwidth) of a laser diode
for analog modulation purposes, and modeling the pulse
response (and from this, the maximum bitrate) for digital
modulation techniques.
Example results are given to illustrate the simulation ap-
proach and to verify its validity. The advantages of the ap-
proach presented here are its relative ease of automation,
graphical model construction, and interactive control of
simulation parameters. Detailed studies of specific devices
can thus be carried out efficiently.
Figure 13: 10mA Model Response to 1ns, 1mA Pulse
X. REFERENCES
[1] R.S. Tucker & D.J. Pope, Large Signal Circuit Model for Simulation of
Injection-Laser Modulation Dynamics, IEEE J.Quantum Electron.,vol.QE-
19, pp1179-1183, July 1983.
[2] S.A. Javro & S.M. Kang, Transforming Tuckers Linearized Laser Rate
Equations to a Form that has a Single Solution Regime, Journal of Lightwave
Technology, vol.13, No.9, pp.1899-1904, September 1995.
[3] L. Bjerkan, A. Ryset, L. Hafskjaer, and D.Myhre, Measurement of Laser
Parameters for Simulation of High-Speed Fibre Optic Systems, Journal of
Lightwave Technology, vol.14, No.5, pp.839-850, May 1996.
[4] G. Keiser, Optical Fiber Communications. New York: McGraw-Hill,
2000, pp. 161-193.
APPENDIX A
LIST OF SYMBOLS
Symbol Value Dimension Description
I(t) - A laser current
S(t) - m
-3
photon density
0.44 optical confinement factor
go 3.0 x 10
-6
cm
-3
/s gain slope
N(t) - m
-3
carrier density
NO 1.2 x 10
18
cm
-3
carrier density at transparency
3.4 x 10
-17
cm
3
gain saturation parameter
p 1.0 x 10
-12
s photon lifetime
4.0 x 10
-4
spontaneous emission factor
n 3.0 x 10
-9
s carrier lifetime
Va 9.0 x 10
-11
cm
3
volume of the active region
- phase of the electric field form the laser
- line-width enhancement factor
P(t) W optical power from the laser
q 1.602 x 10
-19
C electronic charge
0.1 total quantum efficiency
h 6.624 x 10
-34
Js Plank's constant
v - s
-1
instantaneous optical frequency deviation
D - s/m
2
fibre dispersion coefficient
L - m fibre length
APPENDIX B
% smallsig.m : Linearized state-space laser diode model. T.Schaer March 2002
clear all; close all; echo off
% Operating Point Values, obtained from Simulink Model
n_10 = 1.9531e18; s_10 = 1.9221e13;
n_15 = 1.9571e18; s_15 = 1.7141e14;
% Other constants (Javro & Kang)
q = 1.6e-19; Va = 9e-11;
tn = 3e-9; tp = 1e-12;
gamma = 0.44; beta = 4e-4;
N0 = 1.2e18; g0 = 3e-6;
eta = 0.1; h = 6.624e-34;
f = 1.997e14;
% Hand-derived Jacobian matrix
A_10 = [-g0*s_10-1/tn -g0*(n_10 - N0); gamma*g0*s_10+gamma*beta/tn gamma*g0*(n_10 -
N0)-1/tp]
A_15 = [-g0*s_15-1/tn -g0*(n_15 - N0); gamma*g0*s_15+gamma*beta/tn gamma*g0*(n_15 -
N0)-1/tp]
% Remaining Matrices
B = [(1/(q*Va)); 0]
C = [0 1]
D = 0
% Check controllability & observability of model
if(rank(ctrb(A_10,B) == 2))
disp(System is controllable for this operating point);
end
if(rank(obsv(A_10,C) == 2))
disp(System is observable for this operating point);
end
% Generate time vector in a useful range (0-14ns)
t = 0:0.0001:14;
t = t*1e-9;
% Matching step input: step time at 1/8 of sim run length, amplitude = 1mA
t_step = round(length(t)/8);
u = [zeros(1,t_step+1) ones(1,(3*t_step)) zeros(1,4*t_step)];
u = u*1e-3;
% Generate state-space objects
sys = ss(A_10,B,C,D); sys_15 = ss(A_15,B,C,D);
% Simulate response to input
y = lsim(sys,u,t); y2 = lsim(sys_15,u,t);
% recall that laser is already on: add steady-state values
y = s_10*ones(length(y),1) + y; y2 = s_15*ones(length(y2),1) + y2;
plot(t,y)
xlabel(Time (nanoseconds)); ylabel(Photon density (1/m^3));
title(Pulse Response: I=1mA from bias point of I = 10mA);
figure
plot(t,y2)
xlabel(Time (nanoseconds));
title(Pulse Response: I=1mA from bias point of I = 15mA);