# Lyapunov Potential Description for Laser Dynamics

## I Introduction

Even for non-mechanical systems, it is occasionally possible to construct a function (called Lyapunov function or Lyapunov potential) that decreases along trajectories . The usefulness of Lyapunov functions lies on the fact that they allow an easy determination of the fixed points of a dynamical (deterministic) system as the extrema of the Lyapunov function as well as determining the stability of those fixed points. In some cases, the existence of a Lyapunov potential allows an intuitive understanding of the transient and stationary trajectories as movements of test particles in the potential landscape. In the case of non–deterministic dynamics, i.e. in the presence of noise terms, and under some general conditions, the stationary probability distribution can also be governed by the Lyapunov potential and averages can be performed with respect to a known probability density function. The aim of this work is to construct Lyapunov potentials for some laser systems. We start, then, by briefly reviewing the main features of the laser as a dynamical system.

A laser has three basic ingredients: i) a gain medium capable of amplifying the electromagnetic radiation propagating inside the cavity, ii) an optical cavity that provides the necessary feedback, and iii) a pumping mechanism. A complete understanding of laser dynamics is based on a fully quantum-mechanical description of matter-radiation interaction within the laser cavity. However, the laser is a system where the number of photons is much larger than one, thus allowing a semiclassical treatment of the electromagnetic field inside the cavity through the Maxwell equations. This fact was introduced in the semiclassical laser theory, developed by Lamb  and independently by Haken . This model for laser dynamics was constructed from the Maxwell-Bloch equations for a single-mode field interacting with a two-level medium. The semiclassical laser theory ignores the quantum-mechanical nature of the electromagnetic field, and the amplifying medium is modeled quantum-mechanically as a collection of two-level atoms through the Bloch equations. A simpler description can be obtained by deriving rate equations for the temporal change of the electric field (or photons number) inside the cavity and the population inversion (carriers number in the case of semiconductor lasers) . Rate equations, with stochastic terms accounting for spontaneous emission noise, have been extensively used for semiconductor lasers.

Different types of lasers can be classified according to the decay rate of the photons, carriers and material polarization. Arecchi et al.  were the first to use a classification scheme: class C lasers have all the decay rates of the same order, and therefore a set of three nonlinear differential equations is required for a satisfactory description of the electric field, the population inversion and the material polarization. For class B lasers, the polarization decays towards the steady state much faster than the other two variables, and it can be adiabatically eliminated. Class B lasers, of which semiconductor lasers  are an example, are then described by just two rate equations for the atomic population inversion (or carriers number) and the electric field. Another examples of class B lasers are CO<sub>2</sub> lasers and solid state lasers . Finally, in class A lasers population inversion and material polarization decay much faster than the electric field. Both material variables can be adiabatically eliminated, and the equation for the electric field is enough to describe the dynamical evolution of the system. Some properties of class A lasers, like a dye laser, are studied in . In this paper we interpret the dynamics of both class A and class B lasers by using a Lyapunov potential.

The paper is organized as follows. In Section II we present a brief review of the relation of Lyapunov potentials to the dynamical equations and the splitting of those into conservative and dissipative parts. We consider the example of class A lasers. In this case, the Lyapunov potential gives an intuitive understanding of the dynamics observed in the numerical simulations. In the presence of noise, the probability density function obtained from the potential allows the calculation of stationary mean values of interest as, for example, the mean value of the number of photons. We will show that the mean value of the phase of the electric field in the steady state varies linearly with time only when noise is present, in a phenomenon reminiscent of the noise–sustained flows. In Section III, the dynamics of rate equations for class B lasers is presented in terms of the intensity and the carriers number (we will restrict ourselves to the semiconductor laser). In this case we have found a potential which helps to analyze the corresponding dynamics in the absence of noise. By using the conservative part of the equations, one can obtain an expression for the period of the oscillations in the transient regime following the laser switch-on. This expression extends the one obtained in a simpler case by an identification of the laser dynamics with a Toda oscillator in . Here, we have added in the expression for the period the corresponding modifications for the gain saturation term and spontaneous emission noise. Finally, in section IV, we summarize the main results obtained.

## II Potentials and Lyapunov Functions: Class A Lasers

The evolution of a system (dynamical flow) can be classified into different categories according to the relation of the Lyapunov potential to the actual equations of motion . We first consider a deterministic dynamical flow in which the real variables $`(x_1,\mathrm{},x_N)𝐱`$ satisfy the general evolution equations:

$$\frac{dx_i}{dt}=f_i(𝐱),i=1,\mathrm{},N$$
(1)   
In the so–called potential flow, there exists a non–constant function $`V(𝐱)`$ (the potential) in terms of which the above equations can be written as:

$`{\displaystyle \frac{dx_i}{dt}}={\displaystyle \underset{j=1}{\overset{N}{}}}S_{ij}{\displaystyle \frac{V}{x_j}}+v_i`$ (2)   
where $`S(𝐱)`$ is a symmetric and positive definite matrix, and $`v_i(𝐱)`$ satisfy the orthogonality condition:

$$\underset{i=1}{\overset{N}{}}v_i\frac{V}{x_i}=0.$$
(3)   
A non-potential flow, on the other hand, is one for which the splitting (2), satisfying (3), admits only the trivial solution $`V(𝐱)=`$ constant, $`v_i(𝐱)=f_i(𝐱)`$.

Since the above (sufficient) conditions for a potential flow lead to $`dV/dt0`$, one concludes that $`V(𝐱)`$ (when it satisfies the additional condition of being bounded from below) is a Lyapunov potential for the dynamical system. In this case, one can get an intuitive understanding of the dynamics: the fixed points are given by the extrema of $`V(𝐱)`$ and the trajectories relax asymptotically towards the surface of minima of $`V(𝐱)`$. This decay is produced by the only effect of the terms containing the matrix $`S`$ in Eq. (2), since the dynamics induced by $`v_i`$ conserves the potential, and $`v_i(𝐱)`$ represents the residual dynamics on this minima surface. A particular case of potential flow is given when $`v_i(𝐱)`$ can also be derived from the potential, namely:

$$\frac{dx_i}{dt}=\underset{j=1}{\overset{N}{}}D_{ij}\frac{V}{x_j}$$
(4)   
where the matrix $`D(𝐱)=S(𝐱)+A(𝐱)`$, splits into a positive definite symmetric matrix, $`S`$, and an antisymmetric one, $`A`$. In this case, the residual dynamics also ceases after the surface of minima of $`V(𝐱)`$ has been reached.

We now describe the effect of noise on the dynamics of the above systems. The stochastic equations (considered in the Itô sense) are:

$$\frac{dx_i}{dt}=f_i(𝐱)+\underset{j=1}{\overset{N}{}}g_{ij}(𝐱)\xi _j(t)$$
(5)   
where $`g_{ij}(𝐱)`$ are given functions and $`\xi _j(t)`$ are white noise: Gaussian random processes of zero mean and correlations:

$$\xi _i(t)\xi _j(t^{})=2ϵ\delta _{ij}\delta (tt^{})$$
(6)   
$`ϵ`$ is the intensity of the noise.

In the presence of noise terms, it is not adequate to talk about fixed points of the dynamics, but rather consider instead the maxima of the probability density function $`P(𝐱,t)`$, which satisfies the multivariate Fokker-Planck equation  whose general solution is unknown. When the deterministic part of (5) is a potential flow, however, a closed form for the stationary distribution $`P_{st}(𝐱)`$ can be given in terms of the potential $`V(𝐱)`$ if the following (sufficient) conditions are satisfied:

1. The fluctuation–dissipation condition, relating the symmetric matrix $`S`$ to the noise matrix $`g`$:

$$S_{ij}=\underset{k=1}{\overset{N}{}}g_{ik}g_{jk},S=gg^T$$
(7)
2. $`S_{ij}`$ satisfies:

$$\underset{j=1}{\overset{N}{}}\frac{S_{ij}}{x_j}=0,i$$
(8)   
This condition is satisfied, for instance, for a constant matrix $`S`$.
3. $`v_i`$ is divergence free:

$$\underset{i=1}{\overset{N}{}}\frac{v_i}{x_i}=0$$
(9)

this third condition is automatically satisfied for potential flows of the form (4) with a constant matrix A.

Under those circumstances, the stationary probability density function is:

$$P_{st}(𝐱)=Z^1\mathrm{exp}\left(\frac{V(𝐱)}{ϵ}\right)$$
(10)   
where $`Z`$ is a normalization constant. Graham  has shown that if conditions 2 and 3 are not satisfied, then the above expression for $`P_{st}(𝐱)`$ is still valid in the limit $`ϵ0`$.

As an example of the use of Lyapunov potentials in a dynamical system, we consider class A lasers  whose dynamics can be described in terms of the slowly varying complex amplitude $`E`$ of the electric field:

$$\dot{E}=(1+i\alpha )\left(\frac{\mathrm{\Gamma }}{1+\beta |E|^2}\kappa \right)E+\zeta (t)$$
(11)   
where $`\alpha `$, $`\beta `$, $`\mathrm{\Gamma }`$ and $`\kappa `$ are real parameters. $`\kappa `$ is the cavity decay rate; $`\mathrm{\Gamma }`$ the gain parameter; $`\beta `$ the saturation-intensity parameter and $`\alpha `$ is the detuning parameter. Another widely used model expands the non-linear term to give a cubic dependence on the field (third order Lamb theory ), but this is not necessary here. Eq. (11) is written in a reference frame in which the frequency of the on steady state is zero . $`\zeta (t)`$ is a complex Langevin source term accounting for the stochastic nature of spontaneous emission. It is taken as a Gaussian white noise of zero mean and correlations:

$$\zeta (t)\zeta ^{}(t^{})=4\mathrm{\Delta }\delta (tt^{}),$$
(12)   
where $`\mathrm{\Delta }`$ measures the strength of the noise.

By writing the complex variable $`E`$ as $`E=x_1+ix_2`$ and introducing a new dimensionless time such that $`t\kappa t`$, the evolution equations become:

$`\dot{x}_1`$ $`=`$ $`\left({\displaystyle \frac{a}{b+x_1^2+x_2^2}}1\right)(x_1\alpha x_2)+\xi _1(t)`$ (13)   
$`\dot{x}_2`$ $`=`$ $`\left({\displaystyle \frac{a}{b+x_1^2+x_2^2}}1\right)(\alpha x_1+x_2)+\xi _2(t)`$ (14)   
Where $`a=\mathrm{\Gamma }/(\kappa \beta )`$ and $`b=1/\beta `$. $`\xi _1(t)`$ and $`\xi _2(t)`$ are white noise terms with zero mean and correlations given by equation (6) with $`ϵ=\mathrm{\Delta }/\kappa `$.

In the deterministic case ($`ϵ=0`$), these dynamical equations constitute a potential flow of the form (4) where the potential $`V(𝐱)`$ is 

$$V(x_1,x_2)=\frac{1}{2}[x_1^2+x_2^2a\mathrm{ln}(b+x_1^2+x_2^2)]$$
(15)   
and the matrix $`D(𝐱)`$ (split into symmetric and antisymmetric parts) is:

$$D=S+A=\left(\begin{array}{cc}1& 0\\ 0& 1\end{array}\right)+\left(\begin{array}{cc}0& \alpha \\ \alpha & 0\end{array}\right).$$
(16)   
A simpler expression for the potential is given in  and  valid for the case in which the gain term is expanded in Taylor series.

According to our discussion above, the fixed points of the deterministic dynamics are the extrema of the potential $`V(𝐱)`$: for $`a>b`$ there is a maximum at $`(x_1,x_2)=0`$ (corresponding to the laser in the off state) and a line of minima given by $`x_1^2+x_2^2=ab`$ (see Fig. 1). The asymptotic stable situation, then, is that the laser switches to the on state with an intensity $`I|E|^2=x_1^2+x_2^2=ab`$. For $`a<b`$ the only stable fixed point is the off state $`I=0`$.

In the transient dynamics, the symmetric matrix $`S`$ is responsible for driving the system towards the line of minima of $`V`$ following the lines of maximum slope of $`V`$. The antisymmetric part $`A`$ (which is proportional to $`\alpha `$) induces a movement orthogonal to the direction of maximum variation of $`V(𝐱)`$. The combined effects of $`S`$ and $`A`$ produce a spiraling trajectory in the $`(x_1,x_2)`$ plane. The angular velocity of this spiral movement is proportional to $`\alpha `$. Asymptotically, the system tends to one of the minima in the line $`I=ab`$, the exact location depending on the initial conditions. The potential decreases in time until it arrives at its minimum value: $`V(x_1^2+x_2^2=ab)=\frac{1}{2}(a\mathrm{ln}(a)a+b)`$.

In the presence of moderate levels of noise, $`ϵ0`$, the qualitative features of the transient dynamics remain the same as in the deterministic case. The most important differences appear near the stationary situation. As the final value of the intensity is approached and for $`\alpha 0`$, the phase rotation slows down and the mean value of the phase $`\varphi `$ of the electric field $`E`$ changes linearly with time also in the steady state, see Fig. 2. For $`\alpha =0`$ there is only the ordinary phase diffusion around the circumference $`x_1^2+x_2^2=ab`$ that represents the set of all possible deterministic equilibrium states . Therefore, for $`\alpha 0`$ the real and imaginary parts of $`E`$ oscillate not only in the transient dynamics but also in the steady state, and while the frequency of the oscillations still depends on $`\alpha `$ (as well as $`ϵ`$), their amplitude depends on the noise strength $`ϵ`$.

We can understand these aforementioned features of the noisy dynamics using the deterministic Lyapunov potential $`V(x_1,x_2)`$. Since conditions 1-3 above are satisfied, the stationary probability distribution is given by (10) with $`V(x_1,x_2)`$ given by (15). By changing variables to intensity and phase, we find that the probability density functions for $`I`$ and $`\varphi `$ are independent functions, $`P_{st}(\varphi )=1/(2\pi )`$ is a constant and,

$$P_{st}(I)=Z^1\mathrm{e}^{I/(2ϵ)}(b+I)^{a/(2ϵ)}$$
(17)   
where the normalization constant is $`Z=(2ϵ)^{\frac{a}{2ϵ}+1}\mathrm{e}^{\frac{b}{2ϵ}}\mathrm{\Gamma }(\frac{a}{2ϵ}+1,\frac{b}{2ϵ})`$ and $`\mathrm{\Gamma }(x,y)`$ is the incomplete Gamma function. From this expression, we see that, independently of the value for $`ϵ`$, $`P_{st}(I)`$ has its maxima at the deterministic stationary value $`I_m=ab`$. Starting from a given initial condition corresponding, for instance, to the laser in the off state, the intensity fluctuates around a mean value that increases monotonically with time. In the stationary state, the intensity fluctuates around the deterministic value $`I_m=ab`$ but, since the distribution (17) is not symmetric around $`I_m`$, the mean value $`I_{st}`$ is larger than the deterministic value. By using (17) one can easily find that

$$I_{st}=(ab)+2ϵ\left[1+\frac{\mathrm{exp}(b/2ϵ)(b/2ϵ)^{\frac{a}{2ϵ}+1}}{\mathrm{\Gamma }(\frac{a}{2ϵ}+1,\frac{b}{2ϵ})}\right]$$
(18)   
An expression for the mean value of the intensity in the steady state was also given in  in the simpler case of an expansion of the saturation-term parameter in the dynamical equations.

As mentioned before, in the steady state of the stochastic dynamics, the phase $`\varphi `$ of the electric field fluctuates around a mean value that changes linearly with time. Of course, since any value of $`\varphi `$ can be mapped into the interval $`[0,2\pi )`$, this is not inconsistent with the fact that the stationary distribution for $`\varphi `$ is a uniform one. We can easily understand the origin of this noise sustained flow : the rotation inducing terms, those proportional to $`\alpha `$ in the equations of motion, are zero at the line of minima of the potential $`V`$ and, hence, do not act in the steady deterministic state. Fluctuations allow the system to explore regions of the configuration space $`(x_1,x_2)`$ where the potential is not at its minimum value. Since, according to Eq. (18), the mean value of $`I`$ is not at the minimum of the potential, there is, on average, a non-zero contribution of the rotation terms producing the phase drift observed.

The rotation speed can be calculated by writing the evolution equation for the phase of the electric field as:

$$\dot{\varphi }=\left(\frac{a}{b+I}1\right)\alpha +\frac{1}{\sqrt{I}}\xi (t)$$
(19)   
where $`\xi (t)`$ is a white noise term with zero mean value and correlations given by (6). By taking the average value and using the rules of the Itô calculus, one arrives at:

$$\dot{\varphi }=\alpha \frac{a}{b+I}1$$
(20)   
and, by using the distribution (17), one obtains the stochastic frequency shift:

$$\dot{\varphi }_{st}=\alpha \frac{\mathrm{exp}(b/2ϵ)(b/2ϵ)^{\frac{a}{2ϵ}}}{\mathrm{\Gamma }(\frac{a}{2ϵ}+1,\frac{b}{2ϵ})}$$
(21)   
Notice that this average rotation speed is zero in the case of no detuning ($`\alpha =0`$) or for the deterministic dynamics ($`ϵ=0`$) and that, due to the minus sign, the rotation speed is opposite to that of the deterministic transient dynamics when starting from the off state. These results are in excellent agreement with numerical simulations of the rate equations in the presence of noise (see Fig. 2).

## III Class B lasers

The dynamics of a typical class B laser, for instance a single mode semiconductor laser, can be described in terms of two evolution equations, one for the slowly-varying complex amplitude $`E`$ of the electric field inside the laser cavity and the other for the carriers number $`N`$ (or electron-hole pairs) . These equations include noise terms accounting for the stochastic nature of spontaneous emission and random non-radiative carrier recombination due to thermal fluctuations. Both noise sources are usually assumed to be white Gaussian noise.

The equation for the electric field can be written in terms of the optical intensity $`I`$ and the phase $`\varphi `$ by defining $`E=\sqrt{I}e^{\mathrm{i}\varphi }`$. For simplicity, we neglect the explicit random fluctuations terms and retain, as usual , the mean power of the spontaneous emission. The equations are:

$$\frac{dI}{dt}=(G(N,I)\gamma )I+4\beta N$$
(22)   
$$\frac{d\varphi }{dt}=\frac{1}{2}(G(N,I)\gamma )\alpha $$
(23)   
$$\frac{dN}{dt}=C\gamma _eNG(N,I)I$$
(24)   
$`G(N,I)`$ is the material gain given by:

$$G(N,I)=\frac{g(NN_o)}{1+sI}$$
(25)   
The definitions and typical values of the parameters for semiconductor lasers are given in Table 1. The first term of Eq. (22) accounts for the stimulated emission while the second accounts for the mean value of the spontaneous emission power. Eqs. (22 \- 24) are written in the reference frame in which the frequency of the on state is zero when spontaneous emission noise is neglected. The threshold condition is obtained by setting $`G(N,I)=\gamma `$, $`I=0`$ and neglecting spontaneous emission, i.e. $`N_{th}=N_o+\frac{\gamma }{g}`$. The threshold carrier injected per unit time to turn the laser on is given by $`C_{th}=\gamma _eN_{th}`$. Eq. (23) shows that $`\dot{\varphi }`$ is linear with $`N`$ and slightly (due to the smallness of the saturation parameter $`s`$, see Table 1) nonlinear with $`I`$.

Since in the deterministic case considered henceforth, the evolution equations for $`I`$ and $`N`$ do not depend on the phase $`\varphi `$, we can concentrate only in the evolution of $`I`$ and $`N`$. One can obtain a set of simpler dimensionless equations by performing the following change of variables:

$$y=\frac{2g}{\gamma }I,z=\frac{g}{\gamma }(NN_o),\tau =\frac{\gamma }{2}t$$
(26)   
The equations become then:

$`{\displaystyle \frac{dy}{d\tau }}`$ $`=`$ $`2\left({\displaystyle \frac{z}{1+\overline{s}y}}1\right)y+cz+d`$ (27)   
$`{\displaystyle \frac{dz}{d\tau }}`$ $`=`$ $`abz{\displaystyle \frac{zy}{1+\overline{s}y}},`$ (28)   
where we have defined $`a=\frac{2g}{\gamma ^2}(C\gamma _eN_o)`$, $`b=\frac{2\gamma _e}{\gamma }`$, $`c=\frac{16\beta }{\gamma }`$, $`d=\frac{16\beta gN_o}{\gamma ^2}`$ and $`\overline{s}=\frac{s\gamma }{2g}`$. These equations form the basis of our subsequent analysis. The steady states are obtained by setting (27) and (28) equal to zero, i.e.:

$`y_{st}`$ $`=`$ $`{\displaystyle \frac{1}{4(1+b\overline{s})}}[2(ab)+d(1+b\overline{s})+ca\overline{s}+\sqrt{v}]`$ (29)   
$`z_{st}`$ $`=`$ $`{\displaystyle \frac{a(1+\overline{s}y_{st})}{b+y_{st}(1+b\overline{s})}}`$ (30)   
where the constant $`v`$ is given by:

$`v`$ $`=`$ $`4(ab)^2+4d(a+b)(1+b\overline{s})+d^2(1+b\overline{s})^2`$ (31)   
$`+`$ $`c(8a+4a\overline{s}(a+b)+2da\overline{s}(1+b\overline{s}))+c^2a^2\overline{s}^2`$ (32)   
There is another steady state solution for $`y_{st}`$ given by Eq. (29) (with a minus sign in front of $`\sqrt{v}`$) which, however, does not correspond to any possible physical situation, since $`y_{st}<0`$. For a value of the injected carriers per unit time below threshold ($`C<C_{th}`$, equivalent to $`ab<0`$) $`y_{st}`$ is very small. This corresponds to the off solution in which the only emitted light corresponds to the spontaneous emission. Above threshold, stimulated emission occurs and the laser operates in the on state with large $`y_{st}`$. In what follows, we will concentrate in the evolution following the laser switch-on to the on state.

It is known that the dynamical evolution of $`y`$ and $`z`$ is such that they both reach the steady state by performing damped oscillations  whose period decreases with time. This fact is different from the usual relaxation oscillations that are calculated near the steady state by linearizing the dynamical equations. The time evolution of $`y`$ and $`z`$ is shown in Fig. 3a for some parameters (for another values of the parameters equivalent results are obtained), while the corresponding projection in the $`y`$, $`z`$ phase-plane is shown in Fig. 4. We are interested in obtaining a Lyapunov potential that can helps to explain the observed dynamics. This study was done in  without considering neither the saturation term, nor the mean value of the spontaneously emission power, and under those conditions an expression for the period of the transient oscillations was obtained. In our work, we calculate the period of the oscillations by taking into account these two effects. The period is obtained in terms of the potential, by assuming that the latter has a constant value during one period. It will be shown that this assumption works reasonably well and gives a good agreement with numerical calculations. Near the steady state, the relaxation oscillations can also be calculated in this form, but the potential is almost constant and consequently so is the period.

The evolutions equations (27), (28) can be casted in the form of Eq. (4) with the following Lyapunov potential:

$$V(y,z)=A_1y+A_2y^2+A_3\mathrm{ln}(y)+\frac{A_4}{y}+\frac{1}{2}B^2(y,z)$$
(33)   
where

$`A_1`$ $`=`$ $`{\displaystyle \frac{1}{2}}{\displaystyle \frac{1}{2}}a\overline{s}+b\overline{s}{\displaystyle \frac{1}{4}}\overline{s}d(1+b\overline{s}){\displaystyle \frac{1}{4}}a\overline{s}^2c`$ (34)   
$`A_2`$ $`=`$ $`{\displaystyle \frac{\overline{s}}{4}}(1+b\overline{s})`$ (35)   
$`A_3`$ $`=`$ $`{\displaystyle \frac{1}{2}}[ab+(ac+bd)\overline{s}+{\displaystyle \frac{d}{2}}]`$ (36)   
$`A_4`$ $`=`$ $`{\displaystyle \frac{(ac+bd)}{4}}`$ (37)   
$`B(y,z)`$ $`=`$ $`z1\overline{s}y+{\displaystyle \frac{(d+cz)}{2y}}(1+\overline{s}y).`$ (38)   
The corresponding (non-constant) matrix $`D`$ is given by:

$$D=\left(\begin{array}{cc}0& d_{12}\\ d_{12}& d_{22}\end{array}\right).$$
(39)   
being

$`d_{12}`$ $`=`$ $`{\displaystyle \frac{4y^2}{(1+\overline{s}y)[2y+c(1+\overline{s}y)]}}`$ (40)   
$`d_{22}`$ $`=`$ $`{\displaystyle \frac{4y[(1+2\overline{s}+b\overline{s})y^2+by+d+cz]}{(1+\overline{s}y)[2y+c(1+\overline{s}y)]^2}}`$ (41)   
This potential reduces to the one obtained in ref.  when setting $`c=d=\overline{s}=0`$ (which corresponds to set the laser parameters $`\beta =s=0`$). As expected, non-vanishing values for the parameters $`s`$ and $`\beta `$ increase the dissipative part of the potential ($`d_{22}`$), associated with the damping term. This result was pointed out in  when linearizing the rate equations around the steady state.

The equipotential lines of (33) are also plotted in Fig. 4. It is observed that there is only one minimum for V and hence the only stable solution (for this range of parameters) is that the laser switches to the on state and relaxes to the minimum of $`V`$. The movement towards the minimum of $`V`$ has two components: a conservative one that produces closed equipotential trajectories and a damping that decreases the value of the potential. The combined effects drives the system to the minimum following a spiral movement, best observed in Fig. 4.

The time evolution of the potential is also plotted in Fig. 3b. In this figure it can be seen that the Lyapunov potential is approximately constant between two consecutive peaks of the relaxation oscillations (This fact can be also observed with the equipotential lines of Fig. 4). This fact allows us to estimate the relaxation oscillation period by approximating $`V(y,z)=V`$, constant, during this time interval. When the potential is considered to be constant, the period can be evaluated by the standard method of elementary Mechanics: $`z`$ is replaced by its expression obtained from (27) in terms of $`y`$ and $`\dot{y}`$ (the dot stands for the time derivative) in $`V(y,z)`$. Using the condition that $`V(y,z)=V=constant`$, we obtain an equation for $`y`$ of the form: $`F(y,\dot{y})=V`$. From this equation, we can calculate the relaxation oscillation period ($`T`$) by integrating over a cycle. This leads to the expression:

$$T=_{y_0}^{y_1}\frac{1+\overline{s}y}{y}\frac{dy}{[2(VA_1yA_2y^2A_3\mathrm{ln}(y)A_4y^1)]^{1/2}}$$
(42)   
where $`y_0`$ and $`y_1`$ are the values of $`y`$ that cancel the denominator. We stress the fact that the only one approximation used in the derivation of this expression is that the Lyapunov potential is constant during two maxima of the intensity oscillations. The previous equation for the period reduces, in the case $`c=d=\overline{s}=0`$, to the one previously obtained by using the relation between the laser dynamics and the Toda oscillator derived in . Evaluation of the above integral shows that the period $`T`$ decreases as the potential $`V`$ decreases. Since the Lyapunov potential decreases with time, this explains the fact that the period of the oscillations in the transient regime decreases with time. In Fig. 5 we compare the results obtained with the above expression for the period with the one obtained from numerical simulations of the rate equations (27), (28). In the simulations we compute the period as the time between two peaks in the evolution of the variable $`y`$. As can be seen in this figure, the above expression for the period, when using the numerical value of the potential $`V`$, accurately reproduces the simulation results although it is systematically lower than the numerical result. The discrepancy is less than one percent over the whole range of times.

It is possible to quantify the difference between the approximate expression (42) and the exact values near the stationary state. In this case expression (42) reduces to :

$$T=\frac{2\pi }{d_{12,st}\sqrt{EFH^2}}$$
(43)   
where:

$`E`$ $`=`$ $`2\left(A_2{\displaystyle \frac{1}{2}}{\displaystyle \frac{A_3}{y_{st}^2}}+{\displaystyle \frac{A_4}{y_{st}^3}}+{\displaystyle \frac{1}{2}}\left[\overline{s}+{\displaystyle \frac{(d+cz_{st})}{2y_{st}^2}}\right]^2\right)`$ (44)   
$`F`$ $`=`$ $`\left[1+c{\displaystyle \frac{(1+\overline{s}y_{st})}{2y_{st}}}\right]^2`$ (45)   
$`H`$ $`=`$ $`\left[1+{\displaystyle \frac{c(1+\overline{s}y_{st})}{2y_{st}}}\right]\left[\overline{s}+{\displaystyle \frac{(d+cz_{st})}{2y_{st}^2}}\right]`$ (46)   
and $`d_{12,st}`$ is the coefficient $`d_{12}`$ calculated in the steady state. The period of the relaxation oscillations near the steady state can be obtained by linearizing eqs. (27) and (28) after a small perturbation is applied. The frequency of the oscillations in the steady state is the imaginary part of the eigenvalues of the linearized equations. This yields a period:

$$T_{st}=\frac{2\pi }{d_{12,st}\sqrt{EFH^2}}\left[1\frac{d_{22,st}^2}{d_{12,st}^2}\frac{F^2}{4(EFH^2)}\right]^{1/2}$$
(47)   
The difference between (43) and (47) vanishes with $`d_{22,st}`$ (i.e. $`d_{22}`$ in the stationary state). Since $`EFH^2`$ is always a positive quantity, our approximation will give, at least asymptotically, a smaller value for the period.

In order to have a complete understanding of the variation of the period with time, we need to compute the time variation of the potential $`V(\tau )`$ between two consecutive intensity peaks. This variation is induced by the dissipative terms in the equations of motion. We have not been able to derive an expression for the variation of the potential (see  for an approximate expression in a simpler case). However, we have found that a semi-empirical argument can yield a very simple law which is well reproduced by the simulations. We start by studying the decay to the stationary state in the linearized equations. By expanding around the steady state: $`y=y_{st}+\delta y`$, $`z=z_{st}+\delta z`$, the dynamical equations imply that the variables decay to the steady state as: $`\delta y(\tau )`$, $`\delta z(\tau )\mathrm{exp}(\frac{\rho }{2}\tau )`$, where:

$$\rho =d_{22,st}F$$
(48)   
Expanding $`V(y,z)`$ around the steady state and taking an initial condition at $`\tau _0`$ we find an expression for the decay of the potential:

$$\mathrm{ln}[V(\tau )V_{st}]=\mathrm{ln}[V(\tau _0)V_{st}]\rho (\tau \tau _0)$$
(49)   
In Fig. 6 we plot $`\mathrm{ln}[V(\tau )V_{st}]`$ versus time and compare it with the approximation (49). One can see that the latter fits $`\mathrm{ln}[V(\tau )V_{st}]`$ not only near the steady state (where it was derived), but also during the transient dynamics. The value of $`\tau _0`$, being a free parameter, was chosen at the time at which the first peak of the intensity appears. Although other values of $`\tau _0`$ might produce a better fit, the one chosen here has the advantage that it can be calculated analytically by following the technique of ref. . It can be derived from Eq. (42) that the period $`T`$ is linearly related to the potential $`V`$. This, combined with the result of Eq. (49), suggests the semi-empirical law for the evolution of the period of the form:

$$\mathrm{ln}[T(\tau )T_{st}]=\mathrm{ln}[T(\tau _0)T_{st}]\rho (\tau \tau _0)$$
(50)   
This simple expression fits well the calculated period not only near the steady state, but also in the transient regime, see Figs. 5 and 7. The tiny differences observed near the steady state are due to the fact that the semi-empirical law, Eq. (50), is based on the validity of relation Eq. (42) between the period and the potential. As it was already discussed above, that expansion slightly underestimates the asymptotic (stationary) value of the period. By complementing this study with the procedure given in  to describe the switch-on process of a laser, and valid until the first intensity peak is reached, we can obtain a complete description of the variation of the oscillations period in the dynamical evolution following the laser switch-on.

## IV Summary

In this work we have used Lyapunov potentials in the context of laser dynamics. For class A lasers, we have explained qualitatively the observed features of the deterministic dynamics by the movement on the potential landscape. We have identified the relaxational and conservative terms in the dynamical equations of motion. In the stochastic dynamics (when additive noise is added to the equations), we have explained the presence of a “noise sustained flow” for the phase of the electric field as the interaction of the conservative terms with the noise terms. An analytical expression allows the calculation of the phase drift.

In the case of class B lasers, we have obtained a Lyapunov potential valid only in the deterministic case, when noise fluctuations are neglected. We have found that the dynamics is non-relaxational with a non-constant matrix $`D`$. The fixed point corresponding to the laser in the on state is interpreted as a minimum in the potential landscape. By observing that the potential is nearly constant between two consecutive intensity peaks during the relaxation process towards the steady state, but still in a highly non-linear regime, we were able to obtain an approximate expression for the period of the oscillations. Moreover, we have derived a simple exponential approach of the period of the oscillations with time towards the period of the relaxation oscillations near the steady state. This dependence appears to be valid after the first intensity peak following the switch-on of the laser. A possible extension of our work could be to consider the presence of an external field, which is numerically studied in .

## Acknowledgments

We wish to thank Professor M. San Miguel and Professor G.L. Oppo for a careful reading of this manuscript and for useful comments. We acknowledge financial support from DGES (Spain) under Project Nos. PB94-1167 and PB97-0141-C02-01.

Table 1

| PARAMETERS |  | VALUES |
| --- | --- | --- |
| $`C`$ | Carriers injected per unit time. | $`>\mathrm{threshold}`$ |
| $`\gamma `$ | Cavity decay rate. | $`0.5ps^1`$ |
| $`\gamma _e`$ | Carrier decay rate. | $`0.001ps^1`$ |
| $`N_o`$ | Number of carriers at transparency. | $`1.5\times 10^8`$ |
| $`g`$ | Differential gain parameter. | $`1.5\times 10^8ps^1`$ |
| $`s`$ | Saturation parameter. | $`10^810^7`$ |
| $`\beta `$ | Spontaneous emission rate. | $`10^8ps^1`$ |
| $`\alpha `$ | Linewidth enhancement factor. | 3-6 |

## Figure captions

Fig. 1: Potential for a class A laser, Eq. (15) with the parameters: $`a=2`$, $`b=1`$. Dimensionless units. Fig. 2: Time evolution of the mean value of the phase $`\varphi `$ in a class A laser, in the case $`a=2`$, $`b=1`$, $`ϵ=0.1`$. For $`\alpha =0`$ (dashed line) there is only phase diffusion and the average value is $`0`$ for all times. When $`\alpha =5`$ (solid line) there is a linear variation of the mean value of the phase at late times. Error bars are incluted for some values. The dot-dashed line has the slope given by the theoretical prediction Eq. (21). The initial condition is taken as $`x_1=x_2=0`$ and the results were averaged over 10000 trajectories with different realizations of the noise. Dimensionless units. Fig. 3: a) Normalized intensity, $`y`$ (solid line) and normalized carriers number, $`z/40`$ (dot-dashed line) versus time in a class B laser obtained by numerical solution of Eqs. (27) and (28). b) Plot of the potential (33). Parameters: $`a=0.009`$, $`b=0.004`$, $`\overline{s}=0.5`$, $`c=3.2\times 10^9`$, $`d=1.44\times 10^8`$ which correspond to physical parameters in Table 1 with $`C=1.2C_{th}`$. The initial conditions are taken as $`y=5\times 10^8`$ and $`z=0.993`$. Dimensionless units. Fig. 4: Number of carriers versus intensity (scaled variables). The vector field and contour plot (thick lines) are also represented. Same parameters than in Fig. 3. Dimensionless units. Fig. 5: Period versus time in a class B laser. Solid line has been calculated as the distance between two peaks of intensity, with triangles plotted at the begining of each period; dashed line has been calculated using the expression (42), with the value of the potential $`V`$ obtained also from the simulation; dotted line corresponds to the semi-empirical expression (50). Same parameters than in Fig. 3. We have used $`\tau _0=55.55`$, coinciding with the position of the first intensity peak. Dimensionless units. Fig. 6: Logarithm of the potential difference versus time in a class B laser (solid line), compared with the theoretical expression in the steady state (49) (dashed line). Same parameters than in Fig. 3 and $`\tau _0`$ as in Fig. 5. Dimensionless units. Fig. 7: Logarithm of the period difference versus time in a class B laser. Triangles correspond to the period calculated from the simulations as the distance between two consecutive intensity peaks, at the same position than in Fig. 6. The dashed line is the semiempirical expression Eq. (50). Same parameters than in Fig. 3 and $`\tau _0`$ as in Fig. 5. Dimensionless units.