Modeling of the Semiconductor Laser with Optical Feedback

PoJen Huang, Department of Mathematics and Physics, University of Arizona, Arizona 85721

 

Introduction to the Lang-Kabuysashi equation and its application

 

The goal of this semester is to understand the behavior of the semiconductor laser system with feedback. We use Lang-Kabuysashi equation to describe the system:

…(1)

 

 Where  is the line width enhancement factor,  is the laser frequency without feedback,  is the photon lifetime,  is the carrier life time, and  is the injected current density. The feedback parameter  measure the intensity of the light reflected back in to the laser cavity, and the delay time , which is the round trip time of the light in the external cavity of length .

The modal gain per unit time is given by , where  is the gain constant and  is the carrier density at transparency. The modal gain  includes the linear gain and intensity reduction of the fain due to spatial and spectral hole burning and carrier heating, with  being the gain saturation coefficient. Also the electric field is normalized so that  is the total photon number in the laser wave-guide, where  is the volume of the active region. The parameter  is the lasing threshold current density of a solitary laser and  is the threshold currier density. The typical values for the above parameters are given in the table attached to the end of this report.

        For easier numerical convention, we measure the time in the units of the photon lifetime , and normalize excess carrier number density . Then we can rewrite the Lang Kobayashi equation in the normalized form:

…(1)

 

Where , , , and . We start by looking at some special solutions and analysis their stability, but even before we start, we should first review some of the tools that we will use later in the paper.

 

Stability: The sign of the real part of the eigenvalue determines the stability of the dynamical system. Positive real part indicated an unstable solution, and negative real part indicated stable solution. (one can think the real part of the eigenvalue as the stretch/compress of the system) However, solving the characteristic polynomial to find the eigenvalue is often hard. Since the analytic solution is almost impossible to obtain, one must develop some numerical method to find the real part of the eigenvalue.

 

Numerical Method:     To find the engenvluae in the complex plan is nontrivial. Perhaps the easiest way to find the root for a complex function is to calculate the residue. Recall Cauchy’s formula:

…(1)

Where  is the residue of at and  is the winding number of  around . (since we always take simple loops, the winding number is always one)

If we let , then at the singularity of  is where the root of the characteristic polynomial, ie. The value of the eigenvalue.

Since this method is highly depending on the contour we take, there is a potential danger of applying this method. If the contour contains more then one pole, the integral might be zero & we will miss that root. Hence it is very important to find the right size of the contour. (such that only one pole is allowed in one contour) We will illustrate how such analysis can be done in the simple one-dimensional case.

 

Analytical analysis for the steady state

We want to analysis stability of the solution of Lang-Kabuysashi. This is generally hard, so we start by working the system without feedback.

Case 1: Laser is off with no feedback:

…(3)

…(3)

Stability- negative eigenvalue

Consider a small perturbation from steady-state solutions

…(3)

…(3)

Where

Let

…(3)

 

 

Hence in terms of a matrix:

 

Solve for the eigenvalue we find

 

To yield

 

 

 

For the laser to remain off,  must be less then zero, i.e. for the off state to be stable

 

 

Hence, the threshold value of J

 

This is , the laser will be off, for , any small perturbation will lead to the  filed becoming and staying nonzero, i.e. laser turns on.

 

CASE 2: Constant electric field without feedback:

, and

 

…(3)

 

…(3)

Let’s consider a small perturbation from steady-state solutions

…(3)

 

Hence

…(3)

 

 

 

 

Hence in terms of a matrix:

 

Solving for the eigenvalue we find

 

To yield

 

 

 

        It should be understood that Lang-Kabuysashi has a nontrivial behavior even with the simplest condition. To analysis the system better, we will develop our skill by studying the simplest case.

 

The One Dimensional Case

 

Consider a simple delay differential equation

…(1)

Where  is the change of  at time, is the delay time[1], and  is the value of  at time .

 It is easier to view the problem if we rescale the delay time to be one, i.e. define

…(2)

Then from eq(1) we find

…(3)

Or[2]

…(4)

Like the real Lang-Kobayashi equation, it is also difficult to construct the analytic solution, even in the simplest one-dimensional case. The first few solution are below. (with the unit input step function)

 

Time

Analytic Function

Value at right hand endpoint

Value at the left hand endpoint

A more compact way of express the solution is via Laplace transformation[3]:

N/A

The characteristic polynomial is:

…(5)

Since  is a complex number, we can write it as , . Then eq(3) becomes:

Let’s look at some solutions for different value of .

 

Case I

 

Figures: The top left is the solution of the one-dimensional differential equation with , integrate from one to twenty. Top right is the absolute log plot of the top left figure. The slope is determined by the envelop (the maximal value of each hump). Bottom right is the zero contour of the real and imagery part of the characteristic function. Red indicates the real part and blue indicated the imagery part. Notice that the imagery part (the vertical direction) crossing of two contours are evenly spaced. (this is important later when we try to drive the contour size for auto root finder) Bottom left shows the detail around the origin. They cross on the left plan, i.e. negative real part of eignevalue)

 

 

Case II

 

Figures: The top left is the solution of the one-dimensional differential equation with , integrate from one to twenty. Top right is the absolute log plot of the top left figure. The slope is determined by the envelop (the maximal value of each hump). Bottom right is the zero contour of the real and imagery part of the characteristic function. Red indicates the real part and blue indicated the imagery part. Notice that the imagery part (the vertical direction) crossing of two contours are evenly spaced. Bottom left shows the detail around the origin. The crossing is on the imagery axis, i.e. the real part of eignevalue is identically zero)

 

 

Case III

 

Figures: The top left is the solution of the one-dimensional differential equation with , integrate from one to twenty. Top right is the absolute log plot of the top left figure. The slope is determined by the envelop (the maximal value of each hump). Bottom right is the zero contour of the real and imagery part of the characteristic function. Red indicates the real part and blue indicated the imagery part. Notice that the imagery part (the vertical direction) crossing of two contours are evenly spaced. Bottom left shows the detail around the origin. They cross on the right plan, i.e. positive real part of eignevalue)

 

Notice the following two things, (base on observing the three cases):

1.)        That the slope of the log plots is the same as the maximum of the real part of the eigenvalue.

2.)        The imagery part of the eigenvalue is evenly spread, and the distance in between is about .

 

Observation one is not suppressing, since the ratio of the stretching/ compress is determined by the real part of the eigenvalue. Recall the ratio of stretching/ compressing is given by the formula:

        The final solution can be written as:

the  is the envelop of the solution and  is the oscillated solution.

 

PFPFPFPFPF

 PF

 

Adams-Bashford-Moulton (ABM) PC Integration Method

We can rewrite the differential equation in terms of discrete map of dimension  by introducing the step size to be

 

 

Hence the system becomes:

 

 

 

Since the L-K is extremely depend on the integration method we choice. We decide to use the fourth-order Adams-Bashford-Moulton (ABM) PC method.

 

Predictor:

 

Corrector:

 

Where

 

Logarithm derivative

Phase Plot:

, and

 

Let’s rewrite the complex electric field,

 

Since the constant phase shift does not affect the dynamics of the system.

 

Consider the solution has the form:

 

The arctangent only give the value up to pi. (depend on the slice you are on), however the dynamics of the system changes as time varies. We need a way to describe

 

Take the ratio we yield

 

Provided that the electric field is non-zero. (if the field is zero, then the phase is undefined)

 

Integrate both sides we have

 

Then do a little shift on the argument to get the interval we are interested.

 

PHASE DELAY

 

Bifurcation diagram and Sisyphus effect

N/A



[1] For simplicity, the delay time in this paper would be a fix positive number.

[2] By rename the delay time  to be .

[3] For the people who would like to know the detail, there is a great book by **** that has all the calculation in it.