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

>	subs({exp(eta10)=2*b1, exp(eta10s)=2*b1},soliton): simplify(%);

The obtained solution is the so-called first-order soliton  with duration  , amplitude b1  and phase z/2 :

  (z,t)=b1*sech(b1*t)exp(i* z/2)

>	plot(subs({z=0,b1=1},%),t=-5..5,axes=boxed, title=`first-order soliton`);

But what is about different values of 10 ? The choice of the different values of 10  results in the obtaining of the so-called collapsing pulses, i. e. pulses with singularity in the dependence of their parameters on z .

The described here procedure is available for the analysis of the higher-order solitons. For example, the substitution of G1 = +  causes the termination of the series on fifth equation (you can prove this statement by using of the described above procedures). The obtained solution depends on z  and is called as the second order soliton .

>	4*exp(-I*z/2)(cosh(3*t)+3*exp(-4*I*z)*cosh(t))/(cosh(4*t)+4*cosh(t)+3*cos(4*z)):# second-order soliton    plot3d(abs(%)^2,t=-2..2,z=0..5,color=red,axes=boxed, title=`second-order soliton`);

In the case of Schroedinger equation, the main feature of the described above method is the termination of the formal series for an arbitrary order of the solution. Such behavior results from a very rich mathematical structure of the dynamical equation: an existence of the infinitely many nontrivial symmetries and conservation laws, Painleve property and integrability by means of the inverse scattering method. The non-decaying pulse-like solutions of the integrable nonlinear evolutional equations are called as the solitons . But as we shall see later, some nonlinear equations have the soliton-like solutions, but do not belong to integrable class. For example, for nonlinear Landau-Ginzburg equation (see next part) there is the soliton-like solution in the first order of the Hirota's method. But the second-order solution does not lead to the termination of the series. Such soliton-like solutions of the nonintegrable dynamical equation are called as the quasi-solitons  (or solitary waves).

Part VI. Nonlinear Landau-Ginzburg equation: quasi-soliton solution

Here we shall consider a soliton-like pulse, which is generated in the continuous-wave solid-state laser due to power-dependent saturation of the diffraction loss in the presence of the field self-focusing in active medium [ H.A. Haus, "Theory of mode locking with a fast saturable absorber", J. Appl. Phys., 46 , 3049 (1975); H.A. Haus, J.G. Fujimoto, and E.P. Ippen, "Analytic theory of additive pulse and Kerr lens mode locking", IEEE J. Quantum Electronics, 28 , 2086 (1992) ]. The action of the saturable loss can be described by the real cubic nonlinear term. The energy dissipation due to spectral filtering can be introduced by means of the real second-order derivative on t . Then in the absence of SPM and GDD, the dynamical equation is a analog of the Schroedinger equation, but with pure real terms.

 = g  +  +   

Here g  is the net-gain in the laser taking into account the gain and linear loss in the active medium, output loss, and diffraction loss. For sake of simplification, we shall suppose the normalization of time to the inverse bandwidth of the spectral filter  (let = 2.5 fs, that corresponds to the full generation bandwidth of Ti: sapphire laser) and normalization of pulse intensity to the inverse intensity of the loss saturation  (the typical values of  are ~  -    /W). We will consider the steady-state pulse propagation, then  = 0. So, the master equation can be transformed to the well-known Duffing's type equation describing nonlinear oscillations without damping:

>	restart:   with(plots):     with(DEtools):        ode := diff(rho(t),`$`(t,2)) + rho(t)^3 + g*rho(t);

Warning, the name changecoords has been redefined

Its implicit solutions are:

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

The first integral of motion is:

>	numer(diff(lhs(sol[1]),t)):    int_motion := simplify((op(1,%)^2-op(2,%)^2)/2);

These equations describe the motion in the potential:

>	pot := simplify(op(2,int_motion)+op(3,int_motion));

The value of 4*_C1  plays a role of the full energy of system. The dependence of potential on  for the different g  is shown in the next figure:

>	plot3d(subs(rho(t)=rho,pot),g=0.05..-0.1,rho=-0.5..0.5,axes=boxed,title=`Potential of pendulum`);

We can see, that for g  > 0 there is one equilibrium state of pendulum in  = 0 (stable state), and for g  < 0 there are three one (unstable in  = 0 and stable in  = +/- ) . Obviously, that in this system the different oscillation regimes is possible, that is illustrated by the phase portrait on the plane [ y,   z=d /dt ] .

>	sys := convertsys(ode = 0,[],rho(t),t,z); dfieldplot([diff(rho(t),t)=z(t),diff(z(t),t)=-rho(t)^3-subs(g=-0.1,g)*rho(t)], \ [z(t),rho(t)],t=-2..2,rho=-0.5..0.5,z=-0.1..0.1,\ arrows=LARGE,axes=boxed,title=`Nonlinear oscillations`,color=black);

The quasi-soliton solution of the initial equation corresponds to the oscillation around the stable equilibrium state with infinite period. In this case the full energy is equal to 0. Then   _C1  = 0 and the motion begins from  = 0 at t  = - :

>	plot({sqrt(subs({rho(t)=rho,g=-0.1},-pot))/2,-sqrt(subs({rho(t)=rho,g=-0.1},-pot))/2},\ rho=0..0.45,axes=boxed,labels=[`rho(t)`,`drho(t)/dt`],color=red,numpoints=200);

The amplitude of quasi-soliton can be found from the integral of motion (i. e. pulse maximum correspond to d /dt = 0 and _C1 = 0 in the integral of motion).

>	rho0 := solve(factor(pot)/rho(t)^2=0,rho(t));

The explicit integration of the solution produces:

>	assume(g<0):   sol1_a := numer(value(subs(_C1=0,lhs(sol[1]))));     sol1_b := numer(value(subs(_C1=0,lhs(sol[2]))));

Make some transformations:

>	sol2_a := numer(simplify(expand(sol1_a/op(1,sol1_a)),radical,symbolic)); sol2_b := numer(simplify(expand(sol1_b/op(1,sol1_b)),radical,symbolic));

When at the pulse maximum (0)  = , the value of _C2  can be found as:

>	i_C2_a := solve(subs({t=0,rho(t)=rho_max},sol2_a)=0,_C2); i_C2_b := solve(subs({t=0,rho(t)=rho_max},sol2_b)=0,_C2);

Let tanh ( i_C2_a*sqrt ( -g ))  = tanh ( -i_C2_b*sqrt ( -g ) )  = -/+ , then from an expression for the tangents of sum of arguments: tanh ( a + b )  = tanh ( a )  + tanh ( b ) /(1 + tanh ( a ) tanh ( b ) ), the equations can be transformed as:

>	sol3_a := simplify(tanh(op(1,sol2_a)))=(tanh(-op(2,sol2_a))+upsilon)/\ (1+tanh(-op(2,sol2_a))*upsilon); sol3_b := simplify(tanh(op(1,sol2_b)))=(tanh(-op(2,sol2_b))+upsilon)/\ (1+tanh(-op(2,sol2_b))*upsilon);

>	sol4_a := solve(sol3_a, rho(t)); sol4_b := solve(sol3_b, rho(t));

Now, we must note, that the transit  -> 0  (see above) corresponds to the transit  --> . Then the final solutions result from the next operations:

>	sol_fin_1 := limit(sol4_a[1],upsilon=infinity);  sol_fin_2 := limit(sol4_a[2],upsilon=infinity);    sol_fin_1 := limit(sol4_b[1],upsilon=infinity);      sol_fin_2 := limit(sol4_b[2],upsilon=infinity);

The positive roots satisfy to the initial condition and the result is the quasi-soliton pulse with sech  - shape envelope, which has the duration 1/  and the peak amplitude 0 = . The pulse intensity profile is shown in the next figure (  =   /W,  is the time normalized to ):

>	animate(evalf(subs(t=iota/(2.5e-15),sol_fin_1)^2*100),\ iota=-1e-13..1e-13,g=-0.05..-0.01,frames=50,\ axes=boxed,color=red,labels=[`time, fs`,`rho^2, GW/cm^2`],title=`Pulse envelope`);

So, there is the close analogy between an ultrashort pulse generation in the laser with power-dependent loss saturation and an oscillation of nonlinear pendulum. The considered solution has the form of a steady-state quasi-soliton. But in the next part we are going to demonstrate some analogies of the high-order soliton behavior (laser breezers).

Part VII. Autooscillations of quasi-solitons in the laser

The main efforts in the design of the modern femtosecond lasers aim to the stabilization of the ultrashort pulse parameters. For example, as it was shown [ A. M. Sergeev, E. V. Vanin, F. W. Wise, Stability of passively mode-locked lasers with fast saturable absorbers, Opt. Commun., 140, 61 (1997) ], the pulse destabilizations are very important factors limiting the operation of laser with so-called fast saturable absorber (the bleaching of diffraction loss due to Kerr-lensing is the example of such saturable absorber, see previous part). The destabilization of the laser quasi-soliton can produce its destruction or regular/nonregular autooscillations. The last can be considered as the analog of a high-order soliton formation (see part 4), which is reminiscent of the breezers of nonlinear dynamical equation.

We shall consider the master laser equation, which joins the Landau-Ginzburg and Schroedinger equations. The main nonlinear factors now are SPM and power-dependent loss saturation. GDD and spectral filtering will be taken into consideration, too.

>	restart:   with(DEtools):     with(plots):        master_1 := diff(rho(z,t),z) = alpha*rho(z,t)-gamma*rho(z,t)+I*phi*rho(z,t)+tf^2*diff(rho(z,t),`$`(t,2))+I*k_2*diff(rho(z,t),`$`(t,2))+sigma*rho(z,t)^2*conjugate(rho(z,t))-I*beta*rho(z,t)^2*conjugate(rho(z,t));

Warning, the name changecoords has been redefined

Here  is the phase delay on the full cavity round trip,  and  are the gain and loss coefficients, respectively. The general exact solution of this equation is not known, but there is the quasi-soliton solution in the following form:

>	f1 := (z,t)->rho0(z)*sech(t*tau(z))^(1+I*psi(z));

Here 0 is the pulse amplitude,  is the inverse pulse width,  is the chirp. This solution obeys the condition of steady-state propagation, when  = 0, i.e. the pulse parameters are constant. To describe the nonstationary pulse propagation, we shall use the aberrationless approximation: the changes of the pulse parameters do not cause the aberration of the pulse form. Next step is the substitution of f1  into master_1  with following expansion in t  - series. The coefficients of the expansion produce the set of ODE for the evaluating pulse parameters.

>

assume(tau(z),real):   assume(t,real):   expand(lhs(subs(rho(z,t)=f1(z,t),master_1))-rhs(subs(rho(z,t)=f1(z,t),master_1))):   eq := subs({alpha=alpha(z),diff(rho0(z),z)=a,diff(tau(z),z)=b,diff(psi(z),z)=c},\ convert(series(%,t=0,3),polynom)):   assume(rho0(z),real): eq1 := evalc(coeff(eq,t^2)):   eq2 := evalc(coeff(eq,t)):     eq3 := evalc(coeff(eq,t,0)):        eq4 := coeff(eq1,I):           eq5 := coeff(eq1,I,0):              eq6 := coeff(eq3,I):                 eq7 := coeff(eq3,I,0): sol := solve({eq4=0,eq5=0,eq6=0,eq7=0},{a,b,c,phi}): sys := diff(rho0(z),z)=subs(sol,a),diff(tau(z),z)=subs(sol,b),diff(psi(z),z)=subs(sol,c):

The obtained system sys  have to be supplemented with equation for gain evolution:

>	sys := {%,diff(alpha(z),z)=alpha(z)*exp(-2*xi*rho0(z)^2/tau(z)-1/Tr-Pump)+Pump*alphamx*(1-exp(-1/Tr-Pump))/(Pump+1/Tr)-alpha(z)};

Here Tr  is the gain relaxation time normalized to the cavity period, Pump is the dimensionless pump (see part 8 ),  is the inverse gain saturation energy.

At first let's find the parameters of a steady-state quasi-soliton, viz. the solution of sys independent on  z.

>	f[1] :=  -2*tf^2*psi*Tau-4*k_2*psi^2*Tau-4*k_2*Tau-2*tf^2*psi^3*Tau-2*sigma*Phi*psi-2*beta*Phi;# Tau = tau^2, Phi = rho^2 f[2] := alpha+sigma*Phi-gamma-tf^2*Tau+k_2*psi*Tau; f[3] := 3*k_2*Tau*psi+tf^2*psi^2*Tau-2*tf^2*Tau+sigma*Phi;

>	sol1 := solve({f[2]=0,f[3]=0},{Tau,Phi});  subs({Tau=subs(sol1,Tau),Phi=subs(sol1,Phi)},f[1]):   simplify(%):    numer(%)/2/(-alpha+gamma):     sol2 := solve(%=0,psi);

Normalization of the time to tf  and the intensity to  allows simple plotting of the obtained result:  

>	plot({subs( {beta=1,tf=1,sigma=1},sol2[1] ),subs( {beta=1,tf=1,sigma=1},sol2[2] )},k_2=-100..10,axes=boxed,view=-10..1, title=`chirp vs. GDD`);

To plot the pulse duration and its intensity we have to take into account the dependence of  on the pulse energy. However for the sake of simplification we suppose =const<  (see previous section). Then  

>	plot(subs( {beta=1,tf=1,sigma=1,gamma=0.05, alpha=0.04},subs(psi=sol2[2],1/sqrt(subs(sol1,Tau))) ),k_2=-100..10,axes=boxed, title=`pulse width vs. GDD`);

Only solution with >0 in the region of the negative GDD and <0 for the positive GDD (green curve) has a physical meaning because it corresponds to the positive square of the pulse width 1/ . The existence of the pulse duration minimum in the vicinity of zero GDD pays a very important role in the pulse shortening technique based on the creation of appropriate negative GDD in the laser cavity.

Now, we shall consider the evolution of ultrashort pulse parameters on the basis of the obtained system of ODE. Also, this allows to analyse the ultrashort pulse stability in the framework of our model. We suppose to solve the system by the standard operator DEplot. Let's normalize  and  to 1,7*   /W, times to  (2,5 fs for Ti: sapphire laser), then  = 0.0018. The fundamental step is the assumption about saturation of the Kerr nonlinearity:

 = (1 -  ),

 =   (1 -  ),

where  and  are the unsaturated nonlinear parameters.

>

#procedure for solution of the obtained system of ODE ODE_plot := proc(alphamx,gam,sigma0,beta0,Tr,Pump,xi,disp,tg,n)  sigma := sigma0*(1 - sigma0*rho0(z)^2/2 ):   beta := beta0*(1 - beta0*rho0(z)^2/2 ):    sys := [D(alpha)(z) = alpha(z)*exp(-2*xi*rho0(z)^2/tau(z)-1/Tr-Pump)+Pump*alphamx*(1-exp(-1/Tr-Pump))/(Pump+1/Tr)-alpha(z), \ D(psi)(z) = -4*disp*tau(z)^2-2*tg^2*psi(z)*tau(z)^2-4*disp*psi(z)^2*tau(z)^2-2*beta*rho0(z)^2-2*beta*rho0(z)^2*psi(z)^2-2*tg^2*psi(z)^3*tau(z)^2, \ D(rho0)(z) = sigma*rho0(z)^3-gam*rho0(z)-tg^2*rho0(z)*tau(z)^2+alpha(z)*rho0(z)+disp*rho0(z)*psi(z)*tau(z)^2, \ D(tau)(z) = -2*tg^2*tau(z)^3+sigma*rho0(z)^2*tau(z)+3*disp*tau(z)^3*psi(z)+beta*rho0(z)^2*psi(z)*tau(z)+tg^2*psi(z)^2*tau(z)^3]: DEplot(sys,[rho0(z),tau(z),psi(z),alpha(z)],z=0..n,[[rho0(0)=0.001,tau(0)=0.01,alpha(0)=0,psi(0)=0]],stepsize=1,scene=[z,rho0(z)],method=classical[foreuler],axes=FRAME,linecolor=BLACK): end:

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

Warning, `beta` is implicitly declared local to procedure `ODE_plot`

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

The next figure demonstrates the autooscillations of pulse amplitude (quasi-breezer behavior).

>	display(ODE_plot(0.5,0.05,10,1,300,0.0004,0.0018,-10,1,15000));

The character of the pulse evolution strongly depends on the correlations between system's parameters. For example, next figure demonstrates the pulse stabilization due to negative dispersion growth.

>	display(ODE_plot(0.5,0.05,10,1,300,0.0004,0.0018,-20,1,15000));

But the pumping growth produces the irregular autooscillations.

>	display(ODE_plot(0.5,0.05,10,1,300,0.00047,0.0018,-20,1,5000));

The next parameter's set produces the chaotic oscillations:

>	display(ODE_plot(0.5,0.05,1,0,300,0.0004,0.0018,0,1,8000));

The pulse parameters (amplitude and inverse pulse duration) corresponding to irregular oscillations can be shown on the basis of the iteration procedure, which realizes the direct Euler method for the solution of the described above ODE system.

>

# Attention! This block can take a lot of CPU-time iterations := proc(iter)  alphamx := 0.5:    gam := 0.05:      sigma0 := 1:        beta0 := 0:          Tr := 300:             Pump := 0.0004:          xi := 0.0018:        disp := 0:      tg := 1:    rho0n:=0.001;    taun:=0.01;     alphan:=0;       psin:=0;  for m from 1 to iter do     rho0old:=rho0n;           tauold:=taun;        psiold:=psin;        alphaold:=alphan;     sigma := evalhf(sigma0*(1 - sigma0*rho0old^2/2 )):        beta := evalhf(beta0*(1 - beta0*rho0old^2/2 )):     alphan := evalhf( alphaold*exp(-2*xi*rho0old^2/tauold-1/Tr-Pump)+Pump*alphamx*(1-exp(-1/Tr-Pump))/(Pump+1/Tr));     psin := evalhf(psiold-4*disp*tauold^2-2*tg^2*psiold*tauold^2-4*disp*psiold^2*tauold^2-2*beta*rho0old^2-2*beta*rho0old^2*psiold^2-2*tg^2*psiold^3*tauold^2);     rho0n := evalhf(rho0old+sigma*rho0old^3-gam*rho0old-tg^2*rho0old*tauold^2+alphaold*rho0old+disp*rho0old*psiold*tauold^2);     taun := evalhf(tauold-2*tg^2*tauold^3+sigma*rho0old^2*tauold+3*disp*tauold^3*psiold+beta*rho0old^2*psiold*tauold+tg^2*psiold^2*tauold^3); if m = iter then pts := [rho0n, taun] fi;     od;     pts  end: PLOT(seq(POINTS(iterations(i)), i=9800 .. 10000), SYMBOL(POINT));

Warning, `alphamx` is implicitly declared local to procedure `iterations`

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

Warning, `sigma0` is implicitly declared local to procedure `iterations`

Warning, `beta0` is implicitly declared local to procedure `iterations`

Warning, `Tr` is implicitly declared local to procedure `iterations`

Warning, `Pump` is implicitly declared local to procedure `iterations`

Warning, `xi` is implicitly declared local to procedure `iterations`

Warning, `disp` is implicitly declared local to procedure `iterations`

Warning, `tg` is implicitly declared local to procedure `iterations`

Warning, `rho0n` is implicitly declared local to procedure `iterations`

Warning, `taun` is implicitly declared local to procedure `iterations`

Warning, `alphan` is implicitly declared local to procedure `iterations`

Warning, `psin` is implicitly declared local to procedure `iterations`

Warning, `m` is implicitly declared local to procedure `iterations`

Warning, `rho0old` is implicitly declared local to procedure `iterations`

Warning, `tauold` is implicitly declared local to procedure `iterations`

Warning, `psiold` is implicitly declared local to procedure `iterations`

Warning, `alphaold` is implicitly declared local to procedure `iterations`

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

Warning, `beta` is implicitly declared local to procedure `iterations`

Warning, `pts` is implicitly declared local to procedure `iterations`

This is the so-called strange attractor, i. e. the chaotic attracting manifold: the pulse parameters are changed chaotically but within the limited range.

Thus, the pulse oscillations in the Kerr-lens mode-locked laser were analyzed (see J.Jasapara, V.L.Kalashnikov, D.O.Krimer, I.G.Polyko, M.Lenzner, W.Rudolph, "Automodulations in cw Kerr-lens mode-locked solid-state lasers", J. Opt. Soc. Am. B, 17, 319 (2000) , and arXiv: physics/0009020). Such oscillations accompany the negative dispersion decrease and the pump growth and close connect with the nonlinear factors in the system. The regular oscillation is the analogue of the high-order soliton propagation and the irregular one is the analogue of the nonlinear breezer.

Part VIII. Numerical approaches: ultrashort pulse spectrum, stability and multipulsing

Kerr-lens mode-locked Cr:LiSGaF-laser with the Raman self-scattering in active medium

The aberrationless approximation demonstrates the stability loss in the vicinity of zero dispersions. It is necessary to interpret the phenomenon. Moreover, we did not consider a lot of the lasing factors, which affect the ultrashort pulse dynamics in sub-100 fs region, viz. higher-order dispersion, stimulated Raman self-scattering (see Part I), strong fast absorber saturation etc. It is clear, that in order to take into account these phenomena we need the numerical simulations beyond the computer algebra abilities. Hovewer, Maple can help in the preparation of the source code and in the interpretation of the obtained results.

Now we describe the simplest generation model, which is highly useful for the numerical simulations. This model is based on the generalized Landau-Ginzburg equation:

>

restart:  with(codegen,fortran):   with('linalg'):    with(stats):     with(plots):   master1 := diff(a(z,t),z) = (alpha(z) - rho - gamma/(1+sigma*Phi(z,t)))*a(z,t) + diff(a(z,t),t$2) + sum('(-I)^(k+1)*D[k]*diff(a(z,t),t$k)','k'=2..N) - I*Phi(z,t)*a(z,t); # Phi is the field intensity, gamma is the fast absorber modulation depth, sigma is the inverse saturation intensity, we used the normalization of the field to the self-phase modulation coefficient and the time to the inverse gain bandwidth, D[k] are the dispersion coefficients absorbed the (1/factorial(k)) factors

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

Warning, the name changecoords has been redefined

For the gain evolution we have:

>	master2 := diff(alpha(z),z) = P*(a - alpha(z)) -alpha(z)*Int(Phi(z,t),t=-T[cav]/2..T[cav]/2)/E[s] - T[cav]*alpha(z)/T[r];

Here a  is the maximal gain for the full population inversion, P  is the dimensionless pump ( Pump  from the previous section),  and  are the cavity period and the gain relaxation time, respectively,  is the gain saturation energy taking into account the accepted normalization of time and intensity. The best methods for the solution of the system ( master1 , master2 ) is the split-step Fourier method. Then in the Fourier domain we escape to calculate the partial derivatives:

>

op(2,rhs(master1));  inttrans[fourier](%, t, omega)*F(omega); print(`F(w) takes into account the shape factors for the gainband and output coupler spectral profiles`);    subs( N=6,op(3,rhs(master1)) );# we confine the maximal dispersion order     inttrans[fourier](evalc( % ), t, omega):       factor(%);

To define the parameters of the simulation we use the experimental data for Cr:LiSGaF-laser kindly presented by Dr. I.T. Sorokina and Dr. E. Sorokin and written in the corresponding *.txt-files. NOTE! ALL DATA FILES HAVE TO BE PLACED IN YOUR CURRENT DIRECTORY:

>	currentdir();

We used the next experimental setup, which is typical for the Kerr-lens mode-locked lasers (high reflective plane and spherical (focal length f=5 cm) mirrors HR, chirped mirrors Ch, output mirror OC, distances are shown by the numbers in the figure):

The first step is the analysis of the gain band, which allows to define the gain bandwidth and the parameters of the time normalization ( tf  from the previous section) and the frequency normalization (this is the gain bandwidth).

>

gain_x := readdata(`gain_x.dat`,1,float):# wavelength  gain_y := readdata(`gain_y.dat`,1,float):# gain cross-section      n := vectdim(gain_y):     plot([[gain_x[k],gain_y[k]] $k=1..n],color=red,view=0..3e-20,axes=boxed,title=`gain band profile`);    max_cross_section := max(gain_y[k] $k=1..n);#gain band maximum   half_cross_section := evalf(max_cross_section/2);#half of the gain band maximum  P := array(1..n): Q := array(1..n):  for i from 2 to n do   P[i] := evalf(2*Pi*3*1e10/(gain_x[i]*1e-7)):#transition from wavelength to frequency    Q[i] := [P[i],gain_y[i]]:     if gain_y[i]=max_cross_section then X_max := P[i] else fi:    if gain_y[i]>half_cross_section and gain_y[i-1]<half_cross_section then X_half_1 :=  evalf((P[i-1]+P[i])/2) else fi:   if gain_y[i]<half_cross_section and gain_y[i-1]>half_cross_section then X_half_2 := evalf((P[i-1]+P[i])/2) else fi:  od: print(`position of the gain maximum:`); X_max; print(`position of the first half of maximum:`);  X_half_1; print(`position of the second half of maximum:`);   X_half_2; print(`bandwidth:`);    bandwidth := evalf(abs(X_half_2-X_half_1)); print(`inverse bandwidth for Gaussian approximation [s]:`);     minimal_pulse_width:=evalf(4*ln(2)/bandwidth); for i from 1 to n do  P[i] := evalf((2*Pi*3*1e10/(gain_x[i]*1e-7)-X_max)/bandwidth):   Q[i] := [P[i],gain_y[i]]: od:    j := 'j':      plot([Q[j] $j=1..n],axes=BOXED, color=red, title=`gain cross section versus normalized frequency`,view=0..3e-20,color=blue);

The similar manipulations are performed for the output coupler. The fit-approximation gives the spectral shape factor for the output mirror.

>

Out_x := readdata(`out_x.dat`,1,float):  Out_y := readdata(`out_y.dat`,1,float):   n := vectdim(Out_y):    plot([[Out_x[k],Out_y[k]] $k=1..n],color=red,axes=boxed,title=`experimental reflection profile`);     P0 := array(1..n):      Q0 := array(1..n):       for i from 1 to n do      P0[i] := evalf((2*Pi*3*1e10/(Out_x[i]*1e-7)-X_max)/bandwidth):#transition from wavelength to frequency    Q0[i] := evalf(1-Out_y[i]):  od: f1 := plot([[P0[k],Q0[k]] $k=1..n],color=red):#experimental output loss profile in the dimensionless frequency domain  eq0 := fit[leastsquare[[x,y], y=a*x^6+b*x^5+c*x^4+d*x^3+e*x^2+f*x+g]]([[P0[k] $k=1..n], [Q0[k] $k=1..n]]);# fit-approximation   Q0 := [subs(x=P0[k],rhs(eq0)) $k=1..n]:     f2 := plot([[P0[k],Q0[k]] $k=1..n],color=blue):      display(f1,f2,axes=boxed,title=`comparison of the experimental and approximated profiles`);

The first picture presents the experimental data, the second results from the fit-approximation ( eq0 ) and the transition to the dimensionless frequencies .

Now let's find the group-delay dispersion (GDD) in the active medium from the measured data:

>

#crystal length is 0.8 cm for double pass GDD_cry_x := readdata(`gdd_crystal_x.dat`,1,float):  GDD_cry_y := readdata(`gdd_crystal_y.dat`,1,float):   n := vectdim(GDD_cry_y):    fig1 := plot([[GDD_cry_x[k],0.8*GDD_cry_y[k]] $k=1..n],color=red,axes=boxed):     display(fig1,title=`measured GDD in active crystal`);      P1 := array(1..n):       Q1 := array(1..n):      for i from 1 to n do     P1[i] := evalf((2*Pi*3*1e10/(GDD_cry_x[i]*1e-7)-X_max)/bandwidth):#transition from wavelength to frequency   od:  f1 := plot([[P1[k],0.8*(bandwidth*10^(-15))**2*GDD_cry_y[k]]\ $k=1..n],color=red): eq1 := fit[leastsquare[[x,y], y=a*x^2+b*x+c]]([[P1[k] $k=1..n], evalf(0.8*(bandwidth*10^(-15))**2)*GDD_cry_y]);# fit-approximation  Q1 := [subs(x=P1[k],rhs(eq1)) $k=1..n]:   f2 := plot([[P1[k],Q1[k]] $k=1..n],color=blue,axes=boxed):    display(f1,f2,title=`GDD fit-approximation in frequency domain`);

eq1  is the result of the fit-approximation.

The main devices for the GDD manipulation are the chirp mirrors. Their parameters are:

>

#Ch1 (double pass) GDD_Ch1_x := readdata(`gdd_ch1_x.dat`,1,float):  GDD_Ch1_y := readdata(`gdd_ch1_y.dat`,1,float):   n := vectdim(GDD_Ch1_y):    plot([[GDD_Ch1_x[k],2*GDD_Ch1_y[k]] $k=1..n],color=red,axes=boxed,title=`measured GDD for Ch1`);      P2 := array(1..n-4):       Q2 := array(1..n-4):      for i from 1 to n-4 do     P2[i] := evalf((2*Pi*3*1e10/(GDD_Ch1_x[i]*1e-7)-X_max)/bandwidth):#transition from wavelength to frequency    Q2[i] := 2*(bandwidth*10^(-15))**2*GDD_Ch1_y[i] od:   f1 := plot([[P2[k],Q2[k]] $k=1..n-4],color=red):  eq2 := fit[leastsquare[[x,y], y=a*x^6+b*x^5+c*x^4+d*x^3+e*x^2+f*x+g]]([[P2[k] $k=1..n-4], [Q2[k] $k=1..n-4]]); Q2 := [subs(x=P2[k],rhs(eq2)) $k=1..n-4]:  f2 := plot([[P2[k],Q2[k]] $k=1..n-4],color=blue):   display(f1,f2,axes=boxed,title=`fit-approximation in frequency domain`);

>

#Ch2 (four passes) GDD_Ch2_x := readdata(`gdd_ch2_x.dat`,1,float):  GDD_Ch2_y := readdata(`gdd_ch2_y.dat`,1,float):   n := vectdim(GDD_Ch2_y):    plot([[GDD_Ch2_x[k],4*GDD_Ch2_y[k]] $k=1..n],color=red,axes=boxed,title=`measured GDD for Ch2`);      P3 := array(1..n):       Q3 := array(1..n):        for i from 1 to n do       P3[i] := evalf((2*Pi*3*1e10/(GDD_Ch2_x[i]*1e-7)-X_max)/bandwidth):#transition from wavelength to frequency      od:     f1 := plot([[P3[k],4*(bandwidth*10^(-15))**2*GDD_Ch2_y[k]] $k=1..n],color=red):    eq3 := fit[leastsquare[[x,y],\ y=a*x^6+b*x^5+c*x^4+d*x^3+e*x^2+f*x+g]]([[P3[k] $k=1..n],\ evalf(4*(bandwidth*10^(-15))**2)*GDD_Ch2_y]);   Q3 := [subs(x=P3[k],rhs(eq3)) $k=1..n]:  f2 := plot([[P3[k],Q3[k]] $k=1..n],color=blue): display(f1,f2,axes=boxed,title=`fit-approximation in frequency domain`);

And, at the last, the GDD in output coupler and high-reflective mirrors:

>

#Oc (single pass)  GDD_Oc_x := readdata(`gdd_oc_x.dat`,1,float):   GDD_Oc_y := readdata(`gdd_oc_y.dat`,1,float):    n := vectdim(GDD_Oc_y):     plot([[GDD_Oc_x[k],GDD_Oc_y[k]] $k=1..n],color=red,axes=boxed,title=`measured GDD for OC`);       P4 := array(1..n):        Q4 := array(1..n):       for i from 1 to n do      P4[i] := evalf((2*Pi*3*1e10/(GDD_Oc_x[i]*1e-7)-X_max)/bandwidth):#transition from wavelength to frequency     od:    f1 := plot([[P4[k],GDD_Oc_y[k]*(bandwidth*10^(-15))**2] $k=1..n],color=red):   eq4 := fit[leastsquare[[x,y], y=a*x^6+b*x^5+c*x^4+d*x^3+e*x^2+f*x+g]]([[P4[k] $k=1..n], GDD_Oc_y*evalf((bandwidth*10^(-15))**2)]);  Q4 := [subs(x=P4[k],rhs(eq4)) $k=1..n]: f2 := plot([[P4[k],Q4[k]] $k=1..n],color=blue): display(f1,f2,axes=boxed,title=`fit-approximation in frequency domain`);

>

#HR (four passes) GDD_Hr_x := readdata(`gdd_hr_x.dat`,1,float):  GDD_Hr_y := readdata(`gdd_hr_y.dat`,1,float):   n := vectdim(GDD_Hr_y):    plot([[GDD_Hr_x[k],GDD_Hr_y[k]] $k=1..n],color=red,axes=boxed,title=`measured GDD for HR`);      P5 := array(1..n):       Q5 := array(1..n):        for i from 1 to n do       P5[i] := evalf((2*Pi*3*1e10/(GDD_Hr_x[i]*1e-7)-X_max)/bandwidth):#transition from wavelength to frequency      od:     f1 := plot([[P5[k],GDD_Hr_y[k]*4*(bandwidth*10^(-15))**2] $k=1..n],color=red):    eq5 := fit[leastsquare[[x,y], y=a*x^6+b*x^5+c*x^4+d*x^3+e*x^2+f*x+g]]([[P5[k] $k=1..n], GDD_Hr_y*evalf(4*(bandwidth*10^(-15))**2)]);   Q5 := [subs(x=P5[k],rhs(eq5)) $k=1..n]:  f2 := plot([[P5[k],Q5[k]] $k=1..n],color=blue): display(f1,f2,axes=boxed,title=`fit-approximation in frequency domain`);

Now, as a result of the obtained fit-approximations, we have the normalized net-GDD with corresponding FORTRAN-code for simulation:

>	end_res := evalf(rhs(eq1)+rhs(eq2)+rhs(eq3)+rhs(eq4)+rhs(eq5));  plot(%,x=-0.45..0.8,axes=boxed,title=`normalized to tf net-GDD`);   codegen[fortran](end_res);

      t0 = -0.7685781E4*x**4+0.109791E4*x**2-0.2395398E3*x**3-0.6126327E

     #4*x**6+0.1524603E5*x**5-0.5798741E2-0.5689779E2*x

Or in the usual co-ordinates:

>

P6 := array(1..100):  Q6 := array(1..100):   for i from 1 to 100 do    P6[i] := 2*Pi*3e10*1e7/(bandwidth*(-0.4+0.8*i/100)+X_max):# from frequency to wavelength [in nm]     Q6[i] := evalf(subs(x=-0.4+0.8*i/100,end_res/(bandwidth*1e-15)^2)):      od:       plot([[P6[k],Q6[k]] $k=1..100],color=red,axes=BOXED,title=`net-GDD [fs^2] vs. wavelength [nm]`);

The experimentally observed ultrashort-pulse spectrum is:

>

Spec_x := readdata(`LiSGaF_x.dat`,1,float):  Spec_y := readdata(`LiSGaF_y.dat`,1,float):   n := vectdim(Spec_y):    plot([[Spec_x[k],Spec_y[k]] $k=1..n],color=red,axes=boxed,title=`experimental spectrum (a.u. vs. wavelength)`);     P5 := array(1..n):      Q5 := array(1..n):     for i from 1 to n do    P5[i] := evalf((2*Pi*3*1e10/(Spec_x[i]*1e-7)-X_max)/bandwidth):#transition from wavelength to frequency   Q5[i] := Spec_y[i]:  od: plot([[P5[k],Q5[k]] $k=1..n],color=green,axes=boxed,title=`experimental spectrum (a.u. vs/ dimensionless frequency)`);

We can see a pronounced red shift relatively to the gain band center. The narrow line is the so-called side-band or dispersive wave.

The typical results of the numerical simulations based on the described above model are presented in the next figure (power is given in the arbitrary units).

>

Spec1_x := readdata(`num1_x.dat`,1,float):  Spec1_y := readdata(`num1_y.dat`,1,float):   Spec2_x := readdata(`num2_x.dat`,1,float):    Spec2_y := readdata(`num2_y.dat`,1,float):     Spec3_x := readdata(`num3_x.dat`,1,float):      Spec3_y := readdata(`num3_y.dat`,1,float):       Spec4_x := readdata(`num4_x.dat`,1,float):        Spec4_y := readdata(`num4_y.dat`,1,float):         Spec5_x := readdata(`num5_x.dat`,1,float):          Spec5_y := readdata(`num5_y.dat`,1,float):           Spec6_x := readdata(`num6_x.dat`,1,float):            Spec6_y := readdata(`num6_y.dat`,1,float):           n := vectdim(Spec1_y);          p1 := plot([[Spec1_x[k],Spec1_y[k]] $k=1..n/4],color=black):         n := vectdim(Spec2_y);        p2 := plot([[Spec2_x[k],Spec2_y[k]] $k=3*(n-1)/4..n],color=black):       n := vectdim(Spec3_y);      p3 := plot([[Spec3_x[k],2.5*Spec3_y[k]] $k=3*(n-1)/4..n],color=red):     n := vectdim(Spec4_y);    p4 := plot([[Spec4_x[k],2.5*Spec4_y[k]] $k=1..n/4],color=red):   n := vectdim(Spec5_y);  p5 := plot([[Spec5_x[k],25*Spec5_y[k]] $k=3*(n-1)/4..n],color=blue): n := vectdim(Spec6_y);  p6 := plot([[Spec6_x[k],25*Spec6_y[k]] $k=1..n/4],color=blue):   display(p1,p2,p3,p4,p5,p6,view=0..1,axes=BOXED,title=`numerical spectra (a.u.)`);

In the "natural" co-ordinates we have:

>

n := vectdim(Spec1_y);  p1 := plot([[2*Pi*3e10*1e7/(bandwidth*Spec1_x[k]+X_max),Spec1_y[k]] $k=1..n/4],color=black):   n := vectdim(Spec2_y);    p2 := plot([[2*Pi*3e10*1e7/(bandwidth*Spec2_x[k]+X_max),Spec2_y[k]] $k=3*(n-1)/4..n],color=black):     n := vectdim(Spec3_y);      p3 := plot([[2*Pi*3e10*1e7/(bandwidth*Spec3_x[k]+X_max),2.5*Spec3_y[k]] $k=3*(n-1)/4..n],color=red):       n := vectdim(Spec4_y);      p4 := plot([[2*Pi*3e10*1e7/(bandwidth*Spec4_x[k]+X_max),2.5*Spec4_y[k]] $k=1..n/4],color=red):     n := vectdim(Spec5_y);    p5 := plot([[2*Pi*3e10*1e7/(bandwidth*Spec5_x[k]+X_max),25*Spec5_y[k]] $k=3*(n-1)/4..n],color=blue):   n := vectdim(Spec6_y);  p6 := plot([[2*Pi*3e10*1e7/(bandwidth*Spec6_x[k]+X_max),25*Spec6_y[k]] $k=1..n/4],color=blue): display(p1,p2,p3,p4,p5,p6,view=0..1,axes=BOXED,title=`pulse spectrum vs. wavelength [nm]`);

The black curve corresponds to the soliton propagation without lasing factors, i.e. in the absence of gain, loss saturation and spectral filtering. The high-order dispersion slightly transforms the spectrum of the soliton, but there is not the visible frequency shift. Perhaps, the soliton duration (169 fs) is too large for the high-order dispersions' manifestation. But we can not change the pulse duration for fixed GDD. Such possibility is opened by lasing in the presence of the gain and saturable loss. The red and blue curves correspond to the above described laser model with 1.2 W absorbed pump power for 20x30  pump mode and gaussian mode with 25 m diameter, correspondingly. The pulse durations in these cases are 31 and 27 fs, correspondingly. We can see the appearance of spectral spikes (dispersion waves), which locate near from zeros and 2 m - values of GDD. The domain of the large gradient of GDD in the red spectral region is filled by such spikes. But the pulses' spectra have the very small frequency shifts from gain band maximum, although there is the spectral "shoulder" in the red region. The last corresponds to local GDD extremum near from 883 nm (see above). We had performed the various manipulations with the pump and the fast absorber saturation intensity, but as a result, the self-frequency shift can not be obtained only due to high-order dispersion.

So, there is some additional mechanism of the pulse spectrum shift. As the pulse duration is greater than 10 fs, the nonlinear dispersion can not cause such shift. Therefore we have to investigate the stimulated Raman scattering influence on the pulse spectrum.

Let assume, that the frequency shift results from the stimulated Raman scattering in the active crystal. The vibrational amplitude Q  in the dependence on the pump and signal spectral amplitudes A_m  and A_n , respectively,   obeys the following equation (see H.A.Haus, I.Sorokina, E.Sorokin, J. Opt. Soc. Am. B, 15 , 223 (1998) ):

>

restart:  eq1 := diff(Q(t),t$2)+2*diff(Q(t),t)/T+Omega^2*Q(t)=mu*A_m*exp(I*m*omega*t)*conjugate(A_n)*exp(-I*n*omega*t);   Q = subs({_C1=0,_C2=0},rhs(dsolve(eq1,Q(t))));    amp_Q = Sum(Sum(expand(numer(rhs(%))/expand(denom(rhs(%)))/exp(-I*t*omega*(-m+n))),n=-N/2..N/2),m=-N/2..N/2);#here we supposed that the pulse width << T (T is the phonon relaxation time)

Here m  and  n  are the mode numbers,  is the frequency interval between field spectral components,  is the Raman frequency, T  is the relaxation time, N  is the number of the frequency components in the field spectrum (in the case of our numerical simulations N = ),  is the positive real number expressing photon-phonon coupling. Note, that the last expression for the amplitude of the vibrational oscillations (the frequency of these oscillations is   (m-n)* ) can be re-written as

>	amp_Q = Sum(Sum(mu*A_m*conjugate(A_n)/(Omega^2-2*I*omega*(m-n)/Tr-omega^2*(n-m)^2),n = -N/2 .. N/2),m = -N/2 .. N/2); amp_Q = Sum(Sum(mu*A_m*conjugate(A_n)*(Omega^2-omega^2*(n-m)^2+2*I*omega*(m-n)/Tr)/((Omega^2-omega^2*(n-m)^2)^2+4*omega^2*(m-n)^2/(Tr^2)),n = -N/2 .. N/2),m = -N/2 .. N/2);

Now we have to define the values of the used parameters. There are three Raman lines in LiSGaF with following characteristics: =551, 349, 230 , 1/T  = 6.2, 7.6, 4.2

>

Omega1 := evalf(2*Pi*3*10^10*551):#[Hz]  Omega2 := evalf(2*Pi*3*10^10*349):#[Hz]   Omega3 := evalf(2*Pi*3*10^10*230):#[Hz]    bandwidth := .481859640e15:#[Hz] the normalization for frequencies     Omega1 := evalf(Omega1/bandwidth);#normalized Omega1      Omega2 := evalf(Omega2/bandwidth);#normalized Omega2       Omega3 := evalf(Omega3/bandwidth);#normalized Omega3        omega=evalf(2*Pi*bandwidth/N):#[Hz]         print(`normalized w:`);          2*Pi/N;#normalyzed omega          T1 := evalf(1/(2*Pi*3*10^10*6.2)):#[s]           T2 := evalf(1/(2*Pi*3*10^10*7.6)):#[s]          T3 := evalf(1/(2*Pi*3*10^10*4.2)):#[s]         T1 := evalf(T1*bandwidth);#normalized relaxation time        T2 := evalf(T2*bandwidth);#normalized relaxation time       T3 := evalf(T3*bandwidth);#normalized relaxation time      solve(gain_s=6*omega_s*chi/n_s/c,chi):#if permittivity chi in [cm^2/W], gain_s is the Raman signal gain     mu := simplify(3*omega_s*%*2*Omega/n_s/c/T);    mu1 := subs({Omega = Omega1,T = T1,gain_s=evalf(2.5*1.2e-10),beta= .3379192098e-11},mu*0.8/beta);#normalized  mu2 := subs({Omega = Omega2,T = T2,gain_s=evalf(.185*1.2e-10),beta= .3379192098e-11},mu*0.8/beta);#normalized mu3 := subs({Omega = Omega3,T = T3,gain_s=evalf(.15*1.2e-10),beta= .3379192098e-11},mu*0.8/beta);#normalized

So, we can investigate the field evolution on the basis of the following equation:

>

Diff(A_n(z),z) = Sum(mu*conjugate(A_m)*A_n*(I*(Omega^2-omega^2*(n-m)^2)+2*omega*(m-n)/Tr)/((Omega^2-omega^2*(n-m)^2)^2+4*omega^2*(m-n)^2/(Tr^2))*A_m,m = -N/2 .. N/2); Diff(A_n(z),z) = A_n*Sum(mu*abs(A_m)^2*(I*(Omega^2-omega^2*(n-m)^2)+2*omega*(m-n)/Tr)/((Omega^2-omega^2*(n-m)^2)^2+4*omega^2*(m-n)^2/(Tr^2)),m = -N/2 .. N/2);

Here   n  is the "signal" mode, m  is the "pump" mode, degenerate case m = n  corresponds to the pure self-phase modulation, Stokes signal shift m>n  corresponds to an amplification, anti-Stokes shift m < n  corresponds to a loss of signal wave. Because of the Raman line is narrow in the comparison with pulse spectrum (see normalized T ), we can re-write the expressions:

>

Diff(A_n(z),z) = I*A_n*Sum(mu*abs(A_m)^2/(Omega^2+2*I*omega*(m-n)/Tr-omega^2*(n-m)^2),m = -N/2 .. N/2); Diff(A_n(z),z) = I*A_n*Sum(mu*abs(A_m)^2/((Omega-omega*(n-m))*(Omega+omega*(n-m))+2*I*omega*(m-n)/Tr),m = -N/2 .. N/2); Diff(A_n(z),z) = I*A_n*Sum(mu*abs(A_m)^2/(2*Omega*(Omega-omega*(n-m))+2*I*Omega/Tr),m = -N/2 .. N/2);

The last expression results from the assumption T >>1. In this case

>	Diff(A_n(z),z) = I*A_n*mu*abs(A_m)^2*Int((2*Omega*(Omega-x)-2*I*Omega/Tr)/(4*Omega^2*(Omega-x)^2+4*Omega^2/(Tr^2)),x = -infinity .. infinity);

In the last expression x  is the frequency difference, A_m  corresponds to the pump intensity for - (n-m)= . The contribution of the real part in the integral (the self-phase modulation from the both sides of the Raman line) is equal to 0.

>	-(I*Tr/2/Omega)*Int(1/(x^2+1),x=0..infinity);   value(%);

Thus we obtained the simplest expressions for the field evolution:

>	Diff(A_n(z),z) = A_n*gain*Pi*abs(A_m)^2/4; Diff(A_m(z),z) = -A_m*gain*Pi*abs(A_n)^2/4;

The stimulated Raman gain parameters for numerical simulation are:

>	gain_s1 := evalf(0.8*2.5*1.2e-10/.3379192098e-11);  gain_s2 := evalf(0.8*.185*1.2e-10/.3379192098e-11);   gain_s3 := evalf(0.8*.15*1.2e-10/.3379192098e-11);

Additionally we are to take into consideration the spontaneous Raman scattering as a source for the stimulated scattering. The gain coefficients is this case are [ R.L. Sutherland, Handbook of nonlinear optics, New York (1996), p.305 ] :

>

d_sigma := (4*omega_l/3/omega_s)*3*n_s^2*h*omega_s^3*gain_s/(8*Pi^3*c^2*N)/(1-exp(-h*(omega_l-omega_s)/(2*Pi*kb*Tc)));#kb and Tc are the Boltzmann's constant and temperature, correspondingly  d_sigma1 := evalf(sqrt(subs({kb=1.38*1e-23,Tc=300,omega_s=2*Pi*3e10/.9051e-4,omega_l=2*Pi*3e10/.85e-4,n_s=1.4,h=6.62e-34,gain_s=gain_s1,c=3e10,N=1e20},d_sigma*N)));   d_sigma2 := evalf(sqrt(subs({kb=1.38*1e-23,Tc=300,omega_s=2*Pi*3e10/.88e-4,omega_l=2*Pi*3e10/.85e-4,n_s=1.4,h=6.62e-34,gain_s=gain_s2,c=3e10,N=1e20},d_sigma*N)));    d_sigma3 := evalf(sqrt(subs({kb=1.38*1e-23,Tc=300,omega_s=2*Pi*3e10/.87e-4,omega_l=2*Pi*3e10/.85e-4,n_s=1.4,h=6.62e-34,gain_s=gain_s3,c=3e10,N=1e20},d_sigma*N)));

d_sigma  are the increments of the spontaneous Stokes components (i.e. Raman spontaneous seeds) growth .

As a result of the simulations on the basis of this model, we obtained next spectra:

>

with('linalg'):  with(plots):   R_Spec1_x := readdata(`spec1_x.dat`,1,float):    R_Spec1_y := readdata(`spec1_y.dat`,1,float):     R_Spec2_x := readdata(`spec2_x.dat`,1,float):      R_Spec2_y := readdata(`spec2_y.dat`,1,float):       R_Spec3_x := readdata(`spec3_x.dat`,1,float):        R_Spec3_y := readdata(`spec3_y.dat`,1,float):         R_Spec4_x := readdata(`spec4_x.dat`,1,float):          R_Spec4_y := readdata(`spec4_y.dat`,1,float):           R_Spec5_x := readdata(`spec5_x.dat`,1,float):            R_Spec5_y := readdata(`spec5_y.dat`,1,float):             R_Spec6_x := readdata(`spec6_x.dat`,1,float):              R_Spec6_y := readdata(`spec6_y.dat`,1,float):             n := vectdim(R_Spec1_y);            p1 := plot([[R_Spec1_x[k],R_Spec1_y[k]] $k=1..n/16],color=black):           n := vectdim(R_Spec2_y);          p2 := plot([[R_Spec2_x[k],R_Spec2_y[k]] $k=7*(n-1)/8..n],color=black):         n := vectdim(R_Spec3_y);        p3 := plot([[R_Spec3_x[k],R_Spec3_y[k]] $k=1..n/16],color=red):       n := vectdim(R_Spec4_y);      p4 := plot([[R_Spec4_x[k],R_Spec4_y[k]] $k=7*(n-1)/8..n],color=red):     n := vectdim(R_Spec5_y);    p5 := plot([[R_Spec5_x[2*k],R_Spec5_y[2*k]] $k=1..n/32],color=blue):   n := vectdim(R_Spec6_y);  p6 := plot([[R_Spec6_x[2*k],R_Spec6_y[2*k]] $k=7*(n-1)/16..n/2],color=blue): display(p1,p2,p3,p4,p5,p6,axes=boxed,title=`pulse spectrum vs. dimensionless frequency`);

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

Warning, the name changecoords has been redefined

>

X_max := .2245909122e16:  bandwidth := .481859640e15:   n := vectdim(R_Spec1_y):    p1 := plot([[2*Pi*3e10*1e7/(bandwidth*R_Spec1_x[k]+X_max),R_Spec1_y[k]] $k=1..n/16],color=black):     n := vectdim(R_Spec2_y):      p2 := plot([[2*Pi*3e10*1e7/(bandwidth*R_Spec2_x[k]+X_max),R_Spec2_y[k]] $k=7*(n-1)/8..n],color=black):       n := vectdim(R_Spec3_y):        p3 := plot([[2*Pi*3e10*1e7/(bandwidth*R_Spec3_x[k]+X_max),R_Spec3_y[k]] $k=1..n/16],color=red):         n := vectdim(R_Spec4_y):          p4 := plot([[2*Pi*3e10*1e7/(bandwidth*R_Spec4_x[k]+X_max),R_Spec4_y[k]] $k=7*(n-1)/8..n],color=red):           n := vectdim(R_Spec5_y):            p5 := plot([[2*Pi*3e10*1e7/(bandwidth*R_Spec5_x[2*k]+X_max),R_Spec5_y[2*k]] $k=1..n/32],color=blue):             n := vectdim(R_Spec6_y):            p6 := plot([[2*Pi*3e10*1e7/(bandwidth*R_Spec6_x[2*k]+X_max),R_Spec6_y[2*k]] $k=7*(n-1)/16..n/2],color=blue):          display(p1,p2,p3,p4,p5,p6,axes=BOXED,title=`pulse spectrum vs. wavelength [nm]`);

Three different spectra are presented in the figure. Black curve corresponds to the soliton propagation (200 fs pulse duration) in the presence of stimulated Raman scattering. As a result of the propagation, the Stokes component appears at the main Raman frequency. The strong scattering destroys the soliton after 10 000 cavity transitions. In the laser, a balance between all lasing factors stabilizes the ultrashort pulse. But there is a visible red shift of the field spectrum (red and blue curves) due to stimulated Raman scattering. The pulse spectrum can be pushed from the gain band center (red curve, 25 fs pulse duration) as a result of the Raman self-scattering. There is the possibility of the generation of the additional Stokes lines (blue curve, 56 fs). The frequency shift is comparable with the experimental one. Hence, the stimulated Raman scattering is the main source of the Stokes shift of the pulse spectrum in the Cr:LiSGaF Kerr-lens mode-locked laser (see V.L. Kalashnikov, E. Sorokin, I.T. Sorokina, "Mechanisms of spectral shift in ultrashort-pulse laser oscillators", J. Opt. Soc. Am.   18 , 1732 (2001)).

The collaboration between Maple and external numerical simulators (based on the FORTRAN-code in our case) proves to be extremely fruitful: 1) analytical model building (Maple) --> 2) external code generation (Maple) --> 3) calculation of the simulation parameters (Maple) --> 4) external calculations (FORTRAN-code, the external program can be started from the Maple directly through the system  call.) --> 5) data processing (Maple) --> 6) analytical interpretation of the results (Maple).

As an example of the last step, let's consider the problem of the pulse stability in the Kerr-lens mode-locked laser.

Multipulsing and ultrashort pulse stability

We try to shorten the pulse duration and increase its energy. For this aim we tend the net-GDD to zero (see Part VII). As a result, the ultrashort pulse stability can be lost. Let's consider this phenomenon in detail ( V.L. Kalashnikov, E. Sorokin, I.T. Sorokina, "Multipulse Operation and Limits of the Kerr-lens Mode Locking Stability", submitted to IEEE J. Quantum Electron. (2002)). In the framework of the abberationless approximation the stability loss is revealed as the absence of the soliton-like as well as breezer solution (the evolutional equations for the pulse parameters diverge). What is meaning of this divergence?

The answer comes from the numerical simulations based on the above described model. Let's neglect the stimulated Raman scattering and the higher-order dispersions. In this case, the typical results of the simulations in the vicinity of zero GDD are (contour plots of the field intensity logarithm, vertical axes is t , horizontal axes is z ):

Quasi-soliton consideration fails due to the appearance of the regular or irregular multiple pulse generation. The boundary of the soliton-like pulse nonstability obtained from the numerical simulation is shown here (the parameters in question can be found in the work of reference):

>

restart: with(plots): with(stats): points_numer_x := [-160,-75,-40,-25,-18,-16,-20,-40,-45,-75,-150]:# GDD is normalized to tf  points_numer_y := [0.2,0.5,1,2,5,10,20,30,80,200,500]:# sigma is normalized to the self-phase modulation coefficient beta   statplots[scatterplot](points_numer_x, evalf( map(log10,points_numer_y) ),axes=boxed,color=red,symbol=box):     display(%,color=red,TEXT([-120,1],'`stable single pulse`')):   fig1 := display(%,color=blue,TEXT([-40,2.6],'`multipulsing`'),TEXT([-50,-0.6],'`multipulsing`'),title=`logarithm of boundary sigma vs. GDD`): display(fig1);

Warning, the name changecoords has been redefined

There are the upper and lower on  boundaries defining the transition to multipulsing. There is no single pulse in the vicinity of zero GDD. To interpret this picture let's return to the generalized Landau-Ginzburg equation in the week-nonlinear limit, when  <<1 (see Part VII). For the steady-state ultrashort pulse propagation we have (time and intensity  are normalized):

>	master_1 := 0 = alpha - gamma + I*phi + diff(rho(t),`$`(t,2))/rho(t) + I*k_2*diff(rho(t),`$`(t,2))/rho(t) + gamma*sigma*Phi(t) - I*Phi(t);# rho is the filed, Phi is the intensity  f1 := (t)->rho0*sech(t*tau)^(1+I*psi);# quasi-soliton profile   f2 := (t)->rho0^2*sech(t*tau)^2;# pulse intensity

>	simplify( subs({rho(t)=f1(t),Phi(t)=f2(t)},rhs(master_1)) ):  numer(%):   eq1 := collect(%,cosh(t*tau)^2):    eq2 := evalc( coeff(eq1,cosh(t*tau),2) );     eq3 := evalc( coeff(eq1,cosh(t*tau),0) );

>	eq4 := coeff(eq2,I,0);

As it was shown in Parts VI, VII, the pulse exists if    < 0. Simultaneously, it is condition of the cw suppression (the threshold is not exceeded for the noise out of pulse). As a result, the stability against cw oscillation is provided with:

>	eq5 := factor( eq4 - (alpha-gamma) )/tau^2 > 0;

The pulse chirp can be found from the equation eq3 :

>	eq6 := subs({rho0^2=x,tau^2=y},coeff(eq3,I)=0):  eq7 := subs({rho0^2=x,tau^2=y},coeff(eq3,I,0)=0):   simplify (subs( x=solve(eq7,x),eq6 ) ):#we find the intensity    numer(lhs(%))/y = 0:     sol := solve(%,psi);

In the combination with the condition eq5  we have:

>	eq8 := solve( numer( simplify( subs(psi=sol[1],rhs(eq5)) ) ) = 0, sigma ):  eq9 := solve( numer( simplify( subs(psi=sol[2],rhs(eq5)) ) ) = 0, sigma ):

>	plot([log10(subs(gamma=0.01,eq8[1])),log10(subs(gamma=0.01,eq8[2])),log10(subs(gamma=0.01,eq9[1])),log10(subs(gamma=0.01,eq9[2]))],k_2=-160..0,axes=boxed,color=[green,magenta]): display(%,fig1,view=[-160..0,-1..3]);

Hence the lower boundary of the pulse destabilization is good approximated by the condition of the cw excitation (magenta curve) for the large negative GDD. The green curve corresponds to the chirp-free generation in the soliton model. The intersection of this curve with the stability  boundary (red points) gives the system's parameters corresponding to the minimal pulse duration. However, there are the differences between the analytical model and the numerical results: 1) the ultrashort pulse is chirp-free in the wider region than that predicted from the soliton model; 2) the  increase providing the pulse shortening causes the pulse destabilization (upper on the  parameter boundary of the pulse stability, red points); 3) the pulse is unstable in the vicinity of zero GDD.

The numerical simulations demonstrate, that the pulse destabilization on the upper stability boundary occurs for the negative net-gain coefficient <0. This prevents the cw excitation. Let's consider the toy model of the pulse destabilization in the absence of the cw excitation. For the sake of simplicity, we shall consider the totally real Landau-Ginzburg equation (Part VI).

>	master_2 := rho(t)*g + diff(rho(t),`$`(t,2)) + rho(t)^3*Sigma;#g=alpha-gamma, Sigma=gamma*sigma

Let's expand this equation on the small perturbation :

>	expand( subs( rho(t)=rho(t)+mu*zeta(t),master_2 ) ):  limit( diff(%,mu),mu=0 ):# functional Frechet derivative   master_3 := master_2 + %;

and find its steady-state solutions. Thereto we make the following substitution:

>	f3 := (t)->rho0*sech(t*tau);  f4 := (t)->epsilon*diff( sech(t*tau) ,t$2);#rho0

>	simplify( subs({rho(t)=f3(t),zeta(t)=f4(t)},master_3) ):  expand( numer(%)/rho0):   eq10 := collect( numer(%),cosh(t*tau) );

We can see, that this substitution obeys the perturbed steady-state equation (  is the perturbation amplitude).

>	eq11 := expand( coeff(eq10,cosh(t*tau),0)/epsilon/tau^2 );  eq12 := expand( coeff(eq10,cosh(t*tau),2) );   eq13 := simplify( coeff(eq10,cosh(t*tau),4));

>	sol := solve({eq11=0,eq12=0,eq13=0},{rho0,tau,epsilon});

>	sol1_tau := allvalues(subs(sol[1],tau));  sol1_rho := allvalues(subs(sol[1],rho0));   sol1_e := allvalues(subs(sol[1],epsilon)); sol2_tau := allvalues(subs(sol[3],tau));  sol2_rho := allvalues(subs(sol[3],rho0));   sol2_e := allvalues(subs(sol[3],epsilon));

There exist two types of the perturbed solutions. The first one corresponding to unperturbed solution for the arbitrary small  has a form

>	sol1 := subs( {rho0=sol1_rho[1],tau=sol1_tau[1],epsilon=sol1_e[1]},expand( (f3(t) + f4(t))^2 ) ): plot3d(subs(Sigma=1,sol1),t=-10..10,g=-0.05..0,axes=boxed);

The second solution is:

>	sol2 := subs( {rho0=sol2_rho[1],tau=sol2_tau[1],epsilon=sol2_e[1]},expand( (f3(t) + f4(t))^2 ) ): plot3d(subs(Sigma=1,sol2),t=-10..10,g=-0.05..0,axes=boxed);

And the unperturbed solution is:

>	sol3 := subs( {rho0^2=-4*g/Sigma,tau=sqrt(-g),epsilon=0},expand( (f3(t) + f4(t))^2 ) ): plot3d(subs(Sigma=1,sol3),t=-20..20,g=-0.05..0,axes=boxed);

We can see that the perturbations widen the pulse and reduce its intensity. Moreover, the perturbation of the first type splits  the pulse. We suppose, that the excitation of similar perturbations located within the ultrashort pulse dissociates it for the large , when the contribution of the higher-order nonlinear terms is essential.

The important feature of the multiple pulse regimes is the possibility of a strong correlation of the pulse parameters in the multipulse complex. This is evidence of the pulse interaction. An example of such interaction can be illustrated by the following consideration. Let's take the first momentum of the Landau-Ginzburg equation by the multiplication to the conjugated field, adding the conjugated equation and the consequent integration. The evolution of the energy is described by

>	Diff(E(z),z) = 2*g*E(z) + int(conjugate(rho(z,t))*Diff(rho(z,t),`$`(t,2)) + rho(z,t)*Diff(conjugate(rho(z,t)),`$`(t,2)),t=-infinity..infinity) + 2*Sigma*conjugate(rho(z,t))^2*rho(z,t)^2;

The second term by the virtue of

>	with(student): intparts(int(rho(z,t)*diff(rho(z,t),`$`(t,2)),t), rho(z,t));# the first term vanishes

gives:

>	eq14 := Diff(E(z),z) = 2*g*E(z) - 2*int(Diff(conjugate(rho(z,t)),t)*Diff(rho(z,t),t),t=-infinity..infinity) + 2*Sigma*conjugate(rho(z,t))^2*rho(z,t)^2;

Now let's consider the simplest two-pulse complex:

>	fieldr := rho0*(sech((t-delta)*tau)+sech((t+delta)*tau)*cos(phi));# real part fieldim := rho0*sech((t+delta)*tau)*sin(phi);# imaginary part, delta is the distance, phi is the phase difference print(`the spectral term is defined by:`); diff(fieldr,t)^2+diff(fieldim,t)^2;

The spectral loss for this complex is (second term in eq14 ):

>

assume(tau>0): s := 2*rho0^2*tau^2*( int(sech((t-delta)*tau)^2*tanh((t-delta)*tau)^2,t=-infinity..infinity) + int(sech((t+delta)*tau)^2*tanh((t+delta)*tau)^2,t=-infinity..infinity) + 2*int(sech((t-delta)*tau)*tanh((t-delta)*tau)*sech((t+delta)*tau)*tanh((t+delta)*tau)*cos(phi),t=-infinity..infinity)); en := int(fieldr^2 + fieldim^2,t=-infinity..infinity);

>	plot3d(subs(tau=1/10,s/en),delta=0..40,phi=-Pi..Pi,axes=boxed,view=0..0.01,title=`spectral loss vs. phase and distance`);

So, there exists the potential well, which can absorb the pulses. Also, the last term in eq14  contributes to the interpulse attraction due to the loss saturation enhancement produced by the pulse merging.

We can conclude, that the analytical treatment allowed the comprehension of the basic features of the ultrashort pulse dynamics, though they lie out of the quasi-soliton model validity. Such collaboration between numerical and analytical methods realized by means of the Maple faculties has not only technical but also heuristic character.

Part IX. Mode locking due to a "slow" saturable absorber

Analytical theory and linear stability analysis

In two previous parts we considered the ultrashort pulse formation as result of the loss saturation by pulse intensity. This supposes the instant response of the saturable absorber on the signal variation. The nonresonant (phase) nonlinearity obeys this demand even in femtosecond domain. But the resonant nonlinearities are more inertial and the time defining the relaxation of their excitation lies in the wide region from 100 femtosecond to milliseconds and more. Therefore the interaction of the ultrashort pulse with such structures differs essentially from the previously considered.

Let the pulse duration is shorter then the longitudinal relaxation time. Then the loss saturation is caused by pulse energy flux passed through absorber. The similar situation had considered for dynamical gain saturation in part 3 . Here we shall take into consideration simultaneously the dynamical gain and loss saturation. These effects can be taken into consideration by the expansion of the exponential transmission operator    up to second order on pulse energy  (see H.A. Haus, "Theory of mode locking with a slow saturable absorber", IEEE J. Quant. Electr., 11 , 736 (1975 )). Here  is the gain or loss saturation energy. Then the basic differential equation is

>

restart:  with(plots):   master := diff(rho(z,t),z) = (alpha - g - l)*rho(z,t) - chi*alpha*rho(z,t)*int(rho(z,zeta)^2,zeta=0..t) + alpha*rho(z,t)*(chi*int(rho(z,zeta)^2,zeta=0..t))^2 +  g*rho(z,t)*int(rho(z,zeta)^2,zeta=0..t) - g*rho(z,t)*(int(rho(z,zeta)^2,zeta=0..t))^2 + diff(rho(z,t),t$2) + delta*diff(rho(z,t),t);

Warning, the name changecoords has been redefined

Here  is the ratio of the loss saturation energy to gain saturation energy (saturation parameter), l  is the unsaturable loss coefficient, g  and  are the saturable loss and gain coefficients at pulse peak, respectively,  is the pulse delay on the cavity round-trip. The form of this equation supposes the normalization of time on  , pulse energy on loss saturation energy. We shall suppose the soliton-like form of steady-state solution of master :  

>	f1 := (t)->rho0*sech(t*tau);# soliton form  f2 := (zeta)->rho0*sech(zeta*tau);   ss := rhs(master):    subs({rho(z,t)=f1(t),rho(z,zeta)=f2(zeta)},ss):     simplify(%):      expand( numer(%)*2/rho0 ):       eq := collect( collect( combine(%,trig), sinh(2*t*tau)),cosh(2*t*tau) );

Since this equation is valid at any moment, we have the system of the algebraic equations for the coefficients of hyperbolical functions.

>	eq1 := coeff(eq,cosh(2*t*tau));  eq2 := factor(coeff(eq,sinh(2*t*tau)));   eq3 := expand(eq-eq1*cosh(2*t*tau)-eq2*sinh(2*t*tau));

These equations define the inverse pulse duration , pulse intensity , and delay . Let make some manipulations:

>	eq4 := expand( factor(eq1+eq3)/(2*tau^2) );  eq5 := expand((eq1-eq3)/2);

Hence

>	sol1 := solve(eq4=0,tau^2);# solution for tau^2 solve(subs(tau^4=sol1^2,eq5)=0,rho0^4); factor(%);# solution for rho0^4

Note that the last solution needs some consideration. The comparison with second expression gives for intensity:

 

>	sol2 :=  sqrt(2)*(alpha-g-l)/sqrt(g-chi^2*alpha):# solution for intensity

And at last, for delay we have

>	subs( {tau^2=sol1, rho0^2=sol2},expand(eq2/(2*tau)) ):  simplify(%):   numer(%):    expand(%/(-alpha+g+l)):     sol3 := solve(%=0,delta);# solution for delta

The gain coefficient at pulse peak is:

>	assume(tau,positive):  int(f1(t)^2,t=-infinity..0):   subs({tau=sqrt(sol1),rho0^2=sol2},%):# half-pulse energy    numer( simplify( alpha0/(1+%) - alpha ) ) = 0;# equation for saturated gain, alpha0 is nonsaturated gain

>	sol := solve(%,alpha): tau := 'tau':

Let plot the pulse intensity and duration versus nonsaturated gain for different . It should be noted, that  < 1 because for the pulse formation the loss saturation has to leave behind the gain saturation.

>	plot3d(subs({l=0.01,g=0.05},subs(alpha=sol[1],sol2)),alpha0=0.061..1,chi=0..1,axes=boxed,title=`pulse intensity vs nonsaturated gain`);

>	plot3d(subs({l=0.01,g=0.05},subs(alpha=sol[1],log(1/sqrt(sol1)))),alpha0=0.061..1,chi=0..1,axes=boxed,title=`logarithm of pulse width vs nonsaturated gain`);

>	plot3d(subs({l=0.01,g=0.05},subs(alpha=sol[1],sol3)),alpha0=0.061..1,chi=0..1,axes=boxed,title=`logarithm of pulse width vs nonsaturated gain`);

We can see that the pulse width is decreased by gain growth and decrease of . Additionally there is maximum of the dependence of intensity on saturation parameter .

Now we shall analyze the ultrashort pulse stability in framework of linear theory (see part 3 ). In this case the equation for evolution of exponentially growing perturbation ( t ) is

>

(alpha-g-l-lambda)*xi(t) - chi*alpha*xi(t)*int(rho(zeta)^2,zeta = 0 .. t) - chi*alpha*rho(t)*int(rho(zeta)*xi(zeta),zeta = 0 .. t) + alpha*xi(t)*chi^2*int(rho(zeta)^2,zeta = 0 .. t)^2 + g*xi(t)*int(rho(zeta)^2,zeta = 0 .. t) + g*rho(t)*int(rho(zeta)*xi(t),zeta = 0 .. t) - g*xi(t)*int(rho(zeta)^2,zeta = 0 .. t)^2 + diff(xi(t),`$`(t,2)) + delta*diff(xi(t),t);

Here  is the perturbation's growth increment, its positive value corresponds to ultrashort pulse destabilization. An integro-differential character of this equation raises the ODE's order therefore we shall use some assumptions relatively perturbation's envelope.

Let consider a long-wave limit for perturbations. In this case the perturbation's envelope is smooth in compare with pulse envelope. Then we can to exclude the integration over perturbation.

>

f1 := (t)->rho0*sech(t*tau):# soliton form  f2 := (zeta)->rho0*sech(zeta*tau):   eq1 := (alpha-g-l-lambda)*xi(t)-chi*alpha*xi(t)*int(rho(zeta)^2,zeta = 0 .. t)-chi*alpha*xi0*rho(t)*int(rho(zeta),zeta = 0 .. t)+alpha*xi(t)*chi^2*int(rho(zeta)^2,zeta = 0 .. t)^2+g*xi(t)*int(rho(zeta)^2,zeta = 0 .. t)+g*xi0*rho(t)*int(rho(zeta),zeta = 0 .. t)-g*xi(t)*int(rho(zeta)^2,zeta = 0 .. t)^2+diff(xi(t),`$`(t,2))+delta*diff(xi(t),t);# xi0 is the amplitude of perturbation at t=0

The long-wave approximation allows to neglect the second-order derivation in compare with first-order one.

>	value( subs({rho(t)=f1(t), rho(zeta)=f2(zeta)},eq1-diff(xi(t),`$`(t,2))) );

>	subs({op(2,%)=-chi*alpha*xi(t)*rho0^2*tanh(t*tau), op(4,%)=chi^2*alpha*xi(t)*rho0^4*tanh(t*tau)^2, op(5,%)=g*xi(t)*rho0^2*tanh(t*tau), op(7,%)=-g*xi(t)*rho0^4*tanh(t*tau)^2},%);

>	dsolve({%=0, xi(0)=xi0} ,xi(t));

The last term in this expression has not appropriate asymptotic behavior at infinity. As consequence, there are not long-wave excitations in our case.

A short-wave approximation allows to simplify the problem on the basis of Riemann-Lebesgue theorem:  =o(1) ( --> ). That is the integral  is small value if  is quickly oscillating function without steady-state points. In this case (when <<1 ) we have (see V.L.Kalashnikov, V.P.Kalosha, I.G.Poloyko, V.P.Mikhailov, "Mode locking of cw solid-state lasers with a slow saturable absorber", Opt. and Spectroscopy, 81  794 (1996) ):

>

f1 := (t)->rho0*sech(t*tau):# soliton form  f2 := (zeta)->rho0*sech(zeta*tau):   eq1 := (alpha-g-l-lambda)*xi(t) - chi*alpha*xi(t)*int(rho(zeta)^2,zeta = 0 .. t) + alpha*xi(t)*chi^2*int(rho(zeta)^2,zeta = 0 .. t)^2 + g*xi(t)*int(rho(zeta)^2,zeta = 0 .. t) - g*xi(t)*int(rho(zeta)^2,zeta = 0 .. t)^2 + diff(xi(t),`$`(t,2)) + delta*diff(xi(t),t);

>	value( subs({rho(t)=f1(t), rho(zeta)=f2(zeta)},eq1) );

We can rewrite this expression:

>

(alpha-g-l-lambda)*xi(t)-chi*alpha*xi(t)*rho0^2*tanh(t*tau)/tau+alpha*xi(t)*chi^2*rho0^4*tanh(t*tau)^2/tau^2+g*xi(t)*rho0^2*tanh(t*tau)/tau-g*xi(t)*rho0^4*tanh(t*tau)^2/tau^2+diff(xi(t),`$`(t,2))+delta*diff(xi(t),t):  eq2 := collect(%, tanh);   A1 = expand( coeff(eq2, tanh(t*tau)^2)/ xi(t)):    A2 = expand( coeff(eq2, tanh(t*tau))/ xi(t)):     eq3 := subs( alpha-g-l-lambda=A3,A1*xi(t)*tanh(t*tau)^2 + A2*xi(t)*tanh(tau*t) + coeff(eq2,tanh(t*tau),0) );

>	dsolve(eq3=0, xi(t));

>

coeff1 := allvalues( -RootOf(A1-_Z*tau^2+_Z^2*tau^2)+1+RootOf(A3-A2+A1+2*delta*_Z*tau+4*_Z^2*tau^2)+RootOf(4*_Z^2*tau^2-2*delta*_Z*tau+A2+A3+A1) ):# 8 coefficients coeff2 := allvalues( -1/2*(2*RootOf(A1-_Z*tau^2+_Z^2*tau^2)*tau-2*tau+delta-2*RootOf(4*_Z^2*tau^2-2*delta*_Z*tau+A2+A3+A1)*tau+2*RootOf(A3-A2+A1+2*delta*_Z*tau+4*_Z^2*tau^2)*tau)/tau ):# 8 coefficients

We have the set of 16 first coefficients of hypergeometric functions, which can cause the appropriate asymptotic behavior at infinity if they are equal to negative integers (see, in particular, part 3 ). Moreover, they are to be the large negative integers in order to satisfy to our short-wave approximation. So, we investigate the high-level excitations in the "potential well" formed by gain and loss saturation.

>

#first coefficients in hypergeometric functions  c := array(1..16):   for i from 1 to 8 do    c[i] := subs({A1 = -g*rho0^4/(tau^2)+alpha*chi^2*rho0^4/(tau^2),A2 = g*rho0^2/tau-chi*alpha*rho0^2/tau,A3 = alpha-g-l-lambda},coeff1[i]):     od:      for i from 9 to 16 do       c[i] := subs({A1 = -g*rho0^4/(tau^2)+alpha*chi^2*rho0^4/(tau^2),A2 = g*rho0^2/tau-chi*alpha*rho0^2/tau,A3 = alpha-g-l-lambda},coeff2[i-8]):         od:

>	#16 coefficients produce 16 equations for lambda (N is positive integer)  s := array(1..16):   for j from 1 to 16 do    s[j] := solve(c[j] + N = 0, lambda):     od:

>

# the solutions will be evaluated numerically by variation of chi for different N # Attention! This computational block can take a lot time! P := array(1..50,1..16,1..3):  l := 0.01:   g := 0.05:    alpha0 := 0.5:     Lev := [5, 10, 50]: for k from 1 to 3 do for j from 1 to 16 do for i from 1 to 50 do  N := Lev[k]:   chi := i/50:    alpha := evalf( sol[1] ):     tau := evalf( sqrt(sol1) ):      rho0 := evalf( sqrt(sol2) ):              delta := evalf( sol3 ):        P[i,j,k] := evalf( s[j] ):         od: od: print(k); od:

>	# list of plots macro(usercol = COLOR(RGB, 0.8, 5/Lev[k], 0.5)):  for k from 1 to 3 do   for m from 1 to 16 do    p[(k-1)*16+m] := listplot([[n/50,Re(P[n,m,k])] $n=1..50],color=usercol):     od:      od:

>	display({p[ii] $ii=1..40},axes=boxed,title=`stability increment`,view=-2..15);

So, we have only two different solutions for real part of , when "level's number" N  is fixed. The negative value of increment corresponds to decaying perturbations, i. e. stable ultrashort pulse generation. Because of the pulse stabilization results from increase of   there is problem of the shortest pulse generation (note, that increase of  increases pulse width, see corresponding figure above). It is of interest that the region of the pulse stability can have an inhomogeneous character (a "kink" on the curves for small N , that is "deeper level" in potential well). This fact was confirmed by numerical simulations (see V.L.Kalashnikov, D.O.Krimer, I.G.Poloyko, "Soliton generation and picosecond collapse in solid-state lasers with semiconductor saturable absorber", J. Opt. Soc. Am. B. 17 , 519 ( 2000 )) and corresponds to destabilization due to excitation of "deeper levels", that is the pulse envelope splitting. Perturbations with large N  can be interpreted as continuous-wave generation of noise spikes.

Aberrationless approximation

Now we shall consider one practical realization of the described here method of ultrashort pulse generation. The most attractive sort of "slow" absorber is the semiconductor absorber with comparatively short relaxation time and small energy of saturation (see next part). But the absorber, which is based on the impure dielectric crystal, is much more usable in the picosecond time domain due to its simplicity (as rule, the semiconductor shutter has a very complicate structure) and durability. However, the basic disadvantages of the crystalline absorbers are the large relaxation time (10 ns - 1 s) and large saturation energy (up to 1 J / ). Here we shall analyze the possibility of the stable ultrashort pulse generation in the forsterite solid-state laser with YAG: V(3+) crystalline absorber. As the basic method the aberrationless approximation will be used. Now we shall express the loss coefficient in absorber as , that differs from phenomenological expression, which was considered in the beginning of this part, but is close to it when the expansion on energy is made up to small orders of . Such form of the dependence can be obtained from Bloch's equations for two-level absorber (see next part) in noncoherent approximation and in the condition of small pulse duration in compare with longitudinal relaxation time Ta   ( Ta = 22 ns for YAG: V(3+)). The steady-state (there is no dependence on z ) master equation is (see above):

>	restart: master:=alpha*(rho(t)-rho(t)*chi*epsilon(t)+rho(t)*chi^2*epsilon(t)^2/2)-g*(rho(t)-epsilon(t)*rho(t)+rho(t)*epsilon(t)^2/2)+diff(diff(rho(t),t),t)+delta*diff(rho(t),t)-l*rho(t);

Let suppose the soliton-like shape of the pulse. Then the evolution of the ultrashort pulse parameters can be found as result of the expansion at t -series:

>

f1:=(z,t)->rho0(z)*sech(t*tau(z));# pulse amplitude  f2:=(z,t)->rho0(z)^2*(1+tanh(t*tau(z)))/tau(z);# pulse energy   lhs_master:=subs({diff(rho0(z),z)=x,diff(tau(z),z)=y},diff(f1(z,t),z)):# left-hand side of dynamical equation   eq:=collect(series(lhs_master-subs({rho(t)=f1(z,t),epsilon(t)=f2(z,t)},master),t=0,3),t):#dynamical equation eq1 := coeff(eq,t^2): eq2 := coeff(eq,t): eq3 := coeff(eq,t,0): sol := factor(solve({eq1, eq2, eq3},{x, y, delta}));

The equations for the pulse parameters evolution have to be supplemented by the equations for the gain ( --> ( z ) ) and the saturable loss ( g --> g ( z ) ) evolution at the time intervals >>1/  ( 0 and gmx  are the maximal saturable gain and loss, respectively, Ta  and Tr  are the loss and gain relaxation times normalized to the cavity period Tcav , Pump  is the dimensionless pump, see part 8 ):

>	eq4:=-4*rho0(z)^2*g(z)/tau(z)+(gmx-g(z))/Ta; eq5:=Pump*(alpha0-alpha(z))-2*chi*alpha(z)*rho0(z)^2/tau(z)-alpha(z)/Tr;

Then finally:

>	sys := D(g)(z)=eq4,D(a)(z)=eq5,D(rho0)(z)=subs({alpha=alpha(z),g=g(z)},subs(sol,x)), D(tau)(z)=subs({alpha=alpha(z),g=g(z)},subs(sol,y));# basic systems for evolution of the pulse parameters, gain and loss coefficients

We shall change the saturation parameter   and to search the stationary points of the pulse parameter's mapping. These points correspond to the solutions of 'sys'  with zero left-hand sides

>

st_sol1:=solve({eq4,eq5},{g(z),alpha(z)}):  st_g:=subs(st_sol1,g(z));   st_a:=subs(st_sol1,alpha(z));    st_sys1:=[expand(rhs(op(3,[sys]))*2*tau(z)^2/rho0(z)),\     expand(rhs(op(4,[sys]))*2*tau(z))]:      st_sys2:={simplify(op(1,st_sys1)+op(2,st_sys1))=0,op(2,st_sys1)=0};       st_sys3:=subs({rho0(z)=x,tau(z)=y},subs({g(z)=st_g,alpha(z)=st_a},st_sys2));

The next procedure will be used for numerical solution of st_sys3 .

>

num_sol := proc(alpha0,gmx,chi,l,Tr,Ta,Pump)     st_sys := {-Pump*y*Tr*alpha0*x^4*chi^2/(Pump*y*Tr+2*chi*x^2*Tr+y)+y*gmx*x^4/(4*x^2*Ta+y)-4*y^4 = 0, 2*Pump*y^3*Tr*alpha0/(Pump*y*Tr+2*chi*x^2*Tr+y)-2*Pump*y^2*Tr*alpha0*x^2*chi/(Pump*y*Tr+2*chi*x^2*Tr+y)-2*l*y^2-2*y^3*gmx/(4*x^2*Ta+y)+2*y^2*gmx*x^2/(4*x^2*Ta+y)-6*y^4 = 0}:     fsolve(st_sys,{x,y},{x=0..1,y=0..1}): end:

Warning, `st_sys` is implicitly declared local to procedure `num_sol`

>	v := array(1..100):  for i from 1 to 100 do   v[i] := num_sol(0.5,0.05,1/(1.36+9*i/100),0.01,300,2.2,0.001) od:# the normalization of relaxation times to cavity period is supposed (Tcav = 10 ns)

Now we can plot the logarithm of the pulse duration versus .

>

with(plots):  ww:=array(1..100):   for j from 5 to 100 do ww[j]:=evalf(log10(1/subs(v[j],y))) od:    points := {seq([1/(1.36+9*j/100),ww[j]],j=1..100)}:     f1 := plot(points,x=0.1..0.6,style=point,symbol=circle,color=red):      display(f1,TEXT([10,3],'`Up=0.0008`'),view=2.4..3.6,title=`Logarithm of pulse duration versus sigma`,labels=[`saturation parameter`,``]);

Warning, the name changecoords has been redefined

This figure demonstrates the pulse width decrease due to decrease of ,  that corresponds to the predominance of the loss saturation over gain saturation (see previous subsection). For the time normalization to the inverse bandwidth of YAG: V(3+) absorption line the pulse duration at  = 0.3  is about of 50 ps.

Now, we shall consider the ultrashort pulse parameters evolution on the basis of the obtained system of ODE. This procedure solves the systems by the standard operator DEplot

>

with(DEtools):  ODE_plot := proc(alpha0,gmx,chi,l,Tr,Ta,Pump)   sys := [D(g)(z) = -4*rho0(z)^2*g(z)/tau(z)+(gmx-g(z))/Ta, \    D(alpha)(z) = Pump*(alpha0-alpha(z))-2*chi*alpha(z)*rho0(z)^2/tau(z)-alpha(z)/Tr, \     D(rho0)(z) = 1/2*rho0(z)*(-2*l*tau(z)^2+2*alpha(z)*tau(z)^2-2*alpha(z)*rho0(z)^2*chi*tau(z)+alpha(z)*rho0(z)^4*chi^2-2*tau(z)^4-2*g(z)*tau(z)^2+2*g(z)*rho0(z)^2*tau(z)-g(z)*rho0(z)^4)/(tau(z)^2),\       D(tau)(z) = -1/2*(alpha(z)*rho0(z)^4*chi^2-g(z)*rho0(z)^4+4*tau(z)^4)/tau(z)]: DEplot(sys,[rho0(z),tau(z),g(z),alpha(z)],z=0..10000,[[rho0(0)=1e-5,tau(0)=1e-3,g(0)=gmx,alpha(0)=0]],stepsize=1,scene=[z,rho0(z)],axes=FRAME,linecolor=BLACK): end:

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

Let vary the saturation parameter  for the fixed pump:

>	display(ODE_plot(0.5,0.05,0.3,0.01,300,2.2,0.001));

So, we have the decaying oscillations of the ultrashort pulse amplitude (see part 6 ). The decrease of saturation parameter produces

>	display(ODE_plot(0.5,0.05,0.1,0.01,300,2.2,0.001));

One can see the growth of the ultrashort pulse oscillations (autooscillations or so-called Q-switch mode locking). This is a main destabilizing factor for the considered type of the mode locking regime. This result is in agreement with above obtained (note, that now we are outside of framework of linear perturbation theory).

The important obstacle for the pulse generation by slow absorber is the noise growth. In order to investigate pulse stability against noise in framework of the considered model we have to add the equation for the evolution of noise energy n  to sys.

>

ODE_noise := proc(alpha0,gmx,chi,l,Tr,Ta,Pump) sys_noise := [D(g)(z) = -4*rho0(z)^2*g(z)/tau(z)+(gmx-g(z))/Ta, \    D(alpha)(z) = Pump*(alpha0-alpha(z))-2*chi*alpha(z)*rho0(z)^2/tau(z)-alpha(z)/Tr, \     D(rho0)(z) = 1/2*rho0(z)*(-2*l*tau(z)^2+2*alpha(z)*tau(z)^2-2*alpha(z)*rho0(z)^2*chi*tau(z)+alpha(z)*rho0(z)^4*chi^2-2*tau(z)^4-2*g(z)*tau(z)^2+2*g(z)*rho0(z)^2*tau(z)-g(z)*rho0(z)^4)/(tau(z)^2),\       D(tau)(z) = -1/2*(alpha(z)*rho0(z)^4*chi^2-g(z)*rho0(z)^4+4*tau(z)^4)/tau(z),\         D(n)(z) = (alpha(z)-l-(gmx+Ta*(g(z)-gmx)*(1-exp(-1/Ta))))*n(z)]:#see V.L. Kalashnikov et al. Opt. Commun., v.159, 237 (1999)   DEplot3d(sys_noise,[rho0(z),tau(z),g(z),alpha(z),n(z)],z=0..10000,[[rho0(0)=1e-5,tau(0)=1e-3,g(0)=gmx,alpha(0)=0,n(0)=1e-5]],stepsize=1,scene=[z,rho0(z),n(z)],axes=FRAME,linecolor=BLACK): end:

Warning, `sys_noise` is implicitly declared local to procedure `ODE_noise`

>	display(ODE_noise(0.5,0.05,0.3,0.01,300,2.2,0.001),title=`Noise energy evolution`);

As one can see from this picture, the noise energy decays, i. e. 50 ps "auto-stable" pulse is stable against noise, too.