Mathematical ultrashort-pulse laser physics

V.L. Kalashnikov,

Institut fuer Photonik, TU Wien, Gusshausstrasse 27/387, A-1040 Vienna, Austria

vladimir.kalashnikov@tuwien.ac.at 

http://www.geocities.com/optomaplev

Abstract:  The analytical and numerical approaches to the analysis of the ultrashort-pulse solid-state lasers are presented. The unique self-consistent method of the laser dynamics analysis is based on the symbolical, numerical, programming, and graphical capacities of Maple 6. The algorithmization of basic conceptions as well as sophisticated research methods is of interest to both students and experts in the laser physics and the nonlinear dynamics.

Application Areas/Subjects : Optics, Laser Physics, Nonlinear Physics, Differential Equations, and Programming

Keywords:  soliton, ultrashort pulse, mode locking, Q-switching, solid-state laser, nonlinear Schroedinger equation, nonlinear Landau-Ginzburg equation, nonlinear Klein-Gordon equation, harmonic oscillator, nonlinear oscillator, self-phase modulation, group-delay dispersion, self-induced transparency, stimulated Raman scattering

Introduction

The ultrashort laser pulses, i.e. the pulses with the durations ~ -  sec, have a lot of the applications, which range ultrafast spectroscopy,  tracing chemical reactions, precision processing of materials, optical networks and computing, nuclear fusion and X - ray lasing, ophthalmology and surgery (for review see T. Brabec, F. Krausz , "Intense few-cycle laser fields: Frontiers of nonlinear optics", Rev. Mod. Phys. 72 , 545 (2000)). The mechanisms of the ultrashort pulse generation are active or passive loss switching (so-called Q-switching, part I) and locking of the longitudinal laser modes (part II) due to the active (part IV) or passive ultrafast modulation resulting in the laser quasi-soliton formation. Such quasi-soliton is very similar to the well-known Schroedinger soliton, which runs in the optical networks (part V). As a matter of fact, the model describing active mode locking is based on the usual equation of the harmonic oscillator or its nonlinear modifications, while for the passive mode locking description we have primordially nonlinear Landau-Ginzburg equation (part VI). This equation is the dissipative analog of the nonlinear Schroedinger equation and, as a result of the nonlinear dissipation, there exist a lot of nonstationary regimes of the ultrashort pulse generation (part VII). This requires the generalization of the model, which leads to the numerical simulations on the basis of FORTRAN (or C) codes generated by Maple (part VIII). Simultaneously, the obtained numerical results are supported by the analytical modelling in the framework of the computer algebra approach. The last takes into account the main features of the nonlinear dissipation in the mode-locked laser, viz. the power- (part VI) or energy-dependent (part IX) response of the loss to the generation field. In the latter case, there is the possibility of the so-called self-induced transparency formation, which is described by the nonlinear Klein-Gordon equation (part X).

Our considerations are based on the analytical or semi-analytical search of the steady-state soliton-like solutions of the laser dynamical equations and on the investigation of their stability in the framework of linear stability analysis. Also, the breezer-like solutions are considered using the aberrationless approximation. The computer algebra analysis is supported by the numerical simulations on the basis of the Maple generated FORTRAN-code. We present the analysis of these topics by means of the powerful capacities of Maple 6 in the analytical and numerical computations. This worksheet contains some numerical blocks, which can take about of 12 Mb and 18 min of computation (1 GHz Athlon).

Author is Lise Meitner Fellow at TU Vienna and would like to acknowledge the support of Austrian Science Fund's grant #M611. I am grateful to Dr. I.T. Sorokina and Dr. E. Sorokin for their hospitality at the Photonics Institute (TU Vienna), for stimulating discussions and experimental support of Part VIII. Also, I thank D.O. Krimer for his help in programming of part V.

Contents:

1. Nonstationary lasing: passive Q-switching

2. Conception of mode locking

3. Basic model

4. Active mode locking: harmonic oscillator

5. Nonlinear Schroedinger equation: construction of the soliton solution using the direct Hirota's method

6. Nonlinear Landau-Ginzburg equation: quasi-soliton solution

7. Autooscillations of quasi-solitons in the laser

8. Numerical approaches: ultrashort pulse spectrum, stability and multipulsing

9. Mode locking due to a "slow" saturable absorber

10. Coherent pulses: self-induced transparency in the laser

Part I. Nonstationary lasing: passive Q-switching

Continuous-wave oscillation

The basic principle of Q-switching is rather simple, but in the beginning let's consider the steady-state oscillation of laser. The near-steady-state laser containing an active medium and pumped by an external source of the energy (lamp, other laser or diode, for example) obeys the following coupled equations:

>	restart:  with(linalg):    eq1 := diff(Phi(t),t) = (alpha(t) - rho)*Phi(t);# field evolution   eq2 := diff(alpha(t),t) = sigma[p]*(a[m]-alpha(t))*P/(h*nu[p]) - alpha(t)*sigma[g]*Phi(t)/(h*nu[g]) - alpha(t)/T[r];# gain evolution for quasi-two-level active medium

Warning, the protected names norm and trace have been redefined and unprotected

Here  is the time-dependent field intensity,  is the dimensionless gain coefficient, P  is the time-independent (for simplicity sake) pump intensity,  and  are the frequencies of the pump and generation fields, respectively,  and  are the absorption and generation cross-sections, respectively,  is the gain relaxation time,  is the linear loss inclusive the output loss of the laser cavity, and, at last,  is the gain coefficient for the full population inversion in the active medium. The pump increases the gain coefficient (first term in eq2 ), that results in the laser field growth (first term in eq1 ). But the latter causes the gain saturation (second term), which can result in the steady-state operation (so-called continuous-wave, or simply cw, oscillation):

>	rhs( subs({alpha(t)=alpha,Phi(t)=Phi},eq1) ) = 0;  rhs( subs({alpha(t)=alpha,Phi(t)=Phi},eq2) ) = 0;   sol := solve({%,%%},{Phi,alpha});

The second solution defines the cw intensity, which is the linear function of pump intensity:

>

expand( subs( sol[2],Phi ) ):  Phi = collect(%,P);   plot( subs( {lambda[g]=8e-5, lambda[p]=5.6e-5},subs( {h=6.62e-34, sigma[g]= 3e-19, sigma[p]= 1e-19, rho=0.1, nu[g]=3e10/lambda[g], nu[p]=3e10/lambda[p], T[r]=3e-6, a[m]=2.5},rhs(rho*%) ) ), P=4.9e4..1e5, axes=boxed, title=`output intensity vs. pump, [W/cm^2]` );# lambda is the wavelength in [cm]

The pump corresponding to =0 defines so-called generation threshold.
Now let's consider the character of the steady-state points of our system { eq1 , eq2 } presented by sol . The Jacobian of the system { eq1 , eq2 } is:

>	eq3 := subs( {Phi(t)=x,alpha(t)=y},rhs(eq1) ):# x = Phi(t), y = alpha(t)  eq4 := subs( {Phi(t)=x,alpha(t)=y},rhs(eq2) ):   A := vector( [eq3, eq4] );# vector made from the right-hand side of system {eq1, eq2}    B := jacobian(A, [x,y]);

For the cw-solution the eigenvalues of the perturbations are:

>

BB := eigenvalues(B):  BBB := subs( {x=subs( sol[2],Phi ),y=subs( sol[2],alpha )},BB[1] ):   BBBB := subs( {x=subs( sol[2],Phi ),y=subs( sol[2],alpha )},BB[2] ):  plot( [subs( \ {lambda[g]=8e-5, lambda[p]=5.6e-5}\ ,subs( \ {h=6.62e-34, sigma[g]= 3e-19, sigma[p]= 1e-19, rho=0.1, nu[g]=3e10/lambda[g], nu[p]=3e10/lambda[p], T[r]=3e-6, a[m]=2.5}\ ,BBB ) ),\ subs( \ {lambda[g]=8e-5, lambda[p]=5.6e-5}\ ,subs( \ {h=6.62e-34, sigma[g]= 3e-19, sigma[p]= 1e-19, rho=0.1, nu[g]=3e10/lambda[g], nu[p]=3e10/lambda[p], T[r]=3e-6, a[m]=2.5}\ ,BBBB/1e7 ) )], P=4.9e4..1e5, axes=boxed, title=`perturbation eigenvalues vs. pump` );# second eigenvalue is divided by 1e7

So, cw oscillations is stable in our simple case because the eigenvalues are negative. For the zero field solution the perturbation eigenvalues are:

>

BB := eigenvalues(B):  BBB := subs( {x=subs( sol[1],Phi ),y=subs( sol[1],alpha )},BB[1] ):   BBBB := subs( {x=subs( sol[1],Phi ),y=subs( sol[1],alpha )},BB[2] ):  plot( [subs( \ {lambda[g]=8e-5, lambda[p]=5.6e-5}\ ,subs( \ {h=6.62e-34, sigma[g]= 3e-19, sigma[p]= 1e-19, rho=0.1, nu[g]=3e10/lambda[g], nu[p]=3e10/lambda[p], T[r]=3e-6, a[m]=2.5}\ ,BBB ) ),\ subs( \ {lambda[g]=8e-5, lambda[p]=5.6e-5}\ ,subs( \ {h=6.62e-34, sigma[g]= 3e-19, sigma[p]= 1e-19, rho=0.1, nu[g]=3e10/lambda[g], nu[p]=3e10/lambda[p], T[r]=3e-6, a[m]=2.5}\ ,BBBB/1e7 ) )], P=4.9e4..1e5, axes=boxed, title=`perturbation eigenvalues vs. pump` );# second eigenvalue is divided by 1e7

The existence of the positive eigenvalue suggests the instability of the zero-field steady-state solution. Hence there is the spontaneous generation of the cw oscillation above threshold in the model under consideration.

Q-switching

The situation changes radically due to insertion of the saturable absorber into laser cavity. In this case, in addition to the gain saturation, the loss saturation appears. This breaks the steady-state operation and produces the short laser pulses. As a result of the additional absorption, Q-factor of laser is comparatively low (high threshold). This suppresses the generation. When  is small, the gain increases in the absence of the gain saturation (see eq2 from the previous subsection). This causes the field growth. The last saturates the absorption and abruptly increases Q-factor. The absorption "switching off" leads to the explosive generation, when the most part of the energy, which is stored in the active medium during pumping process, converts into laser field. The increased field saturates the gain and this finishes the generation.

As the reference for the model in question see, for example, J.J. Degnan, "Theory of the optimally coupled Q-switched laser", IEEE J. Quant. Elect. 25 , 214   (1989). To formulate the quantitative model of the laser pulse formation let's use the next approximations: 1) the pulse width is much larger than the cavity period, and 2) is less than the relaxation time, 3) the pump action during the stage of the pulse generation is negligible. We shall use the quasi-two level schemes for the gain and loss media (the relaxation from the intermediate levels is fast). Also, the excited-state absorption in absorber will be taken into account.

The system of equation describing the evolution of the photon density  is

>

restart: with(plots): print(`System of basic equations:`); e1 := Diff(n[1](t),t) =\       -sigma[s]*c*phi(t)*n[1](t);# EVOLUTION OF ABSORPTION. The ground level population n[1](t) defines the absorption (quasi-two level scheme). The relaxation is slow in the comparison with the pulse duration, phi(t) is the photon density   e2 := Diff(phi(t),t) =\      (phi(t)/t_cav)*\      (2*sigma[g]*x(t)*l - log(1/R) - L - 2*sigma[s]*n[1](t)*l[s]-2*sigma[esa]*(n[0]-n[1](t))*l[s]);# EVOLUTION OF PHOTON DENSITY. The field variation over cavity round-trip is small, sigma[g] is the gain cross-section, x(t) is the inversion in amplifier defining the gain coefficient, t_cav is the cavity period, l is the active medium length, R is the output coupler refractivity, L is the linear loss, l[s] is the absorber length          e3 := Diff(x(t),t) =\               -gamma*sigma[g]*c*phi(t)*x(t);# EVOLUTION OF GAIN. c is the light velocity, gamma is the parameter of the inversion reduction (2 for the pure three-level scheme and 1 for the pure four-level scheme)

Warning, the name changecoords has been redefined

Now we shall search the ground state population in the absorber as a function of the initial population inversion in amplifier:

>

Diff(n1(n),n) = subs({n[1](t)=n[1](n),x(t)=x},rhs(e1)/rhs(e3));# devision of e1 by e3  Int(1/y,y=n[0]..n[1]) = Int(zeta/z,z=x_i..x);# zeta = sigma[s]/gamma/sigma[g], x_i is the initial inversion in amplifier, n[0] is the concentration of active ions in absorber    solve(value(%),n[1]):     n[1] = expand(%);      print(`Ground state population in absorber:`);       sol_1 := n[1] = n[0]*(x/x_i)^zeta;

The similar manipulation allows to find the photon density as a function of inversion in amplifier:

>

A := Diff(phi(x),x) =\      simplify( subs( n[1](t)=rhs(sol_1),subs( {x(t)=x},rhs(e2)/rhs(e3) ) ) );# devision of e2 by e3   B := numer( rhs(A) );    C := denom( rhs(A) );     BB := subs( {op(4,B)=(x/x_i)^(zeta)*log(1/T[0]^2),2*sigma[esa]*l[s]*n[0]=ln(1/T[s]^2),op(6,B)=-ln(1/T[s]^2)*(x/x_i)^zeta}, B );# T[0] is the initial transmission of the absorber (log(1/T[0]^2)= 2*sigma[s]*l[s]*n[0], l[s] is the absorber length)       CC := 2*l_cav*gamma*sigma[g]*x;# l_cav=t_cav*c/2 is the cavity length, t_cav is the cavity period

>

print(`Evolution of the photon density:`);        e4 := diff(phi(x),x) = subs(ln(1/(T[s]^2))=delta*ln(1/(T[0]^2)),BB)/CC;# delta=sigma[esa]/sigma[s]=ln(T[s])/ln(T[0]) is the parameter defining the contribution of an excited-state absorption with cross-section sigma[esa], T[s] is the fully saturated transmission of absorber

Hence the photon density is:

>	dsolve({e4,phi(x_0)=0},phi(x)):  simplify( subs(x_0=x_i,%) ):   print(`This is the basic dependence for photon density:`);    sol_2 := collect(combine(%,ln),{log(1/T[0]^2),sigma[g],zeta});

So, we have:

( x )=          ( Eq. 1 )  

Now let's define the key Q-switching parameters:

>	subs( n[1](t)=n[0],rhs(e2)*t_cav/phi(t)) = 0;# Q-switching start, n[1](0) = n[0]  print(`Solution for the initial inversion:`);   sol_3 := x_i = subs( 2*sigma[s]*n[0]*l[s]=log(1/T[0]^2), solve(%,x(t)) );

So, the initial inversion defining the gain at Q-switching start is

           ( Eq. 2 )

>	e5 := numer( rhs(e4) ) = 0;# definition of the pulse maximum  print(`The inversion at the pulse maximum:`);   sol_4 := x_t = solve(e5,x);

>	e6 := numer( simplify( rhs(sol_2) ) ) = 0;# definition of the Q-switching finish  print(`The inversion at Q-switching finish`);   sol_5 := x_f = solve(e6,x);

Additionally, we define the inversion at the pulse maximum, when   tends to infinity:

>	subs( (x/x_i)^zeta=0,e5 );# by virtue of x_i > x and zeta -> infinity  print(`Inversion at the pulse maximum for large zeta`);   sol_6 := x_t0 = solve(%,x);

So, we have the expressions for  (initial inversion, sol_3 ),  (the inversion at pulse maximum, sol_4  and e5 ),  (the inversion at pulse maximum when --> , sol_6 ),   (the final inversion, sol_5  and e6 ) and the photon density  as function of inversion x  ( sol_2 ).

As an example, we consider the real situation of Yb/Er-glass laser with the crystalline Co:MALO saturable absorber. The obtained expressions allow to plot the typical dependencies for the pulse parameters:

>	print(`Pulse energy:`);  fun_1 := (h*nu*S)/(2*sigma[g]*gamma)*log(1/R)*log(x_i/x_f);   print(`Pulse power:`);    fun_2 := (h*nu*S*l)/(t_cav*gamma)*log(1/R)*(x_i - x_t - x_t0*log(x_i/x_t) - (x_i-x_t0)*(1-(x_t/x_i)^alpha)/alpha);      print(`Pulse width:`); fun_3 := simplify(fun_1/fun_2);

This numerical procedure plots the dependence of the output pulse energy on the reflectivity of the output mirror:

>

En := proc(gam,L,T_0,l,w0,i) # definition of system's parameters delta := 0.028:  sigma[g] := 7e-21:# the gain cross-section in [cm^2]   sigma[s] := 3.5e-19:# the absorption cross-section in Co:MALO    alpha := sigma[s]/(sigma[g]*gam):     h := 6.63e-34:# J*s      nu := evalf(3e8/1.535e-6):# lasing frequency for 1.54 micrometers       S := Pi*w0^2/2:# area of Gaussian beam in amplifier, w0 is the beam radius      R := 0.5+0.5*i/100:         x_i := 1/2*(ln(1/R) + L + ln(1/(T_0^2)))/(l*sigma[g]):  eq := ln(x)*delta*ln(1/(T_0^2))*alpha-ln(x_i)*delta*ln(1/(T_0^2))*alpha-2*l*x*sigma[g]*alpha+2*l*x_i*sigma[g]*alpha+ln(x)*L*alpha+ln(x)*ln(1/R)*alpha-ln(x_i)*L*alpha-ln(x_i)*ln(1/R)*alpha+(x/x_i)^alpha*ln(1/(T_0^2))-delta*ln(1/(T_0^2))*(x/x_i)^alpha-ln(1/(T_0^2))+delta*ln(1/(T_0^2)) = 0:# e6 sol_f := fsolve(eq, x, avoid={x=0}):  sol_En := evalf( 1/2*h*nu*S*ln(1/R)*ln(x_i/sol_f)/(sigma[g]*gam) * 1e3 ):# [mJ] end: print(`The parameters are:`); print(`1) inversion reduction factor gamma`); print(`2) linear loss L`); print(`3) initial transmission of absorber T[0]`); print(`4) gain medium length l in cm`); print(`5) beam radius in amplifier w0 in cm`);

Warning, `delta` is implicitly declared local to procedure `En`

Warning, `sigma` is implicitly declared local to procedure `En`

Warning, `alpha` is implicitly declared local to procedure `En`

Warning, `h` is implicitly declared local to procedure `En`

Warning, `nu` is implicitly declared local to procedure `En`

Warning, `S` is implicitly declared local to procedure `En`

Warning, `R` is implicitly declared local to procedure `En`

Warning, `x_i` is implicitly declared local to procedure `En`

Warning, `eq` is implicitly declared local to procedure `En`

Warning, `sol_f` is implicitly declared local to procedure `En`

Warning, `sol_En` is implicitly declared local to procedure `En`

For the comparison we use the experimental data (crosses in Figure):

>

points := {seq([0.5+0.5*k/100,En(1.9,0.04,0.886,4.9,0.065,k)],k=1..100)}:  points_exp := {[0.793,10.5],[0.88,9],[0.916,5.5]}:  plot(points,x=0.5..1,style=point,axes=boxed,symbol=circle,color=black,title=`Pulse energy vs. R, [mJ]`):   plot(points_exp,x=0.5..1,style=point,symbol=cross,axes=boxed,color=black):    display({%,%%});

Similarly, for the output power we have:

>

Pow := proc(gam,L,T_0,l,l_cav,w0,i) # definition of system's parameters delta := 0.028:  sigma[g] := 7e-21:# in [cm^2]   sigma[s] := 3.5e-19:# Co:MALO    alpha := sigma[s]/(sigma[g]*gam):     h := 6.63e-34:# J*s      nu := evalf(3e8/1.354e-6):# frequency for 1.54 micrometers       S := Pi*w0^2/2:# area for Gaussian beam      R := 0.5 + 0.5*i/100:     c := 3e10:    refractivity := 1.6:# refractivity coefficient for the active medium   t_cav := 2*(l_cav - l)/c + 2*(l*refractivity)/c:     x_i := 1/2*(ln(1/R) + L + ln(1/(T_0^2)))/(l*sigma[g]): x_t0 := 1/2*(ln(1/R)+L+delta*ln(1/(T_0^2)))/(sigma[g]*l): eq := -2*sigma[g]*x*l+ln(1/R)+L+(x/x_i)^alpha*ln(1/(T_0^2))+delta*ln(1/(T_0^2))-delta*ln(1/(T_0^2))*(x/x_i)^alpha = 0:# e5 sol_t[i] := fsolve(eq, x, avoid={x=0}):  sol_Pow[i] := h*nu*S*l*ln(1/R)* (x_i-sol_t[i]-x_t0*ln(x_i/sol_t[i])-(x_i-x_t0)*(1-(sol_t[i]/x_i)^alpha)/alpha)/(t_cav*gam)/1e3:# [kW] end: print(`The parameters are:`); print(`1) inversion reduction factor gamma`); print(`2) linear loss L`); print(`3) initial transmission of absorber T[0]`); print(`4) gain medium length l in cm`); print(`5) cavity length l_cav in cm`); print(`6) beam radius in amplifier w0 in cm`);

Warning, `delta` is implicitly declared local to procedure `Pow`

Warning, `sigma` is implicitly declared local to procedure `Pow`

Warning, `alpha` is implicitly declared local to procedure `Pow`

Warning, `h` is implicitly declared local to procedure `Pow`

Warning, `nu` is implicitly declared local to procedure `Pow`

Warning, `S` is implicitly declared local to procedure `Pow`

Warning, `R` is implicitly declared local to procedure `Pow`

Warning, `c` is implicitly declared local to procedure `Pow`

Warning, `refractivity` is implicitly declared local to procedure `Pow`

Warning, `t_cav` is implicitly declared local to procedure `Pow`

Warning, `x_i` is implicitly declared local to procedure `Pow`

Warning, `x_t0` is implicitly declared local to procedure `Pow`

Warning, `eq` is implicitly declared local to procedure `Pow`

Warning, `sol_t` is implicitly declared local to procedure `Pow`

Warning, `sol_Pow` is implicitly declared local to procedure `Pow`

>	points := {seq([0.5+0.5*k/100,Pow(1.6,0.01,0.886,4.9,35,0.06,k)],k=1..100)}:  plot(points,x=0.5..1,style=point,symbol=circle,color=black,axes=boxed,title=`Pulse power vs. R, [kW]`);

And, at last, the pulse durations are:

>

Width := proc(gam,L,T_0,l,l_cav,i) # definition of system's parameters delta := 0.028:  sigma[g] := 7e-21:# in [cm^2]   sigma[s] := 3.5e-19:# Co:MALO    alpha := sigma[s]/(sigma[g]*gam):     h := 6.63e-34:# J*s      nu := evalf(3e8/1.354e-6):# frequency for 1.54 micrometers       S := Pi*w0^2/2:# area for Gaussian beam      R := 0.5 + 0.5*i/100:     c := 3e10:     refractivity := 1.5:# active medium    t_cav := 2*(l_cav-l)/c + 2*(l*refractivity)/c:        x_i := 1/2*(ln(1/R) + L + ln(1/(T_0^2)))/(l*sigma[g]): x_t0 := 1/2*(ln(1/R)+L+delta*ln(1/(T_0^2)))/(sigma[g]*l): eq1 := -2*sigma[g]*x*l+ln(1/R)+L+(x/x_i)^alpha*ln(1/(T_0^2))+delta*ln(1/(T_0^2))-delta*ln(1/(T_0^2))*(x/x_i)^alpha = 0:# e5 eq2 := ln(x)*delta*ln(1/(T_0^2))*alpha-ln(x_i)*delta*ln(1/(T_0^2))*alpha-2*l*x*sigma[g]*alpha+2*l*x_i*sigma[g]*alpha+ln(x)*L*alpha+ln(x)*ln(1/R)*alpha-ln(x_i)*L*alpha-ln(x_i)*ln(1/R)*alpha+(x/x_i)^alpha*ln(1/(T_0^2))-delta*ln(1/(T_0^2))*(x/x_i)^alpha-ln(1/(T_0^2))+delta*ln(1/(T_0^2)) = 0:# e6 sol_f := fsolve(eq2, x, avoid={x=0}): sol_t := fsolve(eq1, x, avoid={x=0}):  sol_Width := (-1/2*ln(x_i/sol_f)*t_cav*alpha/(sigma[g]*l*(-x_i*alpha+sol_t*alpha+x_t0*ln(x_i/sol_t)*alpha+x_i-x_i*(sol_t/x_i)^alpha-x_t0+x_t0*(sol_t/x_i)^alpha))) * 1e9:# [ns] end: print(`The parameters are:`); print(`1) reducing parameter gamma`); print(`2) linear loss L`); print(`3) initial transmission of absorber T[0]`); print(`4) gain medium length l in cm`); print(`5) cavity length l_cav in cm`);

Warning, `delta` is implicitly declared local to procedure `Width`

Warning, `sigma` is implicitly declared local to procedure `Width`

Warning, `alpha` is implicitly declared local to procedure `Width`

Warning, `h` is implicitly declared local to procedure `Width`

Warning, `nu` is implicitly declared local to procedure `Width`

Warning, `S` is implicitly declared local to procedure `Width`

Warning, `R` is implicitly declared local to procedure `Width`

Warning, `c` is implicitly declared local to procedure `Width`

Warning, `refractivity` is implicitly declared local to procedure `Width`

Warning, `t_cav` is implicitly declared local to procedure `Width`

Warning, `x_i` is implicitly declared local to procedure `Width`

Warning, `x_t0` is implicitly declared local to procedure `Width`

Warning, `eq1` is implicitly declared local to procedure `Width`

Warning, `eq2` is implicitly declared local to procedure `Width`

Warning, `sol_f` is implicitly declared local to procedure `Width`

Warning, `sol_t` is implicitly declared local to procedure `Width`

Warning, `sol_Width` is implicitly declared local to procedure `Width`

>

points := {seq([0.5+0.5*k/100,Width(1.9,0.04,0.89,4.9,35,k)],k=1..100)}:  points_exp := {[0.793,70],[0.88,70],[0.916,75]}:  plot(points,x=0.5..1,style=point,axes=boxed,symbol=circle,title=`Pulse width vs. R, [ns]`):  plot(points_exp,x=0.5..1,style=point,axes=boxed,color=black):  display({%,%%},view=50..100);

The worse agreement with the experimental data for the pulse durations is caused by the deviation of the pulse shape from the Gaussian profile, which was used for the analytical estimations. More precise consideration is presented in J. Liu, D. Shen, S.-Ch. Tam, and Y.-L. Lam, "Modelling pulse shape of Q-switched lasers", IEEE J. Quantum Electr.   37 , 888 (2001).

The main advantage of the analytical model in question is the potential of the Q-switched laser optimization without any cumbersome numerical simulations. Let's slightly transform the above obtained expressions:

>	e7 := subs( x=x_t,expand( e5/(2*sigma[g]*x_i*l) ) );  e8 := subs( x=x_f,expand( rhs(sol_2)*l_cav*gamma/l ) ) = 0;   sol_6;    sol_3;

>

# Hence  e9 := x_t/x_i = (ln(1/R)+L+delta*ln(1/(T[0]^2)))/(2*l*x_i*sigma[g]) + (x_t/x_i)^zeta*ln(1/T[0]^2)*(1-delta)/(2*l*x_i*sigma[g]);   e10 := x_f - x_i = ln(x_f/x_i)*(ln(1/R) + L + delta*ln(1/(T[0]^2)))/(2*l*sigma[g]) - ln(1/T[0]^2)*(1-delta)*(1-(x_f/x_i)^zeta)/(2*l*sigma[g]*alpha);    simplify( sol_6 - sol_3 );     e11 := subs( ln(1/T[0]^2)=(x_i-x_t0)*2*l*sigma[g]/(1-delta),subs(ln(1/R)+L+delta*ln(1/(T[0]^2))=x_t0*2*l*sigma[g],e9));      e12 := subs( ln(1/T[0]^2)=(x_i-x_t0)*2*l*sigma[g]/(1-delta),subs(ln(1/R)+L+delta*ln(1/(T[0]^2))=x_t0*2*l*sigma[g],e10));

Now, let's plot some typical dependence for Q-switching.

>	eq := subs({x_t0=x,x_i=1},subs(x_t/x_i=y,expand(e11)));# here x=x_t0/x_i, y=x_t/x_i

>	xt := proc(zeta,i) sol := array(1..20):   x := (i-1)/20:    sol[i] := fsolve(y = x+y^zeta-y^zeta*x,y,avoid={y=1},0..infinity): end:

Warning, `sol` is implicitly declared local to procedure `xt`

Warning, `x` is implicitly declared local to procedure `xt`

>

points := {seq([k/20,xt(1.5,k)],k=2..19)}: p1 := plot(points,x=0..1,style=point,axes=boxed):   points := {seq([k/20,xt(3,k)],k=2..19)}:    p2 := plot(points,x=0..1,style=point,axes=boxed):     points := {seq([k/20,xt(6,k)],k=2..19)}:      p3 := plot(points,x=0..1,style=point,axes=boxed):   display({p1,p2,p3},view=0..1,title=`x_t/x_i vs x_t0/x_i for different zeta`);

From the definition of    this dependence for growing  tends to the linear one, which correspond to the maximal efficiency of the population inversion utilization. For the fixed absorption and emission cross-sections, the  parameter increases as a result of the laser beam focusing in the saturable absorber. Then ( --> ), we have:

>	e13 := lhs(e11) = op(1,rhs(e11));  e14 := expand(lhs(e12)/x_i) = expand(op(1,rhs(e12))/x_i);

and the output energy optimization can be realized by this simple way:

>

#                         Energy optimization e15 := x_f/x_i-1 = ln(x_f/x_i)*rhs(sol_6)/x_i;# from e14 and expression for x_t0  e16 := lhs(%) = ln(x_f/x_i)*(a + b + delta*c);# a= ln(1/R)/(2*sigma[g]*x_i*l), b = L/(2*sigma[g]*x_i*l), c= ln(1/T[0]^2)/(2*sigma[g]*x_i*l) are the relative shares of the output, linear and saturable loss in the net-loss   expand(solve(%,a)):    sol_7 := (y/log(y) - 1/log(y) - (1-delta)*b - delta)/(1-delta);# solution for a, y = x_f/x_i, we used a+b+c=1     epsilon := -a*log(y);# dimensionless energy, epsilon = En*2*sigma[g]*gamma/(h*nu*A)/(2*sigma[g]*x_i*l);     simplify( subs(a=sol_7,epsilon) );    diff(%,y) = 0;# maximum of energy    print(`An optimal ratio of final and initial inversions:`);   y_opt := solve(%,y);# optimal y   print(`An optimal mirror:`);  a_opt := subs(y=y_opt,sol_7);# optimal a print(`A maximal dimensionless energy:`); epsilon_opt := simplify( subs({a=a_opt,y=y_opt},epsilon) );# optimal epsilon

>

p1 := plot({[b,subs(delta=0,a_opt),b=0..1],[b,1-subs(delta=0,a_opt)-b,b=0..1]},color=black,axes=boxed,title=`optimal a and c versus b`):     p2 := plot({[b,subs(delta=0.05,a_opt),b=0..1],[b,1-subs(delta=0.05,a_opt)-b,b=0..1]},color=red,axes=boxed,title=`optimal a and c versus b`):      p3 := plot({[b,subs(delta=0.1,a_opt),b=0..1],[b,1-subs(delta=0.1,a_opt)-b,b=0..1]},color=blue,axes=boxed,title=`optimal a and c versus b`):              display(p1,p2,p3,view=0..1);        p1 := plot([b,subs(delta=0,epsilon_opt),b=0..1], axes=boxed,color=black, title=`optimal epsilon versus b`):       p2 := plot([b,subs(delta=0.05,epsilon_opt),b=0..1], axes=boxed,color=red, title=`optimal epsilon versus b`):       p3 := plot([b,subs(delta=0.1,epsilon_opt),b=0..1], axes=boxed,color=blue, title=`optimal epsilon versus b`):        display(p1,p2,p3);

From the first figure the optimization can be performed by the simple graphical way. We have to define the appropriate for our scheme laser's net-loss and to determine (measure or calculate) the intracavity linear loss. This gives the value of b , which is the relation of the linear loss to the net-loss. The upper group of curves gives the value of c , the lower curves give a  (for the different contribution of the excited-state absorption, i. e. ).

Two-color pulsing

Now let consider a more complicated situation, which corresponds to the two-color Q-switching due to presence of the stimulated Raman scattering in the active medium (for example, Yb(3+):KGd(W0(4))(2), see A.A. Lagatsky, A. Abdolvand, N.V. Kuleshov, Opt. Lett. 25 , 616 (2000)). In this case the analytical modelling is not possible, but we can use the numerical capacities of Maple.

Let's the gain medium length is  and the Raman gain coefficient is g . We assume the exact Raman resonance and neglect the phase and group-velocity effects. Then the evolutional equations for the laser and scattered intensities  and , respectively, are:  

>

restart: with(DEtools): with(plots): eq1 := diff(In[p](z),z) = -g*In[s](z)*In[p](z):# evolution of the laser intensity  eq2 := diff(In[s](z),z) = g*In[p](z)*In[s](z):# evolution of the Stokes component intensity   sys := {eq1,eq2};    IC := {In[p](0)=In[p,0],In[s](0)=In[s,0]};# initial conditions   sol := dsolve(sys union IC,{In[s](z),In[p](z)}):

Warning, the name changecoords has been redefined

The integration produces:

>	subs(sol,In[s](z)):# z is the distance in the crystal  sol_1 := In[s](z) = simplify(%); subs(sol,In[p](z)):  sol_2 := In[p](z) = simplify(%);

There are the amplification of the scattered field and the depletion of the laser field, which plays a role of the pump for the Stokes component.

Now let's transit from the intensities to the photon densities:

>

#phi is the photon densities, the laser component at 1033 nm, the Stokes component at 1139 nm sol_3 := simplify( subs({z=2*l[g],In[p,0]=c*h*nu[p]*phi[p,0]/2,In[s,0]=c*h*nu[s]*phi[s,0]/2,In[s](z)=c*h*nu[s]*phi[s](z)/2},2*sol_1) )/(c*h*nu[s]);   sol_4 := simplify( subs({z=2*l[g],In[p,0]=c*h*nu[p]*phi[p,0]/2,In[s,0]=c*h*nu[s]*phi[s,0]/2,In[p](z)=c*h*nu[p]*phi[p](z)/2},2*sol_2) )/(c*h*nu[p]);

Hence the photon densities evolution due to the stimulated Raman scattering obeys:

>	eq3 := diff(phi[p](t),t)[raman] = ( subs( {phi[p,0]=phi[p](t),phi[s,0]=phi[s](t)},rhs(sol_4) ) - phi[p](t) )/t[cav];# t[cav] is the laser cavity period eq4 := diff(phi[s](t),t)[raman] = ( subs( {phi[p,0]=phi[p](t),phi[s,0]=phi[s](t)},rhs(sol_3) )- phi[s](t) )/t[cav];

The inverse and ground-state populations for the gain medium and the absorber (Cr(4+):YAG) obey (see previous subsection):

>	eq5 := diff(n(t),t) = -gamma*sigma[g]*c*phi[p](t)*n(t);  eq6 := diff(n[0](t),t) = -rho*sigma[a]*c*phi[p](t)*n[0](t);# rho=S[g]/S[a] is the ratio of the beam area in the gain medium to that in the absorber

Contribution of gain, saturable, output and linear intracavity loss to the field's evolution results in:

>	eq7 := diff(phi[p](t),t)[gain] = 2*(sigma[g]*n(t)*l[g] - sigma[a]*n[0](t)*l[a])*phi[p](t)/t[cav];  eq8 := diff(phi[p](t),t)[linear] = (-ln(1/R[p]) - L[p])*phi[p](t)/t[cav];   eq9 := diff(phi[s](t),t)[linear] = (-ln(1/R[s]) - L[s] - 2*kappa*sigma[a]*n[0](t)*l[a])*phi[s](t)/t[cav];

where  and  are the linear loss for the laser and Stokes components, respectively,  and  are the output mirror reflectivity at the laser and Stokes wavelengths, respectively,  is the reduction factor taking into account the decrease of the loss cross-section at the Stokes wavelength relatively to the lasing one.

Hence, the field densities evaluate by virtue of:

>	eq10 := diff(phi[p](t),t) = rhs(eq3) + rhs(eq7) + rhs(eq8);  eq11 := diff(phi[s](t),t) = rhs(eq4) + rhs(eq9);

In the agreement with the results of the previous subsection we can transform the equation for the populations in the gain and absorption media:

>	diff(n[0](n),n) = subs( {n(t)=n,n[0](t)=n[0](n)},rhs(eq6)/rhs(eq5) );  diff(n[0](n),n) = zeta*rho*n[0](n)/n;# zeta=sigma[a]/sigma[g]/gamma   dsolve({%},n[0](n));

>	eq12 := n[0] = n[0,i]*(n/n[i])^(zeta*rho);

>	simplify( subs( phi[s](t)=0,subs( {n(t)=n[i],n[0](t)=n[0,i],phi[p](t)=0},expand( rhs(eq10)/phi[p](t) ) ) ) ) = 0:# the generation start  numer(lhs(%)) = 0:   eq13 := n[i] = subs( 2*sigma[a]*n[0,i]*l[a]=ln(1/T[0]^2),solve(%,n[i]) );

Let' introduce new variables:

>

Phi[p](t) = 2*sigma[g]*l[g]*phi[p];# note, that Phi is not the dimensional intensity from the first subsection!  Phi[s](t) = 2*sigma[g]*l[g]*phi[s];   tau = 2*sigma[g]*l[g]*n[i]*t/t[cav];    Nu = nu[s]/nu[p];     G = g*l[g]*c*h*nu[p]/(2*sigma[g]*l[g]);      y(tau) = n(t)/n[i];       Xi = ln(1/(T[0]^2))+ln(1/R[p])+L[p];        kappa = 0.38;

Then from Eqs. ( eq5, eq10, eq11, eq12 ):

>

f[1] := diff(Phi[p](tau),tau) = Phi[p](tau)/Xi*( Xi*y(tau) - ln(1/T[0]^2)*y(tau)^(alpha*rho) - (ln(1/R[p])+L[p]) + ( (Phi[p](tau)+N*Phi[s](tau))/(Phi[p](tau)+exp(G*(Phi[p](tau)+N*Phi[s](tau)))*N*Phi[s](tau)) - 1 ) ); f[2] := diff(Phi[s](tau),tau) = Phi[s](tau)/Xi*(-(ln(1/R[s])+L[s]+kappa*ln(1/(T[0]^2))*y(tau)^(alpha*rho)) + ( Phi[p](tau)+N*Phi[s](tau))/(Phi[p](tau)*exp(-G*(Phi[p](tau)+N*Phi[s](tau)))+N*Phi[s](tau) ) - 1 ); f[3] := diff(y(tau),tau) = -gamma*(l[cav]/l[g])*Phi[p](tau)/Xi;

The next procedures solve the obtained system numerically to obtain the dependences of the normalized photon densities at laser and Stokes wavelengths and the relative inversion vs. cavity roundtrip, n is the duration of the simulation in the cavity roundtrips):

>

ODE_plot1 := proc(T0,Rp,Rs,Lp,Ls,n) gam := 1:# gamma parameter rho := 1: kappa := 0.38: alpha := evalf( 5e-18/2.8e-20 ):#the ratio of the cross-sections x := 29:#  x = l[cav]/l[g] N := evalf(1033/1139): G := subs( {g=4.8e-9,c=3e10,h=6.62e-34,nu_p=3e10/1.033e-4,sigma_g=2.8e-20},g*c*h*nu_p/(2*sigma_g) ): Xi := ln(1/(T0^2))+ln(1/Rp)+Lp: sys := [ D(Phip)(tau) = Phip(tau)*(Xi*y(tau)-ln(1/(T0^2))*y(tau)^(alpha*rho)-ln(1/Rp)-Lp+(Phip(tau)+N*Phis(tau))/(Phip(tau)+exp(G*(Phip(tau)+N*Phis(tau)))*N*Phis(tau))-1)/Xi, \ D(Phis)(tau) = Phis(tau)*(-ln(1/Rs)-Ls-kappa*ln(1/(T0^2))*y(tau)^(alpha*rho)+(Phip(tau)+N*Phis(tau))/(Phip(tau)*exp(-G*(Phip(tau)+N*Phis(tau)))+N*Phis(tau))-1)/Xi, \ D(y)(tau) = -gam*x*Phip(tau)/Xi]: DEplot(sys,[Phip(tau),Phis(tau),y(tau)],tau=0..n,[[Phip(0)=1e-7,Phis(0)=1e-7,y(0)=1]],stepsize=0.1,scene=[tau,Phip(tau)],method=classical[abmoulton],axes=FRAME,linecolor=BLUE): end: ODE_plot2 := proc(T0,Rp,Rs,Lp,Ls,n) gam := 1: rho := 1: kappa := 0.38: alpha := evalf( 5e-18/2.8e-20 ): x := 29:#  x = l[cav]/l[g] N := evalf(1033/1139): G := subs( {g=4.8e-9,c=3e10,h=6.62e-34,nu_p=3e10/1.033e-4,sigma_g=2.8e-20},g*c*h*nu_p/(2*sigma_g) ): Xi := ln(1/(T0^2))+ln(1/Rp)+Lp: sys := [ D(Phip)(tau) = Phip(tau)*(Xi*y(tau)-ln(1/(T0^2))*y(tau)^(alpha*rho)-ln(1/Rp)-Lp+(Phip(tau)+N*Phis(tau))/(Phip(tau)+exp(G*(Phip(tau)+N*Phis(tau)))*N*Phis(tau))-1)/Xi, \ D(Phis)(tau) = Phis(tau)*(-ln(1/Rs)-Ls-kappa*ln(1/(T0^2))*y(tau)^(alpha*rho)+(Phip(tau)+N*Phis(tau))/(Phip(tau)*exp(-G*(Phip(tau)+N*Phis(tau)))+N*Phis(tau))-1)/Xi, \ D(y)(tau) = -gam*x*Phip(tau)/Xi]: DEplot(sys,[Phip(tau),Phis(tau),y(tau)],tau=0..n,[[Phip(0)=1e-7,Phis(0)=1e-7,y(0)=1]],stepsize=0.1,scene=[tau,Phis(tau)],method=classical[abmoulton],axes=FRAME,linecolor=RED): end: ODE_plot3 := proc(T0,Rp,Rs,Lp,Ls,n) gam := 1: rho := 1: kappa := 0.38: alpha := evalf( 5e-18/2.8e-20 ): x := 29:#  x = l[cav]/l[g] N := evalf(1033/1139): G := subs( {g=4.8e-9,c=3e10,h=6.62e-34,nu_p=3e10/1.033e-4,sigma_g=2.8e-20},g*c*h*nu_p/(2*sigma_g) ): Xi := ln(1/(T0^2))+ln(1/Rp)+Lp: sys := [ D(Phip)(tau) = Phip(tau)*(Xi*y(tau)-ln(1/(T0^2))*y(tau)^(alpha*rho)-ln(1/Rp)-Lp+(Phip(tau)+N*Phis(tau))/(Phip(tau)+exp(G*(Phip(tau)+N*Phis(tau)))*N*Phis(tau))-1)/Xi, \ D(Phis)(tau) = Phis(tau)*(-ln(1/Rs)-Ls-kappa*ln(1/(T0^2))*y(tau)^(alpha*rho)+(Phip(tau)+N*Phis(tau))/(Phip(tau)*exp(-G*(Phip(tau)+N*Phis(tau)))+N*Phis(tau))-1)/Xi, \ D(y)(tau) = -gam*x*Phip(tau)/Xi]: DEplot(sys,[Phip(tau),Phis(tau),y(tau)],tau=0..n,[[Phip(0)=1e-7,Phis(0)=1e-7,y(0)=1]],stepsize=0.1,scene=[tau,y(tau)],method=classical[abmoulton],axes=FRAME,linecolor=BLACK): end:

Warning, `gam` is implicitly declared local to procedure `ODE_plot1`

Warning, `rho` is implicitly declared local to procedure `ODE_plot1`

Warning, `kappa` is implicitly declared local to procedure `ODE_plot1`

Warning, `alpha` is implicitly declared local to procedure `ODE_plot1`

Warning, `x` is implicitly declared local to procedure `ODE_plot1`

Warning, `N` is implicitly declared local to procedure `ODE_plot1`

Warning, `G` is implicitly declared local to procedure `ODE_plot1`

Warning, `Xi` is implicitly declared local to procedure `ODE_plot1`

Warning, `sys` is implicitly declared local to procedure `ODE_plot1`

Warning, `gam` is implicitly declared local to procedure `ODE_plot2`

Warning, `rho` is implicitly declared local to procedure `ODE_plot2`

Warning, `kappa` is implicitly declared local to procedure `ODE_plot2`

Warning, `alpha` is implicitly declared local to procedure `ODE_plot2`

Warning, `x` is implicitly declared local to procedure `ODE_plot2`

Warning, `N` is implicitly declared local to procedure `ODE_plot2`

Warning, `G` is implicitly declared local to procedure `ODE_plot2`

Warning, `Xi` is implicitly declared local to procedure `ODE_plot2`

Warning, `sys` is implicitly declared local to procedure `ODE_plot2`

Warning, `gam` is implicitly declared local to procedure `ODE_plot3`

Warning, `rho` is implicitly declared local to procedure `ODE_plot3`

Warning, `kappa` is implicitly declared local to procedure `ODE_plot3`

Warning, `alpha` is implicitly declared local to procedure `ODE_plot3`

Warning, `x` is implicitly declared local to procedure `ODE_plot3`

Warning, `N` is implicitly declared local to procedure `ODE_plot3`

Warning, `G` is implicitly declared local to procedure `ODE_plot3`

Warning, `Xi` is implicitly declared local to procedure `ODE_plot3`

Warning, `sys` is implicitly declared local to procedure `ODE_plot3`

>	display(ODE_plot1(0.65,0.75,0.96,5e-2,5e-2,300),axes=boxed,title=`laser photon density`); display(ODE_plot2(0.65,0.75,0.96,5e-2,5e-2,300),axes=boxed,title=`Stokes photon density`); display(ODE_plot3(0.65,0.75,0.96,5e-2,5e-2,300),axes=boxed,title=`relative inversion`);

So, we can obtain a quite efficient two-color pulsing in the nanosecond time domain without any additional wavelength conversion.

The nonstationary lasing in question produces the pulses with the durations, which are larger than the cavity period. The cavity length shortening, i. e. the use of the microchip lasers, decreases the pulse durations down to ten-hundred picoseconds. But there is a method allowing the fundamental pulse width reduction, viz. mode locking.

Part II. Conception of mode locking

The laser cavity is, in fact, interferometer, which supports the propagation of only defined light waves. Let consider a plane wave, which is reflected from a laser mirror. The initial wave is ( cc  is the complex conjugated term):

>	restart:  with(plots):   with(DEtools):    AI := A0*exp(I*(omega*t+k*x))+cc;

Warning, the name changecoords has been redefined

Then the reflected wave (normal incidence and full reflectivity are supposed) is:

>	AR := A0*exp(I*(omega*t-k*x+Pi))+cc;

where  is the phase shift due to reflection. An interference between incident and reflected waves results in

>	convert(AI+AR,trig):   expand(%): factor(%);

This is the so-called standing wave:

>	animate(sin(x)*sin(t), x=0..2*Pi,t=0..2*Pi, axes=boxed, title=`Standing wave`, color=red);

Note, that a wave node lies on a surface (point x =0). The similar situation takes place in the laser resonator. But the laser resonator consists of two (or more) mirrors and the standing wave is formed due to reflection from the each mirror. So, the wave in the resonator is the standing wave with the nodes placed on the mirrors. Such waves are called as the longitudinal laser modes . The laser resonator can contain a lot of modes with the different frequencies (but its nodes have to lie on the mirrors!) and these modes can interfere.

Let suppose, that the longitudinal modes are numbered by the index m . In fact, we have M  harmonic oscillators with the phase and frequency differences dphi  and domega , correspondingly. Let the amplitude of modes is A0 .

>	mode := 1/2*A0*exp(I*(phi0+m*dphi)+I*(omega0+m*domega)*t)+cc; # amplitude of the mode numbered by index m

Here 0 and 0 are the phase and frequency of the central mode, respectively.

The interference between these modes produces the wave packet:

>	packet := sum(mode-cc,m=-(M-1)/2..(M-1)/2)+cc;# interference of longitudinal modes with constant phases

Now, we can extract the term describing the fast oscillation on the central ("carrier") frequency 0 from the previous expression. The obtained result is the packet's envelope (its slowly varying amplitude):

>	envelope := expand((packet-cc)/exp(I*(t*omega0+phi0)));# slowly varying envelope of the wave packet

It is obviously, that this expression can be converted into following form:

>	envelope := 1/2*A0*sinh(1/2*I*M*(dphi+t*domega))/sinh(1/2*I*(dphi+t*domega));

The squared envelope's amplitude (i. e. a field intensity) is depicted in the next figure for the different M .

>	plot({subs({M=5,A0=1,dphi=0.1,domega=0.1},2*envelope^2),subs({M=20,A0=1,dphi=0.1,domega=0.1},2*envelope^2)},t=-80..80, axes=boxed, title=`result of modes interference`);

One can see, that the interference of modes results in the generation of short pulses. The interval between pulses is equal to 2 / domega .   The growth of  M  decreases the pulse duration 2 /( M*domega )   and to increases the pulse intensity * . The last is the consequence of the following relation:

>	(limit(sin(M/2*x)/sin(x/2),x=0)*A0)^2;# maximal field intensity

In this example, the phase difference between neighboring modes is constant. Such mode locking  causes the generation of short and intense pulses. But in the reality, the laser modes are not locked, i. e. the modes are the oscillations with the independent and accidental phases. In this case:

>

M := 20:# 20 longitudinal modes  A0 := 1:# constant dimensionless amplitude of mode   domega := 0.1:# constant frequency different    phi0 := 0:# phase of the central mode   omega0 := 1:# dimensionless frequency of central mode  mode := 0: for m from -(M-1)/2 to (M-1)/2 do:  die := rand(6):   dphi := die():# accidental phase difference between modes  mode := mode+A0*cos(phi0+m*dphi+(omega0+m*domega)*t): od: plot(mode,t=-80..80, axes=boxed, title=`result of modes interference`);

Thus, the interference of the unlocked modes produces the irregular field beatings, i. e.  the noise spikes with a duration ~ 1/(M*domega) .

What are the methods for the mode locking? Firstly, let consider the simplest model of the harmonic oscillation in the presence of the periodical force.

>	diff_equation := diff(diff(y(t),t),t)+omega^2*y(t)=cos(delta*t+phi);# oscillations in the presence of periodical force (delta and phi are the frequency and phase of the force oscillation, respectively)       dsolve({diff_equation, y(0)=1, D(y)(0)=0},y(t)):             oscill1 := combine(%);

>	animate(subs({omega=1,delta=0.1},subs(oscill1,y(t))),t=0..100,phi=0..2*Pi,color=red,style=point, axes=boxed);

We can see, that the external force causes the oscillations with the additional frequencies: +  and - . If  is equal to the intermode interval, the additional oscillation of mode plays a role of the resonance external force for the neighboring modes. Let consider such resonant oscillation in the presence of the resonant external force:

>	diff_equation := diff(diff(y(t),t),t)+omega^2*y(t)=cos(omega*t+phi);# resonant external force       dsolve({diff_equation,y(0)=1,D(y)(0)=0},y(t)):             oscill2 := combine(%);

The term, which is proportional to t  ("secular" term), equalizes the phase of the oscillations to the phase of the external force. It is the simplest model of a so-called active mode locking . Here the role of the external force can be played by the external amplitude or phase modulator. Main condition for this modulator is the equality of the modulation frequency to the intermode interval that causes the resonant interaction between modes and, as consequence, the mode locking (part 3 ).

The different mechanism, a passive mode locking , is produced by the nonlinear interaction of modes with an optical medium. Such nonlinearity can be caused by saturable absorption, self-focusing etc. (see further parts of worksheet). Now we shall consider the simplest model of the passive mode locking. Let suppose, that there are two modes, which oscillate with different phases and frequencies in the cubic nonlinear medium:

>

omega := 1:  delta := 0.1:   sys := [diff(diff(y(t),t),t)+omega^2*y(t)=-(y(t)^3+z(t)^3),\    diff(diff(z(t),t),t)+(omega+delta)^2*z(t)=-(y(t)^3+z(t)^3)];#two oscillating modes with nonlinear coupling  DEplot3d(sys,[y(t),z(t)],t=0..100,[[y(0)=-1,z(0)=1,D(y)(0)=0,D(z)(0)=0]],stepsize=1,scene=[y(t),z(t),t],axes=boxed,linecolor=BLACK,orientation=[-60,70], title=`mode locking`);

We can see, that the initial oscillations with the different phases are locked due to nonlinear interaction that produces the synchronous oscillations.

As conclusion, we note that the mode locking is resulted from the interference between standing waves with constant and locked phases. Such interference forms a train of the ultrashort pulses. The mechanisms of the mode locking are the external modulation with frequency, which is equal to intermode frequency interval, or the nonlinear response of the optical medium. Later on we shall consider both methods. But firstly we have to obtain the more realistic equations describing the ultrashort pulse generation.

Part III. Basic model

The models describing the laser field evolution are based usually on the so-called semi-classical approximation. In the framework of this approximation the field obeys the classical Maxwell equation and the medium evolution has the quantum-mechanical character. Here we shall consider the wave equation without concretization of the medium evolution.

The Maxwell equation for the light wave can be written as:

>	restart:  with(PDEtools,dchange): maxwell_eq := diff(E(z,t),z,z)-diff(E(z,t),t,t)/c^2=4*Pi*diff(P(t),t,t)/c^2;

where E(z,t)  is the field strength, P(t)  is the medium polarization, z  is the longitudinal coordinate, t  is the time, c  is the light velocity. The change of the variables   z --> z* , t - z/c --> t*  produces

>	macro(zs=`z*`,ts=`t*`): tr := {z = zs, t = ts + zs/c}; maxwell_m := dchange(tr,maxwell_eq,[zs,ts],simplify);

We does not take into account the effects connected with a wave propagation in thin medium layer, that allows to eliminate the second-order derivation on z* . Then

>	int(op(2,expand(lhs(maxwell_m))),ts)-int(rhs(maxwell_m),ts):#integration of both sides of wave equation     numer(%):          wave_1 := expand(%/(-2));

So, we reduced the order of wave equation. The inverse transformation of the coordinates leads to the so-called shortened wave equation:

>	tr := {zs = z, ts = t - z/c}; wave_2 := dchange(tr,wave_1,[z,t],simplify);

Next step is the transition to the slowly-varying amplitude approximation. We shall consider field envelope (z,t) and polarization P(t), which are filled by the fast oscillation with frequency  ( k  is the wave number, N is the concentration of the atoms, d  is the medium length):

>	Theta=omega*t-k*z;# phase E := rho(z,t)*cos(Theta);# field P := N*d*(a(t)*cos(Theta)-b(t)*sin(Theta));# polarization (a and b are the quadrature components)

Then the wave equation can be transformed as:

>

Theta:=omega*t-k*z:   diff(E,z)+diff(E,t)/c+2*Pi*diff(P,t)/c:     combine(%,trig):       collect(%,cos(Theta)):     wave_3 := collect(%,sin(Theta)):   eq_field := op(1,expand(coeff(wave_3,cos(Theta))))+op(3,expand(coeff(wave_3,cos(Theta))))=-op(2,expand(coeff(wave_3,cos(Theta))))-op(4,expand(coeff(wave_3,cos(Theta))));# field equation disp_conditions := Int(coeff(wave_3,sin(Theta)),t)=0;# dispersion condition

The obtained result is the system of the shortened wave equation and the dispersion condition. The right-hand side of the field equation (material part) will be different for the different applications (see below).

Part IV. Active mode locking: harmonic oscillator

Amplitude modulation

We start our consideration with  a relatively simple technique named as active mode locking due to amplitude modulation. The modulator, which is governed by external signal and changes the intracavity loss periodically, plays the role of the external "force" (see part 1). Let consider the situation, when the modulation period is equal and the pulse duration is much less than the cavity period. If the modulation curve is close to cosine, then the master equation for the ultrashort pulse evolution can be written as [ D. Kuizenga, A. Siegman,"FM and AM mode locking of homogeneous laser", IEEE J. Quant. Electr., 6 , 694 (1970) ]:

 = g  + ,

where g  is the net-gain in the laser accounting for the saturated gain  and linear loss l  (including output loss),   is the inverse bandwidth of the spectral filter, M  is the characteristic of the modulation strength taking into account curvature of the modulation curve at the point of maximal loss.

Let try to solve this equation in the case of the steady-state propagation of ultrashort pulse (when =0) and in the absence of detuning of modulation period from cavity period. If the time is normalized to  , then the obtained equation is a well-known equation of harmonic oscillator:

>	restart:  with(plots):   with('linalg'):    ode[1] := diff(rho(t),`$`(t,2)) + g*rho(t) - M*t^2*rho(t);

Warning, the name changecoords has been redefined

Warning, the protected names norm and trace have been redefined and unprotected

>	sol := dsolve(ode[1]=0,rho(t));

The next step is suggested by the asymptotic behavior of the prospective solution:  = 0. But previously it is convenient to transit from Whittaker  functions to hypergeometric  and Kummer  functions:

WhittakerM( , ,z)=  hypergeom( , 1+2 , x ),

 WhittakerW( , ,z)=  KummerU( , 1+2 , x ).

As result, we have ( = , = , x = ):

(t)=   t    hypergeom( , , ) +

  t    KummerU( , , )=

  t    hypergeom(( )+ , , ) +

    hypergeom( , , ).

Now, the asymptotic condition  = 0, which is similar to condition for quantum states in harmonic oscillator, results in (see, for example,  S. Flugge, Practical Quantum Mechanics I, Springer-Verlag (1971) )

(t) =  C ( t )     for n= - ( ) && n is integer,

(t) =  C ( t )     for n= - ( ) && n is integer,

where  and  are Hermite  polynomials. Value of constant C  can be obtained from the energy balance condition, which results from the equation of gain saturation:

=  .

Here  is the inverse flux of the gain saturation energy,  is the gain for small signal defined by gain medium properties and pump intensity (note, that g  = ).

Now let investigate the parameters of the steady-state solution of the master equation.

With that end in view, we have to search the generation field energy:

>	C^2*int( HermiteH(0,t*surd(M, 4))^2 * exp( -t^2*sqrt(M) ), t=-infinity..infinity ); #energy of the first solution ("ground state"), where alpha=l+sqrt(M) (see above)

The use of normalization of intensity  to

and energy balance condition (see above) results in

>	eq := l+sqrt(M) - alpha[0]/(1+C^2*sqrt(Pi)/surd(M, 4)) = 0:  sol := solve(eq, C^2);# pulse peak intensity

Hence the ultrashort pulse obtained as result of active mode locking can be represented as

>	pulse := sol*HermiteH(0,t*surd(M, 4))^2 * exp( -t^2*sqrt(M) );# this is first solution of master equation, which represents a Gaussian pulse profile  plot( subs({l=0.1,M=0.05,alpha[0]=1.2},pulse), t=-6..6, axes=BOXED, title=`field intensity vs time` );

Now we calculate the pulse duration measured on half-level of maximal pulse intensity:

>	eq := exp( -sqrt(M)*t^2 ) = 1/2:  sol := solve(eq, t):   pulse_width := simplify( sol[1] - sol[2] );

>	plot( pulse_width,M=0.3..0.01, axes=BOXED, title=`pulse width vs M` );

One can see, that the increase of modulation parameter decreases the pulse width, but this decrease is slow (~1/ ).

Next step is the taking into account the detuning of cavity and modulation periods. The corresponding normalized parameter  can be introduced in the following form (see corresponding Doppler transformation: t-->t - z  and, as consequence, -->-  for steady-state pulse, here  is, in fact, inverse relative velocity with dimension [time/cavity transit], i. e. time delay on the cavity round-trip):

>	ode[2] := diff(rho(t),`$`(t,2)) + delta*diff(rho(t),t) + g*rho(t) - M*t^2*rho(t);

>	sol := dsolve(ode[2]=0,rho(t));

The comparison with above obtained result leads to

(t) =  C ( t )     for n= - ( ) && n is integer,

(t) =  C ( t )     for n= - ( ) && n is integer,

and we can repeat our previous analysis

>	C^2*int( HermiteH(0,t*surd(M, 4))^2 * exp( -t*(delta+sqrt(M)*t) ), t=-infinity..infinity ); # energy of the first solution ("ground state"), where alpha=l+sqrt(M)+delta^2/4

>	eq := l+sqrt(M)+delta^2/4 - alpha[0]/(1+%) = 0:# energy balance condition  sol_0 := solve(eq, C^2);# pulse peak intensity

Note, that there is the maximal  permitting the "ground state" ultrashort pulse generation:

>	-4*l-4*sqrt(M)-delta^2+4*alpha[0] > 0;  solve(-4*l-4*sqrt(M)-delta^2+4*alpha[0]=0, delta);

We see, that the upper permitted level of detuning parameter is increased due to pump growth (rise of ) and is decreased by M  growth.

It is of interest, that the growth of  does not leads to generation of the "excite state" pulses because of the corresponding limitation for these solutions is more strict:

>

En := C^2*int( expand( HermiteH(1,t*surd(M, 4) ))^2 * exp( -t*(delta+sqrt(M)*t) ), t=-infinity..infinity ):# energy of the second solution ("excite state")  eq1 := 0 = - (3/4+(delta^2-4*(alpha-l))/(16*sqrt(M))):   solve(%,alpha):    eq := % - alpha[0]/(1+En) = 0:# energy balance condition     solve(eq, C^2);# "excite state" pulse peak intensity      -12*sqrt(M)-4*l-delta^2+4*alpha[0] > 0;# condition of generation

The dependence of pulse width on  has following form:

>	eq := exp( -t*(delta+sqrt(M)*t) ) = 1/2:  sol := solve(eq, t):   pulse_width := simplify( sol[1] - sol[2] );    plot( subs(M=0.05,pulse_width),delta=-5..5, axes=BOXED, title=`pulse width vs delta` );

We can see almost linear and symmetric rise of pulse width due to detuning increase. This is an important characteristic of active mode-locked lasers. But in practice (see below), the dependence of the pulse parameters on detuning has more complicated character.

The dependence of the pulse maximum location on detuning parameter can be obtained from the solution for pulse envelope's maximum:  = 0. Hence, the pulse maximum location is:

>	diff(exp(-t*(delta+sqrt(M)*t)/2), t) = 0:  solve(%, t);

>	plot( subs(M=0.05,%),delta=-5..5, axes=BOXED, title=`pulse shift vs delta` );

The increase (decrease) of the pulse round-trip frequency ( <0 and >0, respectively) increases positive (negative) time shift of the pulse maximum relatively modulation curve extremum.

Phase modulation

The external modulation can change not only field amplitude, but its phase. In fact, this regime (phase active modulation) causes the Doppler frequency shift of all field components with the exception of those, which are located in the vicinity of extremum of modulation curve. Hence, there exists the steady-state generation only for field located in vicinity of points, where the phase is stationary. Steady-state regime is described by equation

>	ode[3] := diff(rho(t),`$`(t,2)) + g*rho(t) + I*phi*rho(t) + I*M*t^2*rho(t);

where  is the phase delay on the cavity round-trip.

This equation looks like previous one, but has complex character. This suggests to search its partial solution in the form ( t ) =  :

>

expand( subs(rho(t)=C*exp((a+I*b)*t^2), ode[3]) ):  expand(%/C/exp(t^2*a)/exp(I*t^2*b)):   eq1 := collect (coeff(%,I), t^2 ):    eq2 := collect( coeff(%%,I,0), t^2 ):     eq3 := coeff(eq1, t^2);      eq4 := coeff(eq1, t,0);       eq5 := coeff(eq2, t^2);        eq6 := coeff(eq2, t,0);         sol1 := solve(eq6=0, a);# solution for a        sol2 := solve(subs(a=sol1,eq3=0), b);# solution for b       sol3 := solve(subs(b=sol2,eq4)=0, phi);# solution for phi # but NB that there is eq5 defining, in fact, g:      solve( subs({a=sol1, b=sol2},eq5)=0, g);# solution for g

So, we have ( t ) =    and  =  , g= . The time-dependent parabolic phase of ultrashort pulse is called chirp and is caused by phase modulation. Pulse width for this pulse is:    

>	exp(-sqrt(M/2)*t^2) = 1/2:  sol := solve(%, t):   pulse_width := simplify( sol[1] - sol[2] );

that is 2  . As one can see, that result is equal to one for amplitude modulation. The main difference is the appearance of the chirp.

Let take into account the modulation detuning. In this case we may to suppose the modification of the steady-state solution:

>

ode[4] := diff(rho(t),`$`(t,2)) + g*rho(t) + delta*diff(rho(t),t) + I*phi*rho(t) + I*M*t^2*rho(t);     expand( subs(rho(t)=C*exp((a+I*b)*t^2 + (c+I*d)*t), ode[4]) ):# c is the shift of the pulse profile from extremum of the modulation curve, d is the frequency shift of the pulse from the center of the gain band     expand(%/C/exp(t^2*a)/exp(I*t^2*b)/exp(t*c)/exp(I*t*d)): eq1 := collect (coeff(%,I), t ):  eq2 := collect( coeff(%%,I,0), t ):   eq3 := coeff(eq1, t^2):    eq4 := coeff(eq1, t):     eq5 := coeff(eq1, t, 0):      eq6 := coeff(eq2, t^2);       eq7 := coeff(eq2, t):        eq8 := coeff(eq2, t, 0): allvalues( solve(subs(b=-a,{eq3=0,eq4=0,eq5=0,eq7=0,eq8=0}),{a,c,d,g,phi}) );

We see, that there is not frequency shift of the pulse ( d= 0), but, as it was for amplitude modulation, the time delay appears ( ) that changes the pulse duration as result of modulation detuning. The rise of the detuning prevents from the pulse generation due to saturated net-gain coefficient increase ( ). The pulse parameters behavior coincides with one for amplitude modulation.

Ultrashort pulse stability

Now we shall investigate the ultrashort pulse stability against low perturbation ( t ). The substitution of the perturbed steady-state solution in dynamical equation with subsequent linearization on ( t ) results in

 =  +   - l  +  - ,

where  is the perturbed saturated gain, which is obtained from the assumption about small contribution of perturbation to gain saturation process:   

>	alpha[p] = op(2,convert( series( alpha[0]/(1 + A + B), B=0,2), polynom));

Here A= , B=  and we neglected the high-order terms relatively perturbation amplitude. Note, that this is negative quantity.

Let the dependence of perturbation on z  has exponential form with increment . Then

>	zeta(z,t) := upsilon(t)*exp(lambda*z):# perturbation  ode[4] := expand( diff(zeta(z,t),z)/exp(lambda*z) ) =  expand( (alpha*zeta(z,t) + alpha[p]*rho(t) - l*zeta(z,t) + diff(zeta(z,t),`$`(t,2)) - M*t^2*zeta(z,t))/exp(lambda*z) );

Now we introduce * where B*= , that allows to eliminate the exponent from right-hand side of  (we shall eliminate the asterix below).  

This is equation for eigenvalues  and eigenfunctions  of the perturbed laser operator. Stable generation of the ultrashort pulse corresponds to decaying of perturbations, i.e. <0. For Gaussian pulse:

>	ode[4] := upsilon(t)*lambda = subs(rho(t)=sqrt(surd(M,4)*(alpha[0]-l-sqrt(M))/(sqrt(Pi)*(l+sqrt(M))))*exp(-t^2*sqrt(M)/2), alpha*upsilon(t)+alpha[p]*rho(t)-l*upsilon(t)+diff(upsilon(t),`$`(t,2))-M*t^2*upsilon(t)); sol := dsolve( ode[4],upsilon(t) );

As result (see above), we have

(t) =    + C ( t )     for n= - ( ) && n is integer,

(t) =    + C ( t )     for n= - ( ) && n is integer,

>

eq1 := n = - (1/4-(alpha-lambda-l)/(4*sqrt(M))):  eq2 := n = - (3/4-(alpha-lambda-l)/(4*sqrt(M))):   En := subs( C^2=-surd(M,4)*(l+sqrt(M)-alpha[0])/(sqrt(Pi)*(l+sqrt(M))),  C^2*sqrt(Pi)/(M^(1/4)) ):# ultrashort pulse energy (see above)    lambda[1,n] = solve(eq1, lambda);# increment     lambda[2,n] = solve(eq2, lambda);# increment

The glance on the solution is evidence of absence of solution corresponding to . For others modes we take into account, that = l + . As result, all perturbation modes are unstable because of

    < 0 ( n= 1, 2, ...),

 < 0 ( n= 0, 1, ...).

It should be noted, that for the amplified pulse

,

,

Here we see the decrease of increment as result of n  rise. Therefore, only "ground state"  with  will be amplified predominantly. This fact provides for Gaussian pulse generation (see, F.X. Kartner, D.M. Zumbuhl, and N. Matuschek, "Turbulence in Mode-Locked Lasers", Phys. Rev. Lett. 82 , 4428 (1999)).

In the presence of detuning we have:

>	eq1 := n= - (1/4+(delta^2-4*(alpha-lambda-l))/(16*sqrt(M))):  eq2 := n= - (3/4+(delta^2-4*(alpha-lambda-l))/(16*sqrt(M))):   lambda[1,n] = simplify( solve(eq1, lambda) );    lambda[2,n] = simplify( solve(eq2, lambda) );

It is surprisingly, but, as it was above, our linear analysis predicts the pulse stability regardless of detuning  (because of = l + ). But for the pulse the detuning growth decreases the increment that does not favor the single pulse generation.

Nonlinear processes: self-phase modulation and dynamical gain saturation

Among above considered effects only gain saturation by full pulse energy can be considered as nonlinear process, which, however, does not affect on the pulse envelope, but governs its energy. The time-dependent nonlinear effects, which can transform pulse profile, are self-phase modulation (SPM) and dynamical gain saturation. First one is the dependence of the field phase on its intensity and can play essential role in solid-state lasers [ V.L. Kalashnikov, I.G. Poloyko, V.P. Mikhailov, "Phase modulation of radiation of solid-state lasers in the presence of Kerr optical nonlinearity", Optics and Spectriscopy, 84 , 104 (1998) ]. Second effect is caused by the change of the gain along pulse profile and is essential in lasers with large gain cross-sections and comparatively narrow gain band [ V.L. Kalashnikov, I.G. Poloyko, V.P. Mikhailov, "Generation of ultrashort pulses in lasers with external frequency modulation", Quantum Electronics, 28 , 264 (1998) ].

At first, let analyze the presence of SPM, which can be considered as perturbation of our master equation describing active phase modulation:

>	ode[5] := diff(rho(t),`$`(t,2)) + g*rho(t) + I*phi*rho(t) + I*M*t^2*rho(t) - I*beta*C^2*exp(2*a*t^2)*rho(t);# phase perturbation in the form -I*beta*C^2*abs(rho(t))^2*rho(t), rho(t) has unperturbed profile

>

expand( subs(rho(t)=C*exp((a+I*b)*t^2), ode[5])/exp(a*t^2)/exp(I*t^2*b)/C ):  series(%,t=0,3):# approximate solution as result of expansion   convert(%, polynom):# now we shall collect the coefficients of t and t^2 for real and imaginary parts    eq1 := coeff(%,I):     eq2 := coeff(%%,I,0):      eq3 := coeff(eq1,t^2);       eq4 := coeff(eq1,t,0);        eq5 := coeff(eq2,t^2);         eq6 := coeff(eq2,t,0);        sol1 := solve(eq6=0, a);# solution for a       sol2 := solve(subs(a=sol1,eq3=0), b);# solution for b      sol3 := solve(subs(b=sol2,eq4)=0, phi);# solution for phi # but NB that there is eq5 defining, in fact, g:     solve( subs({a=sol1, b=sol2},eq5)=0, g);# solution for g

Hence we have solution with perturbed parameters:

( t ) =    and   g=  .

As one can see, this solution has enlarged chirp and reduced pulse duration due to SPM.

Now we shall investigate the influence of dynamical gain saturation on the pulse characteristics in the case of active mode locking due to amplitude modulation. Let the contribution of the dynamical gain saturation can be considered as perturbation for Gaussian pulse (see above). The instant energy flux of such pulse is:

>	int(sol_0*exp(-t*(delta+sqrt(M)*t)/2),t):  Energy := simplify( coeff( %, erf(1/2*sqrt(2)*M^(1/4)*t+1/4*delta*sqrt(2)/(M^(1/4))))*(1+erf(1/2*sqrt(2)*M^(1/4)*t+1/4*delta*sqrt(2)/(M^(1/4))) ) );

The approximation of the small contribution of the gain saturation allows the expansion of energy into series on t  up to second order:

>	series(Energy,t=0,2):  convert(%, polynom):   simplify(%);

Hence we have the modified gain coefficient ( t ), where  is the gain coefficient at pulse peak,  is the pulse intensity for unperturbed solution. Note, that the additional term in brackets is resulted from the shift of pulse maximum.

Then the master equation for perturbed solution (in our approximation!) is

>	ode[6] := diff(rho(t),`$`(t,2)) + g*rho(t) + delta*diff(rho(t),t) - (M*t^2+epsilon*t)*rho(t);# here epsilon is alpha*I[0]

>	dsolve(ode[6]=0, rho(t));

The comparison with above obtained result gives

(t) =  C ( t )     for n= - ( ) && n is integer,

(t) =  C ( t )     for n= - ( ) && n is integer.

In future we shall omit -terms.
If we take the unperturbed intensity for calculation of perturbation action, the perturbed pulse energy is

>

C^2*HermiteH(0, t*sqrt(sqrt(M))+epsilon/(2*sqrt(M*sqrt(M))) )^2* exp(-(4*M^2*t^2+4*M*t*epsilon)/(4*M*sqrt(M))-delta*t):# field intensity profile for "ground state" (1/4-(4*g*M-M*delta^2)/(16*M*sqrt(M))) = 0   subs(epsilon=alpha*I[0],%):     int(%,t=-infinity..infinity):        En := simplify(%);# energy of the first "ground state"

>	1/4-(4*(alpha-l)*M - M*delta^2)/(16*M*sqrt(M)):  sol_a := solve(%=0,alpha):   eq := sol_a - alpha[0]/(1+subs( alpha=%,En )) = 0:# energy balance condition    sol_I := solve(eq, C^2);# pulse peak intensity

Now we plot the dependencies of the pulse intensity, width and maximum location versus detuning parameter .

>	plot(subs( {alpha[0]=1.2,l=0.1,M=0.05},subs(I[0]=sol_0,sol_I) ), delta=-2..1.5, axes=BOXED, title=`pulse intensity vs detuning`);

>

eq := exp( -(4*M^2*t^2+4*M*t*epsilon)/(4*M*sqrt(M))-delta*t ) = 1/2:# we take into account pulse profile without epsilon^2  sol := solve(eq, t):   pulse_width := simplify( subs( epsilon^2=0,sol[1] - sol[2] ) ):    subs(epsilon=alpha*I[0],pulse_width):     subs({alpha=sol_a,I[0]=sol_0},%): # pulse width for perturbed solution

>	plot(subs( {alpha[0]=1.2,l=0.1,M=0.05},% ), delta=-2..1.5, axes=BOXED, title=`pulse width vs detuning`);

>	# pulse maximum location diff(exp(-(2*M*t+epsilon)^2/(8*M*sqrt(M))-delta*t/2), t) = 0:  solve(%, t);

>	subs(epsilon=alpha*I[0],%):  subs({alpha=sol_a,I[0]=sol_0},%):   plot(subs( {alpha[0]=1.2,l=0.1,M=0.05},% ), delta=-2..1.5, axes=BOXED, title=`pulse location vs detuning`);

We see, that the main peculiarity here is the asymmetric dependence of the pulse parameters on . The pulse width minimum and intensity maximum don't coincide with =0 and the detuning characteristics have sharper behavior in negative domain of detuning.

Now, as it was made in previous subsection, we estimate the condition of the ultrashort pulse stability.

Here we take into consideration -term, but this does not fail our analysis because of this term contribute only to pulse energy without shift pulse inside modulation window.

For the sake of the simplification, we shall consider the contribution of the destabilizing field to dynamical gain saturation, but only in the form of the unperturbed peak intensity variation  and perturbation of the saturated gain coefficient (see above). Then the stability condition:

>

eq1 := epsilon = subs(I[0]=sol_0,sol_a*(I[0]+zeta)):# zeta is the peak intensity variation -(3/4-(4*(sol_a-lambda-l)*M+epsilon^2-M*delta^2)/(16*M*sqrt(M))):  expand( solve(%=0, lambda) ) < 0;# perturbation increment from eigenvalue problem   subs(epsilon=rhs(eq1),lhs(%)):    animate3d(subs( {alpha[p]=0,alpha[0]=1.2,l=0.1},% ), delta=-2..1.5,  M=0.005..0.05, zeta=-0.15..0.15, axes=BOXED, title=`maximal increment of perturbation`);

We can see, that the perturbation rise can destabilize the pulse as result of    increase (compare with subsection "Amplitude modulation" ). Also, there is the possibility of ultrashort pulse destabilization near =0. So, the presence of dynamical gain saturation gives the behavior of the ultrashort pulse parameters and stability condition, which is close to the experimentally observed and numerically obtained (see, for example, J.M. Dudley, C.M. Loh and J.D. Harvey, Stable and unstable operation of a mode-locked argon laser, Quantum and Semiclass. Opt. , 8 , 1029 ( 1996 )).

Now we try to investigate the influence of dynamical gain saturation in detail by using so-called aberrationless approximation . Let the pulse profile keeps its form with accuracy up to n -order of time-series expansion, but n  pulse parameters are modified as result of pulse propagation. Then the substitution of the expression for pulse profile in master equation with subsequent expansion in t -series produces the system of n  ODE describing the evolution of pulse parameters.

>

f1 := (z,t)->rho0(z)*exp(-a(z)^2*t^2+b(z)*t);# approximate pulse profile      f2 := (z,tau)->rho0(z)*exp(-a(z)^2*tau^2+b(z)*tau):      master := diff(rho(z,t),z) = (alpha - l)*rho(z,t) - chi*alpha*rho(z,t)*int(rho(z,tau)^2,tau=-infinity..t) + diff(rho(z,t),t$2) - M*t^2*rho(z,t) + delta*diff(rho(z,t),t);# approximate master equation with dynamical gain saturation (parameter chi), alpha is the gain coefficient before pulse front expand(lhs(subs({rho(z,t)=f1(z,t),rho(z,tau)=f2(z,tau)},master))*exp(a(z)^2*t^2)/exp(t*b(z)) - rhs(subs({rho(z,t)=f1(z,t),rho(z,tau)=f2(z,tau)},master))*exp(a(z)^2*t^2)/exp(t*b(z))):# substitution of the approximate solution convert( series(%,t=1/2*b(z)/(a(z)^2),3),polynom ):# expansion around  peak at t=1/2*b(z)/(a(z)^2)   eq1 := collect(%,t):    eq2 := subs({diff(rho0(z),z)=u,diff(a(z),z)=v,diff(b(z),z)=w},coeff(eq1,t^2)):     eq3 := subs({diff(rho0(z),z)=u,diff(a(z),z)=v,diff(b(z),z)=w},coeff(eq1,t)):      eq4 := subs({diff(rho0(z),z)=u,diff(a(z),z)=v,diff(b(z),z)=w},coeff(eq1,t,0)):       sol := simplify (solve({eq2=0,eq3=0,eq4=0},{u,v,w}) );

Now try to find the steady-state points of ODE-system, which correspond to stationary pulse parameters. For this aim let introduce substitution  :

>	sol_main := {w = -chi*alpha*rho0(z)^2-2*delta*a(z)^2-4*b(z)*a(z)^2, v = 1/2*(M-4*a(z)^4)/a(z), u = 1/4*rho0(z)*(-chi*alpha*rho0(z)^2*sqrt(Pi)*sqrt(2)*a(z)+4*alpha*a(z)^2-8*a(z)^4-4*l*a(z)^2+2*rho0(z)^2*b(z)*chi*alpha+4*b(z)*delta*a(z)^2+4*b(z)^2*a(z)^2)/(a(z)^2)};

>	eq1 := subs( {rho0(z)^2=x,a(z)=sqrt(y),a(z)^2=y,a(z)^4=y^2},expand( numer( subs(sol_main,u) )/rho0(z) )=0 );# here we eliminated the trivial solution   eq2 := subs( a(z)^4=y^2,numer( subs(sol_main,v) )=0 );     eq3 :=  subs( {rho0(z)^2=x,a(z)^2=y},subs(sol_main,w)=0 );

>	allvalues( solve({eq2,eq3},{y,b(z)}) );

Hence ,  ( x  is the pulse intensity). So, the shift is   that differs from the usual result (see above) as result of dynamical gain saturation (last term), which shifts the pulse maximum in negative side. This additional shift has obvious explanation. The gain at the pulse front is greater than one at pulse tail due to gain saturation. This shifts the pulse forward as hole .

Pulse width is:

>	eq := exp( t*(-delta/2-chi*x*alpha/(2*sqrt(M))) - sqrt(M)*t^2/2 ) = 1/2:  sol := solve(eq, t):   pulse_width := simplify( sol[1] - sol[2] );

We can see, that there is the minimum of the pulse duration in negative domain of  that corresponds to result, which was obtained on the basis of perturbation theory. The pulse intensity :

>	Intensity := simplify( solve( subs({y = 1/2*sqrt(M),b(z) = -1/2*(delta*M+chi*sqrt(M)*x*alpha)/M},eq1), x)[1]);# pulse intensity

So, we have the following dependencies for the pulse duration

>	subs(x=Intensity,pulse_width): animate(subs( {alpha=1.2,l=0.1},% ), delta=-4..4, M=0.01..0.1, numpoints=200, axes=BOXED, title=`pulse width vs detuning`);

and pulse intensity

>	animate(subs( {alpha=1.2,l=0.1,chi=1},Intensity ), delta=-4..4, M=0.01..0.1, numpoints=200, axes=BOXED, view=0..0.45, title=`pulse intensity vs detuning`);

The last step is the stability analysis. The stability of our solutions can be estimated from the eigenvalues of Jacobian of sol .

>

eq1 := subs({rho0(z)=x,a(z)=y,b(z)=z},subs(sol_main,u)):  eq2 := subs({rho0(z)=x,a(z)=y,b(z)=z},subs(sol_main,v)):   eq3 := subs({rho0(z)=x,a(z)=y,b(z)=z},subs(sol_main,w)): m[1,1] := diff( eq1,x ):  m[1,2] := diff( eq1,y ):    m[1,3] := diff( eq1,z ): m[2,1] := diff( eq2,x ):  m[2,2] := diff( eq2,y ):    m[2,3] := diff( eq2,z ): m[3,1] := diff( eq3,x ):  m[3,2] := diff( eq3,y ):    m[3,3] := diff( eq3,z ): A := array([[m[1,1],m[1,2],m[1,3]],[m[2,1],m[2,2],m[2,3]],[m[3,1],m[3,2],m[3,3]]]):# Jacobian

Now we find the eigenvalues  of Jacobian directly by calculation of determinant.

>	evalm(A-[[lambda,0,0],[0,lambda,0],[0,0,lambda]]): numer( simplify(det(%)) ): sol := solve(%=0,lambda):

>

x := sqrt(Intensity):  y := sqrt( sqrt(M)/2 ):   z := -delta/2-x^2*alpha/(2*sqrt(M)):    s1 := evalf( subs({chi=1,l=0.1,alpha=1.2},sol[1]) ):     s2 := evalf( subs({chi=1,l=0.1,alpha=1.2},sol[2]) ):      s3 := evalf( subs({chi=1,l=0.1,alpha=1.2},sol[3]) ): plot3d({Re(s1),Re(s2),Re(s3)}, delta=-4..2, M=0.01..0.1, axes=boxed,title=`Re(lambda) for initial perturbation`);

One can see that the pulse with Gaussian-like form is stable in the region of its existence (see two previous figures). We have to note that the considered here perturbations belong to limited class therefore this criterion is necessary but not sufficient condition of pulse stability (compare with previous consideration on the basis of perturbation theory, where we analyzed not only pulse peak variation but also gain coefficient change).

Part V. Nonlinear Schroedinger equation: construction of the soliton solution using the direct Hirota's method

The Schroedinger equation is the well-known nonlinear equation describing the weak nonlinear waves, in the particular, the laser pulse propagation in fibers. In the last case, a pulse can propagate without decaying over large distance due to  balance between two factors: SPM and group delay dispersion (GDD). These pulses are named as optical solitons [ M. J. Ablowitz, H. Segur, "Solitons and the Inverse Scattering Transform", SIAM Philadelphia, 1980 ]. The ultrashort pulse evolution obeys to the next master equation:

which is the consequence of eq_field  from part 2  in the case of transition to local time t-->t-z/c . The right-hand side terms describe GDD (with coefficient ) and SPM (with coefficient ), respectively .

It is very important to obtain the exact soliton solutions of nonlinear equations. There are the inverse and direct methods to obtain such solutions. One of the direct methods is the so-called Hirota's method. The main steps of this method are: 1) the selection of the suitable substitution instead of the function     (see the master equation),   that allows to obtain the bilinear form of the evolution equation; 2) the consideration of the formal series of perturbation theory for this bilinear equation. In the case of soliton solutions these series are terminated.

The useful substitution for the nonlinear Schroedinger equation is (z,t)=G(z,t)/F(z,t). Let suppose that  F  is the real function. It should be noted that we can make any assumption about   to satisfy the assumptions 1) and 2). Hirota proposed to introduce a new D - operator   in following way:

After substitution of   in the terms of functions G  and F  we obtain two bilinear differential equations with regard to the new operator D :

                      (1)

The functions G  and F  can be expanded into the series of the formal parameter :

  

Let substitute G  and F  into Eq. (1) and treat the terms with powers of  as independent, to get the infinite set of the differential equations relatively G1, G3, ...; F2, F4, ...  . These formal series are terminated only in the case when the master equation has exact N -soliton solution. For instance, the set of first six differential equations in our case is:

For sake of the simplification of the very cumbersome manipulations we introduce the procedure for operator , which acts on the functions a  and b . The lasts are the exponents (or linear combination of the exponents) in the form , where  is linear function.

>	restart: with(plots):

Warning, the name changecoords has been redefined

Procedure  

>

Dt_Dz := proc (a,b,m,n)     local Summa,k,r,result:     Summa := 0:     if  (n>1) and (m<>0) then         for k from 1 to n-1 do             Summa := Summa+binomial(n,k)*(-1)^(n-k+m)*             der( der(b,z,(n-k)),t,m)*der(a,z,k)+             binomial(n,k)*(-1)^(n-k)*der(b,z,(n-k))*             der( der(a,z,k),t,m)         od:     fi:     if  (n>1) and (m>1) then        for r from 1 to (m-1) do            for k from 1 to (n-1) do                Summa := Summa+binomial(m,r)*                binomial(n,k)*(-1)^(n-k+m-r)*                der( der(b,z,(n-k)),t,(m-r))*                der( der(a,z,k),t,r);            od:        od:     fi:     if  (m>1) and (n<>0) then        for r from 1 to (m-1) do            Summa := Summa+binomial(m,r)*(-1)^(m-r+n)*            der( der(b,z,n),t,(m-r))*der(a,t,r)+            binomial(m,r)*(-1)^(m-r)*der(b,t,(m-r))*            der( der(a,z,n),t,r);        od;     fi:     if (m<>0) and (n<>0) then        Summa := Summa+(-1)^(m+n)*der(der(b,z,n),t,m)*a+        (-1)^m*der(a,z,n)*der(b,t,m)+(-1)^n*der(a,t,m)*        der(b,z,n)+der(der(a,z,n),t,m)*b;       fi:        if (m=0) and (n>1) then        Summa := Summa+(-1)^(n)*der(b,z,n)*a+der(a,z,n)*b;            for k from 1 to (n-1) do                Summa := Summa+binomial(n,k)*(-1)^(n-k)*                der(b,z,(n-k))*der(a,z,k);            od:     fi:     if (m=0) and (n=1) then        Summa := der(a,z,1)*b-der(b,z,1)*a:     fi:     if (n=0) and (m>1) then        Summa := Summa+(-1)^(m)*der(b,t,m)*a+der(a,t,m)*b;            for r from 1 to (m-1) do                Summa := Summa+binomial(m,r)*(-1)^(m-r)*                der(b,t,(m-r))*der(a,t,r);            od:     fi:       if (n=0) and (m=1) then        Summa := der(a,t,1)*b-der(b,t,1)*a:     fi:     if (n=0) and (m=0) then        Summa := a*b     fi:     result := combine(Summa,exp):  end:

The next procedure will be used for calculation of the derivative of  (or combination of exponents) on t  or z   with further simplification of the obtained expression.

Procedure der

>

der := proc (f,x,m)    local  difF,i,function:    subs(eta1=eta1(x),eta1s=eta1s(x),    eta2=eta2(x),eta2s=eta2s(x),f):    difF := diff(%,x$m):    if (x=t) then       function := subs({diff(eta1(x),x)=b1,       diff(eta2(x),x)=b2,diff(eta1s(x),x)=b1s,       diff(eta2s(x),x)=b2s},difF)             else         function := subs({diff(eta1(x),x)=a1,       diff(eta2(x),x)=a2,diff(eta1s(x),x)=a1s,       diff(eta2s(x),x)=a2s},difF)    fi;    subs(eta1(x)=eta1,eta1s(x)=eta1s,    eta2(x)=eta2,eta2s(x)=eta2s, function):    if (m>1) then       difF := subs({diff(a1,x)=0,diff(a2,x)=0,       diff(b1,x)=0,diff(b2,x)=0,diff(a1s,x)=0,       diff(a2s,x)=0,diff(b1s,x)=0,diff(b2s,x)=0},%)             else       combine(%)    fi:    collect(%,exp): end:

The next procedure is used to calculate an integral of  (or combination of exponents) on t  or z with further simplification of the expression.

Procedure Integr

>

integr := proc (f,x,m)      local intF,i,g1,g1s,g2,g2s,function:      intF := subs(eta1=g1*x,eta1s=g1s*x,eta2=g2*x,      eta2s=g2s*x,f):      for i from 1 to m do          intF := int(intF,x);      od:      if (x=t) then          x := t; g1 := b1; g1s := b1s; g2 := b2; g2s := b2s;               else          x := z; g1 := a1; g1s := a1s; g2 := a2; g2s := a2s;      fi:      intF;      collect(%,exp):      subs(b1*t=eta1,b1s*t=eta1s,b2*t=eta2,      b2s*t=eta2s,a1*z=eta1,a1s*z=eta1s,a2*z=eta2,      a2s*z=eta2s,%); end:

Now, let try to obtain a first-order soliton  for nonlinear Schroedinger equation.

>	macro(Gs=`G*`,G1s=`G1*`,G3s=`G3*`,G5s=`G5*`, as=`a*`,bs=`b*`,eta0s = `eta0*`):

>	G1 := exp(eta1):# successful substitution! I*der(%,z,1)+k_2*der(%,t,2): #first from the equations set    factor(%);

>	#as result we can find the parameter a1    a1 := I*k_2*b1^2:       a1s := -I*k_2*b1s^2:

>	G1s := exp(eta1s):# conjugated to G1    G1G1s := combine(G1*G1s):        F2 := beta/(2*k_2)*integr(%,t,2);#F2 from the second equation of the set

>	# But the next equation of the set results in    I_Dz_G1_F2 := I*factor(Dt_Dz(G1,F2,0,1)):        d_Dt2_G1_F2 := k_2*factor(Dt_Dz(G1,F2,2,0)):            factor(I_Dz_G1_F2+d_Dt2_G1_F2);                Dt_Dz(F2,F2,2,0);

To obtain this result we use the trivial relationships:

>	eta1=a1*z+b1*t+eta10:    eta2=a2*z+b2*t+eta20:        Dz*exp(eta1)*exp(eta2)=(a1-a2)*exp(eta1+eta2);            Dt^2*exp(eta1)*exp(eta2)=(b1-b2)^2*exp(eta1+eta2);

As was shown above a= - i , hence the last term in third equation of set is equal to 0 . So, we are to choose G3  = 0 to satisfy third equation. Furthermore F2  from fourth equation of the set is ( ) exp( + _s). But in concordance with above obtained relationships this expression is equal to zero. So we can choose F4  = 0. Thus to satisfy other equations we can keep in the expansion of functions G  and F  only G1 , G3  and F2 .

So, the formal series are terminated. Since  is independent parameter we can take =1.

>	rho := G1/(1+F2);

>	subs(eta1=a1*z+b1*t+eta10,eta1s=a1s*z+b1s*t+eta10s, rho): soliton := expand(subs({b1s=b1,beta=1,k_2=1/2},%));#the choice of k_2 and beta is only normalization of the values in equation

>	# Now we check the obtained solution by substitution of one in dynamical equation I*der(rho,z,1)/rho+beta*exp(eta1+eta1s)/((1+1/2*beta*exp(eta1+eta1s)/(k_2*(b1+b1s)^2)))^2+k_2*der(rho,t,2)/rho: simplify(%);

All right! This is the exact solution of the Schroedinger equation. Physically b1  has a sense of the inverse pulse duration. So it is real parameter. But what is the free parameter 10 ? Let 10 is real and  = b1. Then