Propagation of the Gaussian ultrashort pulse

selff.mws

Propagation of the Gaussian ultrashort pulse in a nonlinear laser medium

V.L. Kalashnikov

Photonics Institute, TU Vienna, Gusshausstr. 27/387, A-1040, Vienna, Austria

Abstract . Propagation of the Gaussian ultrashort pulse in a nonlinear medium is considered on the basis of the variation and momentum methods. The combined action of self-phase modulation, self-focusing, dispersion, diffraction, spectral filtering, radial varying gain and thermo-lensing is considered. The presented approaches are of interest for an analysis of the Kerr-lens mode-locking ability and the ultrashort pulse stability of the continuous-wave pumped solid-state lasers.

Nonlinear propagation of the ultrashort pulse is governed by the different factors affecting both temporal and spatial field characteristics. Usually, we can segregate the spatial factors from the temporal ones. Then the dynamical equation does not contain the mixed derivatives and the propagating pulse can be approximated by the Gaussian spatial-time profile, which form remains invariable ( aberrationless approximation ). The Gaussian spatial profile is most interesting because of it describes the fundamental mode of a laser resonator. Such beam can be considered both within and out of the resonator. The former case involves the dissipative processes in the active laser medium whereas the latter can be considered as the non-dissipative process (for example, intra-fiber propagation). For the time profile of the propagating pulse the Gaussian approximation is not necessary and can be replaced by the sech -approximation, which describes the laser quasi-soliton. This replacement only slightly modifies the below obtained results and will be not considered here. As a rule, the analysis is performed in the framework of the paraxial (parabolic) approximation (see, for example, Y.R. Shen, "The principles of nonlinear optics", John Wiley and Sons, New-York, 1984 and J.Jasapara, V.L.Kalashnikov, D.O.Krimer, I.G.Poloyko, M.Lenzner, W.Rudolph, J. Opt. Soc. Am. B 17 , 319 (2000) ) that imposes the limitations on the validity of the obtained analytical results. Here we shall consider the approaches, which are based only on the aberrationless approximation beyond the bounds of the parabolic approximation and thereby agree with the numerical and experimental data within the wider parametric range.

References:

1. D. Anderson, "Variational approach to nonlinear pulse propagation in optical fibers", Phys. Rev. A 27, 3135 (1983).

2. M. Karlsson, D. Anderson, M. Desaix, M. Lisak, "Dynamic effects of Kerr nonlinearity and spatial diffraction on self-phase modulation of optical pulses", Optics Lett. 16, 1373 (1991).

3. M. Desaix, D. Anderson, M. Lisak, "Variational approach to collapse of optical pulses", J. Opt. Soc. Am. B 8, 2082 (1991).

4. S.N. Vlasov, V.A. Petrishchev, V.I. Talanov, "Averaged description of wave beams in linear and nonlinear media", Izvestija vuzov (radiophizika) 14, 1353 (1971), in russian.

5. J. Herrmann, "Theory of Kerr-lens mode locking: role of self-focusing and radially varying gain", J. Opt. Soc. Am. B 11, 498 (1994).

As we shall consider the variation approach therefore it is convenient to use the QFT-package , which can be downloaded from the University of Waterloo site.

> restart:

interface( showassumed=0 ):

libname := libname, `c:\\Program Files\\Maple 8\\QFT`:# ! QFT package has to be installed !

with(QFT):
with(student):

dimension(3, `+`):# we shall consider 1+2 dimensional dynamical equations

$`Default dimension of the Tensor-Space is set to: 4 ,\` - \` `$

${Dc, d_, gamma, sigma, epsilon, delta, g_}, {}$

$`for \\partial_{mu} (differentiation) use d_[mu]; it will be displayed as`, `ð`[mu]$

Warning, the protected name Dirac has been redefined and unprotected

$`Default dimension of the Tensor-Space is set to: 3 ,\` + \` `$

${Dc, d_, gamma, sigma, epsilon, delta, g_}, {}$

As the Gaussian beams have the cylindrical symmetry, we consider only longitudinal, radial and time coordinates, i.e. , and , respectively.

> coordinates();

The considered nonlinearity is of the Kerr type and transforms the field phase. Hence the propagation equation is complex. The slowly varying field amplitude is and its complex conjugation is .

>    with(PDEtools):
declare(quiet):
  ON;
   interface(imaginaryunit=i);# imaginary unit
    macro(df=diff):

Define(A);
declare(A(X)):# X[1] is the propagation distance , X[2] is the radial coordinate (X[2]=0 corresponds to the beam axis), X[3] is the local time (X[3]=0 corresponds to the pulse maximum)

`print/A1` := () -> if args=(X) then (A^`*`) else (A^`*`)(args) fi;

${Dc, d_, gamma, sigma, A, epsilon, delta, g_}, {}$

> A=A(X),A1=A1(X);# field envelope and its complex conjugation

> # Henceforth the derivatives look as follows:

diff(A(X),X[1]);
diff(A1(X),X[1]);

Variational method

At first let us consider the non-dissipative propagation of the Gaussian pulse on the basis of the variational approach . The master equation is the generalized nonlinear Schroedinger equation in the cylindrical coordinates:

> eq1 := df(A(X),X[1]) - 2*i*K*X[2]^2/wp^2*A(X) + i*(1/X[2]*df( (X[2]*df(A(X),X[2])),X[2] ) )/2/k + i*D/2*df(A(X),X[3]$2) = - i*kappa*abs(A(X))^2*A(X);# generalyzed nonlinear Schroedinger equation

eq2 := df(A1(X),X[1]) + 2*i*K*X[2]^2/wp^2*A1(X) - i*(1/X[2]*df( (X[2]*df(A1(X),X[2])),X[2] ) )/2/k - i*D/2*df(A1(X),X[3]$2) = i*kappa*abs(A(X))^2*A1(X);# complex conjugate equation

Here the spatial derivatives form the transverse Laplacian in the cylindrical coordinates describing the beam diffraction ( k is the wave number), the time-derivatives describe the group-velocity dispersion ( D is the group-velocity dispersion coefficient and its positive value corresponds to the anomalous dispersion). The right-hand sides describe the effect of the field self-focusing and self-phase modulation due to the Kerr nonlinearity, is the nonlinear coefficient ( and are the nonlinear and linear refraction indexes, respectively). The thermo-lensing induced by the pump beam is considered in the framework of the aberrationless approximation, ( a is the loss coefficient, P is the pump power, is the thermal conductivity, see J.Frauchiger, P.Albers, H.P.Weber, IEEE J. Quantum Electron. 28, 1046 (1992)). These non-dissipative equations can be obtained from the Lagrangian, whose form is suggested by the following transformation:

> expand( value( eq1*X[2]*A1(X) - eq2*X[2]*A(X) ) ):# transformation

eq3 := Int( lhs(%) - rhs(%),X[2] ) = 0;

The integration by parts transforms the radial derivatives:

>    intparts(Int(1/2*i*X[2]*A1(X)/k*df(A(X),`$`(X[2],2)), X[2]), X[2]*A1(X)):
  expand(%);

intparts(Int(1/2*i*X[2]*A(X)/k*df(A1(X),`$`(X[2],2)), X[2]), X[2]*A(X)):
  expand(%);

This cancels the first-order derivatives and transforms the second-order derivative into . The similar transformation is valid also for the dispersion term. Hence the Langrangian density is:

> # Lagrangian density

L[1] := X[2]*(A1(X)*df(A(X),X[1]) - A(X)*df(A1(X),X[1])) - i*X[2]*df(A(X),X[2])*df(A1(X),X[2])/k + i*X[2]*kappa*A(X)^2*A1(X)^2 - 4*i*X[2]*K*X[2]^2/wp^2*A(X)*A1(X) - i*X[2]*D*df(A(X),X[3])*df(A1(X),X[3]);

The Schroedinger equations ( eq1 and eq2 ) result from the Euler-Lagrange equations:

> # Euler-Lagrange equations

Leq := (L,q) -> pdiff(L,q) - pdiff(L, df(q,X[1]), X[1]) - pdiff(L, df(q,X[2]), X[2]) - pdiff(L, df(q,X[3]), X[3]);# general form of the Euler-Lagrange equations

eq4 := collect( expand( Leq(L[1],A(X))/X[2]/2 ),K ) = 0;# generalyzed nonlinear Schroedinger equation
eq5 := collect( expand( Leq(L[1],A1(X))/X[2]/2 ),K ) = 0;# generalyzed nonlinear Schroedinger equation

Now let us consider the trial solution in the form of the propagating Gaussian pulse. This allows integrating out the radial and time dependences from the Lagrangian.

>    assume(b(X[1]),real):
assume(a(X[1]),real):

eq6 := A(X) = W(X[1])*exp( i*b(X[1])*X[2]^2 - X[2]^2/a(X[1])^2/2 + i*psi(X[1])*X[3]^2 - X[3]^2/tau(X[1])^2 );# Gaussian pulse with complex amplitude W, inverse square of the wavefront curvature radius b, beam size a, pulse width tau and chirp psi
eq7 := A1(X) = `W*`(X[1])*exp( - i*b(X[1])*X[2]^2 - X[2]^2/a(X[1])^2/2 - i*psi(X[1])*X[3]^2 - X[3]^2/tau(X[1])^2);# complex conjugate solution

Int( Int( L[1], X[2]=0..infinity ), X[3]=-infinity..infinity );# integration of the Lagrangian density over X[2] and X[3]
subs({A(X)=rhs(eq6),A1(X)=rhs(eq7)},%):
  simplify(value(%)):
   f := expand( i*numer(%)/2/sqrt(2)/(csgn(conjugate(tau(x1)))*sqrt(Pi))/tau(X[1])^2/wp^2 );

The last expression is the Lagrangian expressed through the trial Gaussian solution. This Lagrangian can be expressed in the following form:

> L[2] := i*tau(X[1])*a(X[1])^2*( `W*`(X[1])*df(W(X[1]),X[1]) - W(X[1])*df(`W*`(X[1]),X[1]) ) - (W(X[1])*`W*`(X[1]))*( 2*tau(X[1])*a(X[1])^4*df(b(X[1]),X[1]) - tau(X[1])/k - 4*a(X[1])^4*b(X[1])^2*tau(X[1])/k + sqrt(2)*kappa*tau(X[1])*a(X[1])^2*(W(X[1])*`W*`(X[1]))/4 - 4*K*tau(X[1])*a(X[1])^4/wp^2 + tau(X[1])^3*a(X[1])^2*df(psi(X[1]),X[1])/2 - D*tau(X[1])^3*a(X[1])^2*psi(X[1])^2 - D*a(X[1])^2/tau(X[1]));

simplify(f*tau(X[1])/k - L[2]);# check up

The Euler-Lagrange equations for the evaluating pulse parameters are:

>    Leq := (L,q) -> pdiff(L,q) - pdiff(L, df(q,X[1]), X[1]):# Euler-Lagrange equation for the pulse parameters

eq8 := Leq(L[2],W(X[1])) = 0;
eq9 := Leq(L[2],`W*`(X[1])) = 0;
  eq10 := Leq(L[2],a(X[1])) = 0;
  eq11 := Leq(L[2],b(X[1])) = 0;
eq12 := Leq(L[2],psi(X[1])) = 0;
eq13 := Leq(L[2],tau(X[1])) = 0;

Thus, we reduced the initial generalized nonlinear Schroedinger equation to the system of the first-order ordinary differential equations describing the evolution of the pulse parameters.

eq14 gives the energy conservation law :

> expand( -(eq8*W(X[1]) - eq9*`W*`(X[1]))/2/i );

eq14 := Diff( W(X[1])*`W*`(X[1])*tau(X[1])*a(X[1])^2,X[1] ) = 0;

There is the equation for the evolution of the field amplitude :

> simplify( eq8*W(X[1]) + eq9*`W*`(X[1]) ):
expand( numer(lhs(%))/2/wp^2/k ) = 0;

eq15 := i*a(X[1])^2*tau(X[1])^2*( `W*`(X[1])*df(W(X[1]),X[1]) - W(X[1])*df(`W*`(X[1]),X[1]) ) = -simplify(lhs(%)-i*a(X[1])^2*tau(X[1])^2*( `W*`(X[1])*df(W(X[1]),X[1]) - W(X[1])*df(`W*`(X[1]),X[1]) ));

Taking into account this equation and eq10 we obtain the equation for the wave-front curvature :

>    eq10/2;
rhs(eq15) + tau(X[1])*a(X[1])*op(2,lhs(%)):
  simplify(%):
   numer(%)/(W(X[1])*`W*`(X[1])*tau(X[1])^2):
    fun[1] := df(b(X[1]),X[1]) =expand( solve(%,df(b(X[1]),X[1])) );

The energy conservation law and eq11 allow the equation for the Gaussian beam size :

>    eq11/2;
lhs(%) - lhs( expand( a(X[1])^2*value(eq14) ) ):
  simplify(%):
   fun[2] := df(a(X[1]),X[1]) = solve(%,df(a(X[1]),X[1]));

And, at last, the similar manipulations result in the equations for the chirp and the pulse duration :

> rhs(eq15) + tau(X[1])^2*op(2,lhs(eq13)):
simplify(%);
fun[3] := df(psi(X[1]),X[1]) = expand( solve(%,df(psi(X[1]),X[1])) );

>    2*eq12;
lhs( % );
  lhs( expand( tau(X[1])^2*value(eq14) ) );
   simplify(%% - %);
    fun[4] := df(tau(X[1]),X[1]) = solve(%,df(tau(X[1]),X[1]));

The resulting system is:

>    eq[1] := abs(W(X[1]))^2 = abs(W0)^2*exp(2*G(X[1]))*tau0*a0^2/(tau(X[1])*a(X[1])^2);# from eq14
eq[2] := fun[4];
  eq[3] := fun[2];
   eq[4] := lhs(fun[3]) = 2*D*(psi(X[1])^2 - 1/tau(X[1])^4) + sqrt(2)/4*kappa*rhs(eq[1])/tau(X[1])^2;
    eq[5] := lhs(fun[1]) = 2*K/wp^2 - 1/2*1/(a(X[1])^4*k) + 2/k*b(X[1])^2 + sqrt(2)/8*kappa*rhs(eq[1])/a(X[1])^2;

subs({W(X[1])=sqrt( rhs(eq[1]) )*exp(i*phi(X[1])),`W*`(X[1])=sqrt( rhs(eq[1]) )*exp(-i*phi(X[1]))},eq15):# phi is the radially independent phase
value(%):
  simplify(%):
   solve(%,df(phi(X[1]),X[1])):# equation for df(phi(X[1]),X[1])
    subs({eq[4],eq[5]},%):
      eq[6] := df(phi(X[1]),X[1]) = expand(%);

Here we supposed that .

Momentums method

The considered variation approach is very powerful for the analysis of the non-dissipative systems allowing the Lagrangian description. However, the dissipative factors inherent in the laser systems require the analysis on the basis of the so-called momentums method . It should be noted, that this approach needs some caution and can result in the nonlinear terms, which numerical coefficients slightly differ from those obtained on the basis of the variation approach. Hence we have to control the validity of the obtained results by the numerical simulation of the initial dynamical equation (such MatLab-based control simulation will be presented elsewhere). Also, we don't agree with the approach of Ref. 5, where the analysis is based on the consideration of the arbitrary field momentums . As a result, the spatial parameters are affected directly by the dispersion ( V.P.Kalosha, M.Mueller, J.Herrmann, S.Gatz, J. Opt. Soc. Am. B 15, 535 (1998)). This is not possible in the framework of the considered model (see above; there exists certainly the spatial-time connection through the nolinearity).

Our approach is based on the analysis of the momentums implied by the symmetries of the master dynamical equation. The --> -invariance, the transverse and time translating invariances suggest the analysis of the following momentums (see also N.N.Akhmediev, A.Ankiewicz, " Solitons: Nonlinear Pulses and Beams " (Chapman&Hall, London, 1997)): , and .

Now let us consider the zero-order momentum . The substitution of the trial functions eq6, eq7 for the pulse profile gives:

>    # T[0,0]

subs( {eq6,eq7},lhs(eq1) - rhs(eq1) )*rhs(eq7) + subs( {eq6,eq7},lhs(eq2) - rhs(eq2) )*rhs(eq6):
simplify( value(%) ):
  Int( Int( numer(%),X[2]=0..infinity ),X[3]=-infinity..infinity );
  value(%);

The obtained expression results in:

> simplify( -1/4*a(X[1])^3*Pi*tau(X[1])^4*(2*W(X[1])*`W*`(X[1])*b(X[1])*a(X[1])-W(X[1])*df(`W*`(X[1]),X[1])*a(X[1])*k+2*`W*`(X[1])*D*W(X[1])*a(X[1])*k*psi(X[1])-`W*`(X[1])*df(W(X[1]),X[1])*a(X[1])*k-`W*`(X[1])*W(X[1])*k*df(a(X[1]),X[1]))*2^(1/2)+1/4*a(X[1])^4*Pi*`W*`(X[1])*W(X[1])*k*(2*D*tau(X[1])*psi(X[1])+df(tau(X[1]),X[1]))*2^(1/2)*tau(X[1])^3 ):

expand(%/(1/4*a(X[1])^3*Pi*tau(X[1])^3*sqrt(2)*k));

This gives:

> eq16 := Diff(`W*`(X[1])*W(X[1])*a(X[1])*tau(X[1]),X[1]) = 2/k*tau(X[1])*W(X[1])*`W*`(X[1])*b(X[1])*a(X[1]);

The spatial first-order momentum is:

>    # T[1,0]

subs( {eq6,eq7},lhs(eq1) - rhs(eq1) )*rhs(eq7)*X[2] + subs( {eq6,eq7},lhs(eq2) - rhs(eq2) )*rhs(eq6)*X[2]:
simplify( value(%) ):
  Int( Int( numer(%),X[2]=0..infinity ),X[3]=-infinity..infinity );
  value(%);

> simplify( 1/4*k*a(X[1])^4*tau(X[1])^4*(2*`W*`(X[1])*W(X[1])*df(a(X[1]),X[1])-2*`W*`(X[1])*D*W(X[1])*a(X[1])*psi(X[1])+W(X[1])*df(`W*`(X[1]),X[1])*a(X[1])+`W*`(X[1])*df(W(X[1]),X[1])*a(X[1]))*2^(1/2)*Pi^(1/2)+1/4*k*a(X[1])^5*W(X[1])*`W*`(X[1])*(df(tau(X[1]),X[1])+2*D*tau(X[1])*psi(X[1]))*2^(1/2)*tau(X[1])^3*Pi^(1/2) ):

expand(%/(1/4*k*a(X[1])^3*tau(X[1])^3*sqrt(2*Pi)));

And we have the energy conservation law:

> eq17 := Diff(`W*`(X[1])*W(X[1])*a(X[1])^2*tau(X[1]),X[1]) = 0;# this gives eq[1]

In the combination with eq16 this gives the equation for the beam size:

> simplify( value(eq16 - eq17/a(X[1])) );
print(`This is eq[3]`);
df(a(X[1]),X[1]) = solve(%,df(a(X[1]),X[1]));

The spatial second-order momentum is:

>    # T[2,0]

subs( {eq6,eq7},lhs(eq1) - rhs(eq1) )*rhs(eq7)*X[2]^2 + subs( {eq6,eq7},lhs(eq2) - rhs(eq2) )*rhs(eq6)*X[2]^2:
simplify( value(%) ):
  Int( Int( numer(%),X[2]=0..infinity ),X[3]=-infinity..infinity );
  value(%);

> simplify(
-1/8*a(X[1])^5*Pi*tau(X[1])^4*(-2*`W*`(X[1])*W(X[1])*b(X[1])*a(X[1])-W(X[1])*df(`W*`(X[1]),X[1])*a(X[1])*k+2*`W*`(X[1])*D*W(X[1])*a(X[1])*k*psi(X[1])-`W*`(X[1])*df(W(X[1]),X[1])*a(X[1])*k-3*`W*`(X[1])*W(X[1])*k*df(a(X[1]),X[1]))*2^(1/2)+1/8*a(X[1])^6*Pi*`W*`(X[1])*W(X[1])*k*(2*D*tau(X[1])*psi(X[1])+df(tau(X[1]),X[1]))*2^(1/2)*tau(X[1])^3 ):

eq18 := expand(%/(1/8*a(X[1])^5*Pi*tau(X[1])^3*sqrt(2)*k));

This is the energy conservation law:

> value( eq16 ):
expand(a(X[1])*(lhs(%) - rhs(%) + eq18)/2);# this is energy conservation law, too

In fact, the expression for the beam size allows transforming eq16 into the energy conservation law:

> subs( b(X[1]) = -k*df(a(X[1]),X[1])/2,value(eq16) );# energy conservation law

Now let us consider the time-momentums. =0. However, gives:

>    # T[0,2]

subs( {eq6,eq7},lhs(eq1) - rhs(eq1) )*rhs(eq7)*X[3]^2 + subs( {eq6,eq7},lhs(eq2) - rhs(eq2) )*rhs(eq6)*X[3]^2:
simplify( value(%) ):
  Int( Int( numer(%),X[2]=0..infinity ),X[3]=-infinity..infinity );
  value(%);

Hence

> simplify( -1/16*a(X[1])^3*Pi*tau(X[1])^6*(2*W(X[1])*`W*`(X[1])*b(X[1])*a(X[1])-W(X[1])*df(`W*`(X[1]),X[1])*a(X[1])*k+2*`W*`(X[1])*D*W(X[1])*a(X[1])*k*psi(X[1])-`W*`(X[1])*df(W(X[1]),X[1])*a(X[1])*k-`W*`(X[1])*W(X[1])*k*df(a(X[1]),X[1]))*2^(1/2)+3/16*a(X[1])^4*Pi*`W*`(X[1])*W(X[1])*k*(2*D*tau(X[1])*psi(X[1])+df(tau(X[1]),X[1]))*2^(1/2)*tau(X[1])^5 ):

eq19 := expand(%/(1/16*a(X[1])^3*Pi*tau(X[1])^5*sqrt(2)*k));

In the combination with the energy conservation law this leads to the equation for the pulse duration:

>    value( eq16 ):
expand(a(X[1])*(lhs(%) - rhs(%) - eq19)/2);
   print(`This is eq[2]`);
    df(tau(X[1]),X[1]) = solve(%,df(tau(X[1]),X[1]));

To obtain the equations for chirp and wave-front curvature we need and -momentums.

> # J-momentum

eq20 := Diff( (df(A1(X),X[3])*A(X) - df(A(X),X[3])*A1(X)),X[1] );
eq21 := df(A(X),X[1])=subs( abs(A(X))^2=A(X)*A1(X),expand(solve(eq1,df(A(X),X[1]))) );
eq22 := df(A1(X),X[1])=subs( abs(A(X))^2=A(X)*A1(X),expand(solve(eq2,df(A1(X),X[1]))) );

subs({eq21,eq22},value(eq20)):
eq23 := expand( simplify(%) );

eq20 gives J -momentum and eq23 gives its rate of change. The substitution of the trial functions eq6, eq7 for the pulse profile results in:

>    subs({eq6,eq7},eq23):
simplify(%):
  eq24 := expand(%);

subs({eq6,eq7},eq20):
simplify(value(%)):
  eq25 := expand(%);

> eq26 := expand((tau(x1)^4*a(x1)^3*k)*(eq25 - eq24)/4/i);

Since let's try :

> # this group requires a double evoluation!

simplify(int(eq26*X[3],X[2]=0..infinity)):
int(%,X[3]=-infinity..infinity):
eq27 := expand( simplify(%)*16/tau(X[1])^4/Pi/sqrt(2)/csgn(conjugate(tau(X[1])))/abs(a(X[1]))^3/k );

Hence and from the energy conservation law we have:

>    expand( subs({eq[2],eq[3]}, value(lhs(eq17)) )/tau(X[1])/a(X[1])^2 ):# energy conservation law
expand( solve(%=0,`W*`(X[1])*df(W(X[1]),X[1])) ):

expand( 2*k*simplify( subs({eq[2],eq[3]},eq27) )/a(X[1])/k/2/tau(X[1])^3/psi(X[1]) ):# from eq27
simplify( subs(`W*`(X[1])*df(W(X[1]),X[1])=%%,%) ):
  numer(%)/`W*`(X[1])/W(X[1]) = 0;
   df(psi(X[1]),X[1]) = collect( expand( solve(%,df(psi(X[1]),X[1])) ),D );
    eq[4];# from the variational method

We face with the difference of the numerical coefficients before the nonlinear term:

> evalf(sqrt(2)/4/(1/2));# there exists some difference in the coefficients before nonlinear terms in both approaches

Such difference requires the numerical investigation, which will be given elsewhere. However, we can suppose that the difference can result from the Gaussian approximation for the time profile of the pulse. Probably, we need sech -form of -dependence for the trial functions . This guess is suggested also by the difference of the coefficients obtained for the Gaussian approximation and for the time-independent beam in the framework of the variation method only.

Now let us consider -momentum:

> eq28 := Diff( (df(A1(X),X[2])*A(X) - df(A(X),X[2])*A1(X)),X[1] );
eq29 := df(A(X),X[1])=subs( abs(A(X))^2=A(X)*A1(X),expand(solve(eq1,df(A(X),X[1]))) );
eq30 := df(A1(X),X[1])=subs( abs(A(X))^2=A(X)*A1(X),expand(solve(eq2,df(A1(X),X[1]))) );

subs({eq29,eq30},value(eq28)):
eq31 := expand( simplify(%) );

For the considered trial functions we have:

>    subs({eq6,eq7},eq31):
simplify(%):
  eq32 := expand(%);

subs({eq6,eq7},eq28):
simplify(value(%)):
  eq33 := expand(%);

And -momentum is:

> eq34 := expand( a(X[1])^4*simplify((eq33 - eq32))/2/i );

> simplify(int(eq34,X[2]=0..infinity)):
int(%,X[3]=-infinity..infinity);

Hence

> simplify( 1/4*W(X[1])^2*`W*`(X[1])^2*kappa*a(X[1])^4*tau(X[1])*Pi^(1/2) + 1/4*a(X[1])^2*tau(X[1])*(-2*`W*`(X[1])*df(W(X[1]),X[1])*b(X[1])*a(X[1])^4*k*wp^2-4*a(X[1])^3*`W*`(X[1])*W(X[1])*b(X[1])*df(a(X[1]),X[1])*k*wp^2+4*W(X[1])*`W*`(X[1])*K*a(X[1])^4*k-2*`W*`(X[1])*W(X[1])*df(b(X[1]),X[1])*a(X[1])^4*k*wp^2-2*df(`W*`(X[1]),X[1])*W(X[1])*b(X[1])*a(X[1])^4*k*wp^2-W(X[1])*`W*`(X[1])*wp^2+4*a(X[1])^4*W(X[1])*`W*`(X[1])*b(X[1])^2*wp^2+4*`W*`(X[1])*W(X[1])*D*b(X[1])*psi(X[1])*a(X[1])^4*k*wp^2)*2^(1/2)*Pi^(1/2)/k/wp^2-1/2*a(X[1])^6*`W*`(X[1])*W(X[1])*b(X[1])*(df(tau(X[1]),X[1])+2*D*tau(X[1])*psi(X[1]))*2^(1/2)*Pi^(1/2) ):

eq35 := expand( 2*%/a(X[1])^4/sqrt(2*Pi)/b(X[1]) );

In the combination with the energy conservation law this results in:

>    lhs( value(eq17) ):
simplify(% + eq35):
  expand(b(X[1])*%/(tau(X[1])*W(X[1])*`W*`(X[1]))):
   df(b(X[1]),X[1])=expand( solve(%=0,df(b(X[1]),X[1])) );
    eq[5];# from the variational method

Here also there exists some numerical difference of the coefficients before the nonlinear term:

> evalf(sqrt(2)/8/(sqrt(2)/4));

Dissipative system

Now let us consider the ultrashort pulse propagation within the real active medium possessing the radially varying gain (g is the saturated gain coefficient of the active medium, G is the gain of the propagating pulse) and the gain band ( is the inverse gain bandwidth):

>    assume(g,real):
assume(tg,real):
assume(wp,real):
assume(wp,positive):

eq36 := lhs(eq1) = subs(abs(A(X))^2=A(X)*A1(X),rhs(eq1)) + g*exp(-2*X[2]^2/wp^2)*A(X) + tg^2*df(A(X),X[3]$2);
eq37 := lhs(eq2) = subs(abs(A(X))^2=A(X)*A1(X),rhs(eq2)) + g*exp(-2*X[2]^2/wp^2)*A1(X) + tg^2*df(A1(X),X[3]$2);
  eq38 := A(X) = W(X[1])*exp( G(X[1]) + i*b(X[1])*X[2]^2 - X[2]^2/a(X[1])^2/2 + i*psi(X[1])*X[3]^2 - X[3]^2/tau(X[1])^2 );# Gaussian pulse with complex amplitude W, inverse square of the wavefront curvature radius b, beam size a, pulse width tau and chirp psi
   eq39 := A1(X) = `W*`(X[1])*exp( G(X[1]) - i*b(X[1])*X[2]^2 - X[2]^2/a(X[1])^2/2 - i*psi(X[1])*X[3]^2 - X[3]^2/tau(X[1])^2);# complex conjugate solution

-momentum is:

>    # T[0,0]

subs( {eq6,eq7}, lhs(eq36) - rhs(eq36) )*rhs(eq7) + subs( {eq6,eq7}, lhs(eq37) - rhs(eq37) )*rhs(eq6):
simplify( value(%) ):
  Int( Int( numer(%),X[2]=0..infinity ),X[3]=-infinity..infinity ):
   value(%):
    value( subs({W(X[1])=W(X[1])*exp(G(X[1])),`W*`(X[1])=`W*`(X[1])*exp(G(X[1]))},%) );# from eq38, eq39

The last is:

> eq40 := expand( 4*( -1/4*a(X[1])^3*Pi*tau(X[1])^3*(-W(X[1])*exp(G(X[1]))*(df(`W*`(X[1]),X[1])*exp(G(X[1]))+`W*`(X[1])*df(G(X[1]),X[1])*exp(G(X[1])))*tau(X[1])^2*k*(wp^2+2*a(X[1])^2)^(1/2)*a(X[1])-4*`W*`(X[1])*exp(G(X[1]))^2*tg^2*W(X[1])*k*(wp^2+2*a(X[1])^2)^(1/2)*a(X[1])+2*`W*`(X[1])*exp(G(X[1]))^2*W(X[1])*tau(X[1])^2*b(X[1])*(wp^2+2*a(X[1])^2)^(1/2)*a(X[1])+2*`W*`(X[1])*exp(G(X[1]))^2*g*W(X[1])*tau(X[1])^2*k*wp*a(X[1])+2*`W*`(X[1])*exp(G(X[1]))^2*D*W(X[1])*k*psi(X[1])*tau(X[1])^2*(wp^2+2*a(X[1])^2)^(1/2)*a(X[1])-`W*`(X[1])*exp(G(X[1]))^2*W(X[1])*tau(X[1])^2*k*df(a(X[1]),X[1])*(wp^2+2*a(X[1])^2)^(1/2)-`W*`(X[1])*exp(G(X[1]))*(df(W(X[1]),X[1])*exp(G(X[1]))+W(X[1])*df(G(X[1]),X[1])*exp(G(X[1])))*tau(X[1])^2*k*(wp^2+2*a(X[1])^2)^(1/2)*a(X[1]))*2^(1/2)/(wp^2+2*a(X[1])^2)^(1/2)+1/4*2^(1/2)*a(X[1])^4*Pi*`W*`(X[1])*exp(G(X[1]))^2*W(X[1])*k*(2*D*psi(X[1])*tau(X[1])^2+tau(X[1])*df(tau(X[1]),X[1])-2*tg^2+2*tg^2*psi(X[1])^2*tau(X[1])^4)*tau(X[1])^3 ) /exp(2*G(X[1]))/sqrt(2)/Pi/tau(X[1])^4/a(X[1])^3/k );

Hence

> eq41 := Diff(W(X[1])*`W*`(X[1])*a(X[1])*tau(X[1]),X[1]) = -expand(eq40 - df(W(X[1])*`W*`(X[1])*a(X[1])*tau(X[1]),X[1]));

-momentum is:

>    # T[1,0]

subs( {eq6,eq7}, lhs(eq36) - rhs(eq36) )*rhs(eq7)*X[2] + subs( {eq6,eq7}, lhs(eq37) - rhs(eq37) )*rhs(eq6)*X[2]:
simplify( value(%) ):
  Int( Int( numer(%),X[2]=0..infinity ),X[3]=-infinity..infinity ):
   value(%):
    value( subs({W(X[1])=W(X[1])*exp(G(X[1])),`W*`(X[1])=`W*`(X[1])*exp(G(X[1]))},%) );# from eq38, eq39

Thus, we have:

> eq42 := expand( numer( simplify( -1/4*a(X[1])^4*k*tau(X[1])^3*(-a(X[1])*`W*`(X[1])*exp(G(X[1]))*(df(W(X[1]),X[1])*exp(G(X[1]))+W(X[1])*df(G(X[1]),X[1])*exp(G(X[1])))*tau(X[1])^2*wp^2-2*`W*`(X[1])*exp(G(X[1]))^2*W(X[1])*tau(X[1])^2*df(a(X[1]),X[1])*wp^2+4*a(X[1])^3*`W*`(X[1])*exp(G(X[1]))^2*D*W(X[1])*psi(X[1])*tau(X[1])^2+2*a(X[1])*`W*`(X[1])*exp(G(X[1]))^2*D*W(X[1])*psi(X[1])*tau(X[1])^2*wp^2+2*a(X[1])*`W*`(X[1])*exp(G(X[1]))^2*g*W(X[1])*tau(X[1])^2*wp^2-4*a(X[1])^2*`W*`(X[1])*exp(G(X[1]))^2*W(X[1])*tau(X[1])^2*df(a(X[1]),X[1])-2*a(X[1])^3*`W*`(X[1])*exp(G(X[1]))*(df(W(X[1]),X[1])*exp(G(X[1]))+W(X[1])*df(G(X[1]),X[1])*exp(G(X[1])))*tau(X[1])^2-a(X[1])*W(X[1])*exp(G(X[1]))*(df(`W*`(X[1]),X[1])*exp(G(X[1]))+`W*`(X[1])*df(G(X[1]),X[1])*exp(G(X[1])))*tau(X[1])^2*wp^2-8*a(X[1])^3*`W*`(X[1])*exp(G(X[1]))^2*tg^2*W(X[1])-4*a(X[1])*`W*`(X[1])*exp(G(X[1]))^2*tg^2*W(X[1])*wp^2-2*a(X[1])^3*W(X[1])*exp(G(X[1]))*(df(`W*`(X[1]),X[1])*exp(G(X[1]))+`W*`(X[1])*df(G(X[1]),X[1])*exp(G(X[1])))*tau(X[1])^2)*2^(1/2)*Pi^(1/2)/(wp^2+2*a(X[1])^2)+1/4*Pi^(1/2)*tau(X[1])^3*2^(1/2)*(2*D*psi(X[1])*tau(X[1])^2+tau(X[1])*df(tau(X[1]),X[1])-2*tg^2+2*tg^2*psi(X[1])^2*tau(X[1])^4)*W(X[1])*exp(G(X[1]))^2*`W*`(X[1])*k*a(X[1])^5 ))/k/sqrt(2*Pi)/exp(2*G(X[1]))/a(X[1])^3/tau(X[1])^4/wp^2 );

Hence (the law of the energy change):

> eq43 := Diff( `W*`(X[1])*W(X[1])*a(X[1])^2*tau(X[1]),X[1] ) = -expand(eq42 - (df( `W*`(X[1])*W(X[1])*a(X[1])^2*tau(X[1]),X[1] ) + 2*a(X[1])^2*df(`W*`(X[1])*W(X[1])*a(X[1])^2*tau(X[1]),X[1])/wp^2 ))/(1 + 2*a(X[1])^2/wp^2);

Eqs. 41,43 allow the following simplification:

>    expand( value( lhs(eq43) )/a(X[1]) ):
value( lhs(eq41) ):
  expand(%% - %):
   simplify( rhs(eq43)/a(X[1]) - rhs(eq41) ):
    %% = %:
     expand( solve(%,df(a(X[1]),X[1])) ):
    eq44 := df(a(X[1]),X[1]) = collect(%,{g,b(X[1])});

The right-hand side can be transformed as:

> simplify(op(2,rhs(eq44)));
simplify(op(1,rhs(eq44)),symbolic);
simplify(% - (-2*g*a(X[1])*(sqrt(1+2*a(X[1])^2/wp^2)-1)/(1+2*a(X[1])^2/wp^2)));

Hence:

> sol[1] := df(a(X[1]),X[1]) = -2*g*a(X[1])*(sqrt(1+2*a(X[1])^2/wp^2)-1)/(1+2*a(X[1])^2/wp^2) - 2*a(X[1])*b(X[1])/k;

The next figure demonstrates that in our case (red curve) the soft aperture induced by the radially varying gain is twice as that in Ref. 5 (green curve) and decays more slowly with the growth of the laser beam size.

>    with(plots):

print(`increment of the soft aperture (our result):`);
op(1,rhs(sol[1]))/a(X[1])/g;
print(`increment of the soft aperture (Ref. 5):`);
-(wp^2/a(X[1])^2)/(1+wp^2/a(X[1])^2/2)^2/2;

# wp = 100 mkm

p1 := plot(subs({wp=1e-2,a(X[1])=x},%%),x=0..1e-1,color=red):
p2 := plot(subs({wp=1e-2,a(X[1])=x},%%),x=0..1e-1,color=green):
  display(p1,p2,axes=boxed,title=`increment of the soft aperture vs. laser beam size`);

Warning, the name changecoords has been redefined

The analysis of the next moment smoothes away this disagreement and gains the difference between an unambiguity of the variation approach and some ambiguity of the momentums method. We suppose that the last is rooted in the nonintegrability of the considered dissipative equation.

>    # T[2,0]

subs( {eq6,eq7}, lhs(eq36) - rhs(eq36) )*rhs(eq7)*X[2]^2 + subs( {eq6,eq7}, lhs(eq37) - rhs(eq37) )*rhs(eq6)*X[2]^2:
simplify( value(%) ):
  Int( Int( numer(%),X[2]=0..infinity ),X[3]=-infinity..infinity ):
   value(%):
    value( subs({W(X[1])=W(X[1])*exp(G(X[1])),`W*`(X[1])=`W*`(X[1])*exp(G(X[1]))},%) );# from eq38, eq39

Hence:

> numer( simplify( -1/8*a(X[1])^5*Pi*tau(X[1])^3*(-W(X[1])*exp(G(X[1]))*(df(`W*`(X[1]),X[1])*exp(G(X[1]))+`W*`(X[1])*df(G(X[1]),X[1])*exp(G(X[1])))*tau(X[1])^2*k*(wp^2+2*a(X[1])^2)^(1/2)*wp^2*a(X[1])-2*a(X[1])^3*W(X[1])*exp(G(X[1]))*(df(`W*`(X[1]),X[1])*exp(G(X[1]))+`W*`(X[1])*df(G(X[1]),X[1])*exp(G(X[1])))*tau(X[1])^2*k*(wp^2+2*a(X[1])^2)^(1/2)-4*`W*`(X[1])*exp(G(X[1]))^2*tg^2*W(X[1])*k*(wp^2+2*a(X[1])^2)^(1/2)*wp^2*a(X[1])-8*a(X[1])^3*`W*`(X[1])*exp(G(X[1]))^2*tg^2*W(X[1])*k*(wp^2+2*a(X[1])^2)^(1/2)-3*`W*`(X[1])*exp(G(X[1]))^2*W(X[1])*tau(X[1])^2*k*df(a(X[1]),X[1])*(wp^2+2*a(X[1])^2)^(1/2)*wp^2-6*`W*`(X[1])*exp(G(X[1]))^2*W(X[1])*tau(X[1])^2*k*df(a(X[1]),X[1])*a(X[1])^2*(wp^2+2*a(X[1])^2)^(1/2)+2*wp^3*`W*`(X[1])*exp(G(X[1]))^2*g*W(X[1])*tau(X[1])^2*k*a(X[1])-2*a(X[1])^3*`W*`(X[1])*exp(G(X[1]))*(df(W(X[1]),X[1])*exp(G(X[1]))+W(X[1])*df(G(X[1]),X[1])*exp(G(X[1])))*tau(X[1])^2*k*(wp^2+2*a(X[1])^2)^(1/2)-4*a(X[1])^3*`W*`(X[1])*exp(G(X[1]))^2*W(X[1])*tau(X[1])^2*b(X[1])*(wp^2+2*a(X[1])^2)^(1/2)-2*`W*`(X[1])*exp(G(X[1]))^2*W(X[1])*tau(X[1])^2*b(X[1])*(wp^2+2*a(X[1])^2)^(1/2)*wp^2*a(X[1])+4*a(X[1])^3*`W*`(X[1])*exp(G(X[1]))^2*D*W(X[1])*k*psi(X[1])*tau(X[1])^2*(wp^2+2*a(X[1])^2)^(1/2)+2*`W*`(X[1])*exp(G(X[1]))^2*D*W(X[1])*k*psi(X[1])*tau(X[1])^2*(wp^2+2*a(X[1])^2)^(1/2)*wp^2*a(X[1])-`W*`(x1)*exp(G(X[1]))*(df(W(X[1]),X[1])*exp(G(X[1]))+W(X[1])*df(G(X[1]),X[1])*exp(G(X[1])))*tau(X[1])^2*k*(wp^2+2*a(X[1])^2)^(1/2)*wp^2*a(X[1]))*2^(1/2)/(wp^2+2*a(X[1])^2)^(3/2)+1/8*tau(X[1])^3*(2*D*psi(X[1])*tau(X[1])^2+tau(X[1])*df(tau(X[1]),X[1])-2*tg^2+2*tg^2*psi(X[1])^2*tau(X[1])^4)*k*W(X[1])*exp(G(X[1]))^2*`W*`(X[1])*Pi*a(X[1])^6*2^(1/2) )*8/(a(X[1])^5*Pi*tau(X[1])^3*exp(2*G(X[1]))*sqrt(2)) ):

eq45 := expand(%/sqrt(wp^2+2*a(X[1])^2)/wp^2/k);

Combination with eq42 results in:

>    simplify(eq45 - tau(X[1])*eq42/a(X[1])):
expand( numer(%)/(`W*`(X[1])*W(X[1])*tau(X[1])^2) );
  collect( expand( solve(%=0,df(a(X[1]),X[1])) ),{g,b(X[1])} ):
   simplify(op(2,%));
    simplify(op(1,%%),symbolic);

Thus

> sol[2] := df(a(X[1]),X[1]) = -2*g*a(X[1])*(sqrt(1+2*a(X[1])^2/wp^2) - 1)/(1+2*a(X[1])^2/wp^2)^(3/2) -2*a(X[1])*b(X[1])/k;

>    op(1,rhs(sol[2]))/a(X[1])/g;# new approximation
op(1,rhs(sol[1]))/a(X[1])/g;# previous approximation

p1 := plot(subs({wp=1e-2,a(X[1])=x},%%),x=0..1e-1,color=blue):
p2 := plot(subs({wp=1e-2,a(X[1])=x},%%),x=0..1e-1,color=red):
-(wp^2/a(X[1])^2)/(1+wp^2/a(X[1])^2/2)^2/2;# Ref. 5
  p3 := plot(subs({wp=1e-2,a(X[1])=x},%),x=0..1e-1,color=green):
  display(p1,p2,p3,axes=boxed,title=`increment of the soft aperture vs. laser beam size`);

Thus, we approached to the solution of Ref. 5. However, we can not consider this law as the limit corresponding to the real law for beam size change. For example, eq41 in the combination with eq45 produce the next approximation:

> eq46 := tau(X[1])*W(X[1])*df(`W*`(X[1]),X[1])*a(X[1]) + tau(X[1])*`W*`(X[1])*df(W(X[1]),X[1])*a(X[1]) + `W*`(X[1])*W(X[1])*a(X[1])*df(tau(X[1]),X[1]) = rhs(eq41) - tau(X[1])*W(X[1])*`W*`(X[1])*df(a(X[1]),X[1]);# from eq41

>    eq47 := expand(eq45/tau(X[1])):
simplify(eq47 - (lhs(eq46) - rhs(eq46))):
  eq48 := expand(numer(%)/2/tau(X[1])/a(X[1])^2/k):
   simplify(eq48 - (lhs(eq46) - rhs(eq46))):
    expand(numer(%)/(tau(X[1])*`W*`(X[1])*W(X[1]))):
     collect( expand( solve(%=0,df(a(X[1]),X[1])) ),{g,b(X[1])} ):
      simplify(op(2,%)):
       simplify(op(1,%%),symbolic):
sol[3] := df(a(X[1]),X[1]) = %% + %;

> op(2,rhs(sol[3]))/a(X[1])/g;# new approximation
p1 := plot(subs({wp=1e-2,a(X[1])=x},%),x=0..1e-1,color=black):
op(1,rhs(sol[2]))/a(X[1])/g;# previous approximation
p2 := plot(subs({wp=1e-2,a(X[1])=x},%),x=0..1e-1,color=red):
op(1,rhs(sol[1]))/a(X[1])/g;# previous approximation
p3 := plot(subs({wp=1e-2,a(X[1])=x},%),x=0..1e-1,color=red):
-(wp^2/a(X[1])^2)/(1+wp^2/a(X[1])^2/2)^2/2;# Ref. 5
p4 := plot(subs({wp=1e-2,a(X[1])=x},%),x=0..1e-1,color=green):

display(p1,p2,p3,p4,axes=boxed,title=`increment of the soft aperture vs. laser beam size`);

Later on we shall use in our analysis.

Now let us consider the next momentum ( ).

>    # T[0,2]

subs( {eq6,eq7}, lhs(eq36) - rhs(eq36) )*rhs(eq7)*X[3]^2 + subs( {eq6,eq7}, lhs(eq37) - rhs(eq37) )*rhs(eq6)*X[3]^2:
simplify( value(%) ):
  Int( Int( numer(%),X[2]=0..infinity ),X[3]=-infinity..infinity ):
   value(%):
    value( subs({W(X[1])=W(X[1])*exp(G(X[1])),`W*`(X[1])=`W*`(X[1])*exp(G(X[1]))},%) );# from eq38, eq39

> simplify( -1/16*a(X[1])^3*Pi*tau(X[1])^5*(-W(X[1])*exp(G(X[1]))*(df(`W*`(X[1]),X[1])*exp(G(X[1]))+`W*`(X[1])*df(G(X[1]),X[1])*exp(G(X[1])))*tau(X[1])^2*k*(wp^2+2*a(X[1])^2)^(1/2)*a(X[1])-4*`W*`(X[1])*exp(G(X[1]))^2*tg^2*W(X[1])*k*(wp^2+2*a(X[1])^2)^(1/2)*a(X[1])+2*`W*`(X[1])*exp(G(X[1]))^2*W(X[1])*tau(X[1])^2*b(X[1])*(wp^2+2*a(X[1])^2)^(1/2)*a(X[1])+2*`W*`(X[1])*exp(G(X[1]))^2*g*W(X[1])*tau(X[1])^2*k*wp*a(X[1])+2*`W*`(X[1])*exp(G(X[1]))^2*D*W(X[1])*k*psi(X[1])*tau(X[1])^2*(wp^2+2*a(X[1])^2)^(1/2)*a(X[1])-`W*`(X[1])*exp(G(X[1]))^2*W(X[1])*tau(X[1])^2*k*df(a(X[1]),X[1])*(wp^2+2*a(X[1])^2)^(1/2)-`W*`(X[1])*exp(G(X[1]))*(df(W(X[1]),X[1])*exp(G(X[1]))+W(X[1])*df(G(X[1]),X[1])*exp(G(X[1])))*tau(X[1])^2*k*(wp^2+2*a(X[1])^2)^(1/2)*a(X[1]))*2^(1/2)/(wp^2+2*a(X[1])^2)^(1/2)+3/16*2^(1/2)*a(X[1])^4*Pi*`W*`(X[1])*exp(G(X[1]))^2*W(X[1])*k*(2*D*psi(X[1])*tau(X[1])^2+tau(X[1])*df(tau(X[1]),X[1])-2*tg^2+2*tg^2*psi(X[1])^2*tau(X[1])^4)*tau(X[1])^5 ):

eq49 := expand(numer(16*%)/(wp^2*k*a(X[1])^3*tau(X[1])^6*exp(2*G(X[1]))*2^(1/2)*Pi)*sqrt(wp^2+2*a(X[1])^2));

In the combination with eq40 we have the solution for the pulse width:

>    expand(wp^2*(eq49 - eq40)/2/a(X[1])^2):
simplify(% - eq40):
  expand( numer(%)/(`W*`(X[1])*W(X[1])) ):
   sol[4] := df(tau(X[1]),X[1]) = collect( expand( solve(%=0,df(tau(X[1]),X[1])) ),tg^2 );

Note, that eq41 gives the law for the amplitude change:

> expand( value( eq41 )/(tau(X[1])*a(X[1])) ):

lhs(%) - 1/a(X[1])*`W*`(X[1])*W(X[1])*df(a(X[1]),X[1]) - 1/tau(X[1])*`W*`(X[1])*W(X[1])*df(tau(X[1]),X[1]) = rhs(%) - 1/a(X[1])*`W*`(X[1])*W(X[1])*df(a(X[1]),X[1]) - 1/tau(X[1])*`W*`(X[1])*W(X[1])*df(tau(X[1]),X[1]);

Diff(abs(W(X[1]))^2,X[1]) + 2*abs(W(X[1]))^2*df(G(X[1]),X[1]) = expand( subs({sol[1],sol[4]},rhs(%)) + 2*W(X[1])*`W*`(X[1])*df(G(X[1]),X[1]));

Hence

> eq50 := Diff(abs(W(X[1]))^2,X[1]) + 2*abs(W(X[1]))^2*df(G(X[1]),X[1]) = -4*tg^2*abs(W(X[1]))^2/tau(X[1])^2 + 4*abs(W(X[1]))^2*b(X[1])/k + 4*g*abs(W(X[1]))^2/sqrt(1+2*a(X[1])^2/wp^2) - 2*g*abs(W(X[1]))^2/(1+2*a(X[1])^2/wp^2) + 2*D*psi(X[1])*abs(W(X[1]))^2;

If we take into account the geometrical "magnification" of the field intensity, the law of the gain change can be expressed as:

>    value( subs(abs(W(X[1]))^2=abs(W0)^2*a0^2/a(X[1])^2,eq50) ):
subs(sol[1],%):
  simplify(lhs(%) - rhs(%)):
   expand( numer(%)/(2*abs(W0)^2*a0^2) ):
    solve(%=0,df(G(X[1]),X[1])):
     expand(%);

> sol[5] := df(G(X[1]),X[1]) = g/(1 + 2*a(X[1])^2/wp^2) - 2*tg^2/tau(X[1])^2 + D*psi(X[1]);

Now we turn to the M -momentums.

> eq51 := Diff( (df(A1(X),X[2])*A(X) - df(A(X),X[2])*A1(X)),X[1] );
eq52 := df(A(X),X[1])=subs( abs(A(X))^2=A(X)*A1(X),expand(solve(eq36,df(A(X),X[1]))) );
eq53 := df(A1(X),X[1])=subs( abs(A(X))^2=A(X)*A1(X),expand(solve(eq37,df(A1(X),X[1]))) );

subs({eq52,eq53},value(eq51)):
eq54 := expand( simplify(%) );

>    subs({eq38,eq39},eq54):
simplify(%):
  eq55 := expand(%);

subs({eq38,eq39},eq51):
simplify(value(%)):
  eq56 := expand(%);

Hence -momentum is:

> eq57 := expand( k*a(X[1])^4*tau(X[1])^4*simplify((eq56 - eq55))/2/i );

> # this group requires the double evoluation!

simplify(int(eq57,X[2]=0..infinity)):
simplify( int(%,X[3]=-infinity..infinity) ):
eq58 := simplify( 4*numer(%)/(a(X[1])^2*tau(X[1])^2*exp(2*G(X[1])))*Pi );

> eq59 := expand( -eq58/4/tau(X[1])/Pi^(3/2)/csgn(conjugate(tau(X[1])))/wp^4/sqrt(2)/b(X[1])/2/a(X[1])^2/k/tau(X[1]) );

In the combination with eq42 this results in:

> simplify(eq59 - eq42):
expand( -4*numer(%)/tau(X[1])/(`W*`(X[1])*W(X[1])) ):
expand( solve(%=0,df(b(X[1]),X[1])) );

Hence (taking into account the geometrical magnification):

> sol[6] := df(b(X[1]),X[1]) = 2/wp^2*K+2/k*b(X[1])^2+1/4*1/a(X[1])^4*2^(1/2)*kappa*abs(W0)^2*a0^2*exp(G(X[1]))^2-1/2*1/(a(X[1])^4*k);

Now, as it was done above, we analyze J -momentum:

> # J-momentum

eq60 := Diff( (df(A1(X),X[3])*A(X) - df(A(X),X[3])*A1(X)),X[1] );
eq61 := df(A(X),X[1])=subs( abs(A(X))^2=A(X)*A1(X),expand(solve(eq36,df(A(X),X[1]))) );
eq62 := df(A1(X),X[1])=subs( abs(A(X))^2=A(X)*A1(X),expand(solve(eq37,df(A1(X),X[1]))) );

subs({eq61,eq62},value(eq60)):
eq63 := expand( simplify(%) );

>    subs({eq38,eq39},eq63):
simplify(%):
  eq64 := expand(%);

subs({eq38,eq39},eq60):
simplify(value(%)):
  eq65 := expand(%);

> eq66 := expand((tau(x1)^4*a(x1)^3*k)*(eq65 - eq64)/4/i);

-momentum is:

> simplify(int(eq66*X[3],X[2]=0..infinity)):
int(%,X[3]=-infinity..infinity):
eq67 := numer(simplify(%));

> eq68 := expand( eq67/(-2*k*psi(X[1])*sqrt(wp^2+2*a(x1)^2)*exp(2*G(X[1]))*tau(X[1])^6*2^(1/2)*abs(a(X[1]))^3*Pi*csgn(conjugate(tau(X[1])))) );

and in the combination with eq40 we have

> expand( 2*numer( simplify( eq68 - eq40 ) )/(W(X[1])*`W*`(X[1])*a(X[1])) ):
simplify( subs(sol[4],%) ):
expand( solve( %=0,df(psi(X[1]),X[1]) ) );

sol[7] := df(psi(X[1]),X[1]) = 2*D*(psi(X[1])^2 - 1/1/tau(X[1])^4) - 8*psi(X[1])/tau(X[1])^2*tg^2 + 1/2*1/tau(X[1])^3*abs(W0)^2*a0^2*tau0*kappa*exp(G(X[1]))^2/a(X[1])^2;

Thus, we have the final system , which describes the ultrashort pulse propagation within the gain medium:

>    sol[3];
  sol[4];
   sol[5];
    sol[6];
     sol[7];
      abs(W(X[1]))^2 = abs(W0)^2*a0^2*exp(2*G(X[1]))/a(X[1])^2;

Conclusion

Thus, the analysis of the ultrashort pulse propagation within the gain medium can be reduced to the solution of the first-order differential equations without parabolic approximations and expansions in the series of the time- and spatil-variables. This allows reducing the analysis of the real-world Kerr-lens mode-locked lasers to the simple iteration procedure based on the above obtained equations (for the intra-medium propagation) and usual ABCD-law (for the extra-medium propagation). Also, we can take into account the frequency filtering on the out-put mirror and the real configuration of the intra-cavity group-velocity dispersion. This approach can be generalized on the systems possessing the beam astigmatism. Moreover, we can analyze the competition of the continuous-wave and pulsing regimes (including multi-pulsing). Such generalizations will be given elsewhere.

Acknowledgements: This work has been performed under Austrian Fonds zur Fï¿½rderung der Wissenschaflischen Forschung, Grant # M688

ï¿½ Vladimir L. Kalashnikov, 05.2003, Wien

[email protected]

Hosted by www.Geocities.ws