Multi-Parameter Laser Modes in Paraxial Optics

©   Christoph Koutschan, Erwin Suazo, and Sergei K. Suslov

(Last modified on 20 September 2015.)

The HolonomicFunctions package can be downloaded from
http://www.risc.jku.at/research/combinat/software/HolonomicFunctions/

<<RISC`HolonomicFunctions`

HolonomicFunctions Package version 1.7.1 (09-Oct-2013)
written by Christoph Koutschan
Copyright 2007-2013, Research Institute for Symbolic Computation (RISC),
Johannes Kepler University, Linz, Austria

2.1.1 Green’s Function and Generalized Fresnel Integrals

(* This is Equation (2.1). *)
schroedinger= I*D[ψ[x,t],t]+ a[t]*D[ψ[x,t],x,x]+(I*x*c[t] - I*g[t])*D[ψ[x,t],x]+(-x^2*b[t] + I*d[t] + x*f[t])*ψ[x,t]

Derive the differential equations for α, ..., κ,  starting from the Schrödinger equation, using the Green’s function. We end up with a Riccati-type system, compare Equations (2.41) - (2.47) with (Subscript[c, 0])=0.

(* Plug in the Green's function G(x,y,t) given in (2.3). *)
(* Divide by the exponential term. *)
Numerator[Together[
With[{exp=Exp[I*(α[t]*x^2+δ[t]*x+κ[t]+β[t]*x*y+ϵ[t]*y+γ[t]*y^2)]},
(schroedinger/.ψ->Function@@{{x,t},1/Sqrt[2Pi*μ[t]]*exp})/exp/Sqrt[2]]]]

(* Coefficient comparison w.r.t. x, y, and the imaginary unit. *)
riccati=DeleteCases[Flatten[CoefficientList[%/.Complex[a_,b_]:>a+i*b,{x,y,i}]],0]

(* This gives the Riccati-type system. *)
riccati=Sort[Solve[#==0,First[Cases[#,Derivative[_][_][_],Infinity]]][[1,1]]&/@riccati];TableForm[riccati/.a1_[t]->a1]

(* For initialization only. *)
riccati={α'[t]->-b[t]-2 c[t] α[t]-4 a[t] α[t]^2,β'[t]->-c[t] β[t]-4 a[t] α[t] β[t],γ'[t]->-a[t] β[t]^2,δ'[t]->f[t]+2 g[t] α[t]-c[t] δ[t]-4 a[t] α[t] δ[t],ϵ'[t]->g[t] β[t]-2 a[t] β[t] δ[t],κ'[t]->g[t] δ[t]-a[t] δ[t]^2,μ'[t]->2 (d[t] μ[t]+2 a[t] α[t] μ[t])};

(* This is the fundamental solution (2.4) - (2.9) of the Riccati-type system. *)
eqns24to29={
(Subscript[α, 0])->Function[t,(4*a[t])^(-1)*D[(Subscript[μ, 0])[t],t]/(Subscript[μ, 0])[t]-d[t]/(2*a[t])],
(Subscript[β, 0])->Function[t,-λ[t]/(Subscript[μ, 0])[t]],
λ->Function[t,Exp[-Integrate[c[ζ]-2*d[ζ],{ζ,0,t}]]],
(Subscript[γ, 0])->Function[t,(2*(Subscript[μ, 1])[0])^(-1)*(Subscript[μ, 1])[t]/(Subscript[μ, 0])[t]+d[0]/(2*a[0])],
(Subscript[δ, 0])->Function[t,λ[t]/(Subscript[μ, 0])[t]*Integrate[((f[ξ]-d[ξ]/a[ξ]*g[ξ])*(Subscript[μ, 0])[ξ]+g[ξ]/(2*a[ξ])*(Subscript[μ, 0])'[ξ])/λ[ξ],{ξ,0,t}]],
(Subscript[ϵ, 0])->Function[t,-2*a[t]*λ[t]*(Subscript[δ, 0])[t]/(Subscript[μ, 0])'[t]+8*Integrate[a[s]*σ[s]*λ[s]*(Subscript[μ, 0])[s]*(Subscript[δ, 0])[s]/(Subscript[μ, 0])'[s]^2,{s,0,t}]
+2*Integrate[a[s]*λ[s]/(Subscript[μ, 0])'[s]*(f[s]-d[s]/a[s]*g[s]),{s,0,t}]],
(Subscript[κ, 0])->Function[t,a[t]*(Subscript[μ, 0])[t]/(Subscript[μ, 0])'[t]*(Subscript[δ, 0])[t]^2-4*Integrate[a[s]*σ[s]/(Subscript[μ, 0])'[s]^2*((Subscript[μ, 0])[s]*(Subscript[δ, 0])[s])^2,{s,0,t}]-2*Integrate[a[s]/(Subscript[μ, 0])'[s]*(Subscript[μ, 0])[s]*(Subscript[δ, 0])[s]*(f[s]-d[s]/a[s]*g[s]),{s,0,t}]]};
TableForm[(#1[t]==#2[t]&@@@eqns24to29)]//TraditionalForm

Now we derive the second-order ODE for μ.

(* Start with the Riccati equation for alpha, and substitute (2.4). *)
Together[#1-#2&@@riccati[[1]]/.(eqns24to29[[1]]/.Subscript[x_,0]->x)]

(* We obtain precisely Equation (2.10). *)
(eqn210=Collect[%/Coefficient[%,μ''[t]],μ[t]|μ'[t],Simplify])/.x_[t]->x

Store some relations for the function μ; the coefficients of its ODE are denoted σ and τ, see Equations (2.10) and (2.11).

eqns210to211={
σ[t_]->Coefficient[eqn210,μ[t]]/4,
τ[t_]->-Coefficient[eqn210,μ'[t]],
(* ODE (16) for mu. *)
Derivative[2][μ:((Subscript[μ, 0])|(Subscript[μ, 1]))][t_]->τ[t]*μ'[t]-4*σ[t]*μ[t],
(* Use Wronskian of mu0 and mu1. *)
Derivative[1][(Subscript[μ, 1])][t]->((Subscript[μ, 1])[t] (Subscript[μ, 0])'[t]-2 a[t] (Subscript[μ, 1])[0]*λ[t]^2)/(Subscript[μ, 0])[t]};

(* For initialization only. *)
eqns210to211={σ[t_]->(1/4) (4 a[t] b[t]-4 c[t] d[t]+4 ((d[t])^2)+((2 d[t] a'[t])/(a[t]))-2 d'[t]),τ[t_]->-2 c[t]+4 d[t]+((a'[t])/(a[t])),((μ:(Subscript[μ, 0])|(Subscript[μ, 1])))''[t_]->-4 μ[t] σ[t]+τ[t] μ'[t],(Subscript[μ, 1])'[t]->((-2 a[t] ((λ[t])^2) (Subscript[μ, 1])[0]+(Subscript[μ, 1])[t] (Subscript[μ, 0])'[t])/((Subscript[μ, 0])[t]))}

The initial values for (Subscript[μ, 0]) and (Subscript[μ, 1]):

initsMu={(Subscript[μ, 0])[0]->0,(Subscript[μ, 0])'[0]->2a[0],(Subscript[μ, 1])'[0]->0};

Verify that the fundamental solution (2.4) - (2.9) satisfies the Riccati-type system:

Timing[Together[(#1-#2&@@@riccati)/.x:(α|β|γ|δ|ϵ|κ|μ)->(Subscript[x, 0])//.Join[eqns24to29,eqns210to211]]]

Compute the fundamental solution (2.4) - (2.9) of the Riccati-type system

TableForm[eqns=Equal@@@(Take[riccati,6]/.x:(α|β|γ|δ|ϵ|κ)->(Subscript[x, 0]))]

(* My simplification procedure for integrals. *)
MySimp[expr_]:=FullSimplify[FullSimplify[expr]//.{
(* Rewrite integrals in terms of lambda if possible. *)
Integrate[a1_.*c[x_]+a2_.*d[x_]+a3_.,{x_,a4_,a5_}]/;a2/a1==-2&&FreeQ[a1,x]:>a1*(Log[λ[a4]]-Log[λ[a5]])+Integrate[a3,{x,a4,a5}],
(* Extract constant factors. *)
Integrate[a1_*a2_,{x_,a3__}]/;FreeQ[a1,x]:>a1*Integrate[a2,{x,a3}],
(* Combine integrals with different ranges. *)
Integrate[a1_,{x_,a2_,a3_}]+Integrate[a1_,{x_,a3_,a4_}]:>Integrate[a1,{x,a2,a4}]}/.
K[_]->s]

(* Start with some substitutions. *)
TableForm[
subs=Join[{#1[t]->#2[t]&@@eqns24to29[[1]]},eqns210to211[[{1,2}]],{λ'[t]->(-c[t]+2 d[t]) λ[t],(Subscript[μ, 0])''[t]->τ[t]*(Subscript[μ, 0])'[t]-4*σ[t]*(Subscript[μ, 0])[t]}]/.Pattern->(#&)]

MySimp[DSolve[eqns[[2]]/.subs,(Subscript[β, 0])[t],t]]

(* By choosing C[1] appropriately, we get our solution for beta_0, see Equation (2.5). *)
%[[1,1]]/.C[1]->-λ[1]/(Subscript[μ, 0])[1]
AppendTo[subs,%];

(* Now gamma_0: this doesn't give what we want. *)
MySimp[DSolve[eqns[[3]]/.subs,(Subscript[γ, 0])[t],t]]

(* Trick: use the Wronskian. *)
eqns[[3]]/.subs/.a[t]*λ[t]^2->((Subscript[μ, 1])[t] (Subscript[μ, 0])'[t]-(Subscript[μ, 0])[t]*(Subscript[μ, 1])'[t])/(2  (Subscript[μ, 1])[0])

(* Appropriate choice of C[1] gives the solution gamma_0, see Equation (2.6). *)
DSolve[%,(Subscript[γ, 0])[t],t]/.C[1]->d[0]/(2a[0])

(* Solve the equation for delta_0. *)
MySimp[DSolve[eqns[[4]]/.subs,(Subscript[δ, 0])[t],t]]

(* Choose C[1]. This is the solution for delta_0 given in (2.7). *)
MySimp[%/.C[1]->λ[1]/(2(Subscript[μ, 0])[1])*Integrate[(2 (a[s] f[s]-d[s] g[s]) (Subscript[μ, 0])[s]+g[s] (Subscript[μ, 0])'[s])/(a[s] λ[s]),{s,0,1}]]

(* Add the equation for delta to our list of substitutions. *)
AppendTo[subs,eqns[[4]]/.Equal->Rule];

(* Now epsilon_0: this is not what we want (has mu0 in the denominator). *)
MySimp[DSolve[eqns[[5]]/.subs,(Subscript[ϵ, 0])[t],t]]

(* The solution for epsilon_0. *)
soleps=eqns24to29[[6,2]][t]

(* Check again that it satisfies the equation. *)
Factor[{eqns[[5,2]],D[soleps,t]}//.subs]

(* The solution for epsilon'[t] presented in the paper, but put under a single integral sign. *)
(* This is only the integrand. *)
result=Together[D[eqns24to29[[6,2]][t],t]]

(* Test it. *)
Together[eqns[[5,2]]-result//.subs]

(* The right-hand side of the equation epsilon'[t] = ... *)
eqn5=eqns[[5,2]]/.x_[t]->x

(* The substitutions generate an ideal in some multivariate polynomial ring. *)
ideal=Numerator[Together[#1-#2&@@@subs/.x_[t]->x]]

(* We can prove correctness of the given solution by reduction with Groebner basis. *)
vars=#[[1,0]]&/@subs;
Last[PolynomialReduce[Together[((Subscript[μ, 0])')^2*(eqn5-result/.x_[t]->x)],GroebnerBasis[ideal,vars],vars]]

(* Similarly, we can find a different, but equivalent, expression for epsilon'[t]. *)
(* Play around with vars to influence which variables appear in the result. *)
vars={(Subscript[δ, 0])',(Subscript[β, 0]),b,c};
Factor[Last[PolynomialReduce[(Subscript[μ, 0])'*eqn5,GroebnerBasis[ideal,vars],vars]]/(Subscript[μ, 0])']

(* Together with the desired initial condition, we get this solution for epsilon_0. *)
(Subscript[ϵ, 0])->Function[t,-g[0]/(2a[0])-2Integrate[λ[s]/(Subscript[μ, 0])'[s]*(d[s]+2a[s] (Subscript[α, 0])[s])*(g[s]-2a[s] (Subscript[δ, 0])[s]),{s,0,t}]]

(* Test it. *)
Together[(#1-#2&@@eqns[[5]])/.%//.subs]

(* Another different, still equivalent, expression for epsilon_0'[t], free of alpha_0. *)
vars={(Subscript[α, 0]),(Subscript[β, 0]),(Subscript[δ, 0])',b,c};
Factor[Last[PolynomialReduce[(Subscript[μ, 0])'*eqn5,GroebnerBasis[ideal,vars],vars]]/(Subscript[μ, 0])']

(* In order to get rid of delta', we can do integration by parts. *)
Factor[%-(-2a λ (Subscript[δ, 0])'/(Subscript[μ, 0])')-(Subscript[δ, 0])*(-2(a' λ+a λ')(Subscript[μ, 0])'-(-2a λ)(Subscript[μ, 0])'')/((Subscript[μ, 0])')^2]

(* Perform some of our substitutions and get the integrand of Equation (2.9). *)
Factor[%//.(subs/.x_[t]->x)/.First[Solve[(#1-#2&@@subs[[2]]/.x_[t]->x)==0,d']]]

(* Together with the desired initial condition, we get this solution for epsilon_0. *)
(Subscript[ϵ, 0])->Function[t,-g[0]/(2a[0])+2Integrate[(λ[s]*(a[s]*(f[s]-(Subscript[δ, 0])'[s])-d[s]*g[s])+a[s]*(Subscript[δ, 0])[s]*λ'[s])/(Subscript[μ, 0])'[s],{s,0,t}]]

(* Test it. *)
Together[(#1-#2&@@eqns[[5]])/.%//.subs]

(* Finally, consider kappa_0. *)
MySimp[DSolve[eqns[[6]]/.subs,(Subscript[κ, 0])[t],t]]

(* C[1] is chosen such that we get 0 as lower integration bound. *)
MySimp[%/.C[1]->Integrate[(Subscript[δ, 0])[s] (g[s]-a[s] (Subscript[δ, 0])[s]),{s,0,1}]]

Derive the formula (2.12) for the Wronskian

(* Define the Wronskian and make an ansatz for its differential equation. *)
W[t_]:=Module[{z},Det[Table[D[(Subscript[μ, i])[z],{z,j}],{i,0,1},{j,0,1}]]/.z->t];
ansatz=c0*w[t]+c1*D[w[t],t];
ansatz/.w->W/.Derivative[2][μ_][t]->τ[t]*D[μ[t],t]-4*σ[t]*μ[t]

(* Conditions on the unknown coefficients. *)
DeleteCases[Flatten[CoefficientList[%,Variables[W[t]]]],0]

(* This is the ODE for the Wronskian. *)
deW=ansatz/.c0->-c1*τ[t]/.c1->1/.τ[t]->D[a[t],t]/a[t]-2*c[t]+4*d[t]

(* Its solution. *)
wronsk=ExpandAll[DSolve[deW==0,w[t],t][[1,1,2]]/.K[1]->z]

(* Expand the integral. *)wronsk=wronsk//.HoldPattern[Integrate[a_+b_,{z,a1_,a2_}]]:>Integrate[a,{z,a1,a2}]+Integrate[b,{z,a1,a2}]/.
HoldPattern[Integrate[a1_.*a2_,{z,1,t}]/;FreeQ[a1,z]]:>a1*(Integrate[a2,{z,0,t}]-Integrate[a2,{z,0,1}])

(* Use definition of lambda and simplify. *)
wronsk=ExpandAll[wronsk/.HoldPattern[Integrate[c[z],{z,0,t}]]:>2*Integrate[d[z],{z,0,t}]-Log[λ[t]]]

(* Delete all factors free of t and subsume them into the new constant C[1]. *)
wronsk=C[1]*DeleteCases[wronsk,a_/;FreeQ[a,t]]

Solve[(W[t]-wronsk==0)/.t->0/.initsMu/.λ[0]->1,C[1]]

(* This is the right-hand side of (2.12), used above. *)
wronsk/.First[%]

Asymptotic Expansions (2.13)

Now compute the asymptotics of the fundamental solution. The result corresponds exactly to Equation (2.13):

(* This takes some time since Mathematica tries to simplify the integrals; *)
(* the trick with Hold speeds it up. *)
Timing[
exprs=ReleaseHold[Hold@@{Last[#][t]&/@eqns24to29}//.eqns24to29/.initsMu/.
(* Substitute the Taylor expansion of mu. *)
(μ:(Subscript[μ, 0])|(Subscript[μ, 1]))[t_]:>μ[0]+μ'[0]*t+μ''[0]/2*t^2/.initsMu];
(* Now the series expansions of the above expressions. *)
ser=Series[exprs,{t,0,0}]//.eqns210to211/.initsMu;
]
TableForm[MapThread[Rule,{First[#][t]&/@eqns24to29,ExpandAll[ser]}]]

2.1.2 Special Beam Modes in Weakly Inhomogeneous Media

Generalized Hermite-Gaussian beams and Ermakov-type system

(* This is equation (2.15). *)
HG[x_,t_]:=Exp[I*(α[t]*x^2+δ[t]*x+κ[t]+(2*n+1)*γ[t])]/Sqrt[μ[t]*2^n*n!*Sqrt[Pi]]*Exp[-(β[t]*x+ϵ[t])^2/2]*HermiteH[n,β[t]*x+ϵ[t]];
TraditionalForm[HG[x,t]/.x_[t]->x]

(* These are equations (2.16) - (2.22). *)
eqns216to222={
μ->Function[t,μ[0]*(Subscript[μ, 0])[t]*Sqrt[β[0]^4+4*(α[0]+(Subscript[γ, 0])[t])^2]],
α->Function[t,(Subscript[α, 0])[t]-(Subscript[β, 0])[t]^2*((α[0]+(Subscript[γ, 0])[t])/(β[0]^4+4*(α[0]+(Subscript[γ, 0])[t])^2))],
β->Function[t,-β[0]*(Subscript[β, 0])[t]/Sqrt[β[0]^4+4*(α[0]+(Subscript[γ, 0])[t])^2]],
γ->Function[t,γ[0]-1/2*ArcTan[β[0]^2/(2*(α[0]+(Subscript[γ, 0])[t]))]],
δ->Function[t,(Subscript[δ, 0])[t]-(Subscript[β, 0])[t]*(ϵ[0]*β[0]^3+2*(α[0]+(Subscript[γ, 0])[t])*(δ[0]+(Subscript[ϵ, 0])[t]))/(β[0]^4+4*(α[0]+(Subscript[γ, 0])[t])^2)],
ϵ->Function[t,(2*ϵ[0]*(α[0]+(Subscript[γ, 0])[t])-β[0]*(δ[0]+(Subscript[ϵ, 0])[t]))/Sqrt[β[0]^4+4*(α[0]+(Subscript[γ, 0])[t])^2]],
κ->Function[t,κ[0]+(Subscript[κ, 0])[t]-ϵ[0]*β[0]^3*(δ[0]+(Subscript[ϵ, 0])[t])/(β[0]^4+4*(α[0]+(Subscript[γ, 0])[t])^2)+(α[0]+(Subscript[γ, 0])[t])*((ϵ[0]^2*β[0]^2-(δ[0]+(Subscript[ϵ, 0])[t])^2)/(β[0]^4+4*(α[0]+(Subscript[γ, 0])[t])^2))]};
TableForm[(#1[t]==#2[t]&@@@eqns216to222)/.x_[t]->x]//TraditionalForm

First we verify that (2.15) - (2.22) satisfies the Schrödinger equation (2.1), using the fact that (Subscript[α, 0]),...,(Subscript[κ, 0]) are the fundamental solution of the Riccati-type system.

Timing[Together[
(* Plug (2.15) into (2.1). Use equations (2.16) - (2.22) and the Riccati-type system. *)
schroedinger/.ψ->HG/.eqns216to222/.(riccati/.x:(α|β|γ|δ|ϵ|κ|μ)->(Subscript[x, 0]))/.
(* Rewrite H_{n-2} in terms of H_{n-1} and H_n. *)
HermiteH[n-2,x_]:>(2x*HermiteH[n-1,x]-HermiteH[n,x])/(2n-2)
]]

(* This is the Ermakov-type system, basically Equations (2.41) - (2.47), see below. *)
ermakov=MapThread[#1[[1]]->#1[[2]]+#2&,{riccati,(Subscript[c, 0])*a[t]*β[t]^2*{β[t]^2,0,0,2β[t]*ϵ[t],0,ϵ[t]^2,0}}];
TableForm[ermakov/.x_[t]->x]

Now verify that the expressions (2.16) - (2.22) satisfy the Ermakov-type system with (Subscript[c, 0])=1. Use the fact that (Subscript[α, 0]),...,(Subscript[κ, 0]) are the fundamental solution of the Riccati-type system with (Subscript[c, 0])=0.

Together[(#1-#2&@@@Take[ermakov,6])/.(Subscript[c, 0])->1/.eqns216to222/.(riccati/.x:(α|β|γ|δ|ϵ|κ|μ)->(Subscript[x, 0]))]

Generalized Airy beams and Riccati-type system

(* This is equation (2.23). *)
AB[x_,t_]:=Exp[I*(x^2*α[t]+x*δ[t]+κ[t]-γ[t]*(x*β[t]+ϵ[t]-2/3*γ[t]^2))]/Sqrt[μ[t]]*AiryAi[x*β[t]+ϵ[t]-γ[t]^2];
TraditionalForm[AB[x,t]/.x_[t]->x]

(* These are equations (2.24) - (2.30). *)
eqns224to230={
μ->Function[t,2*μ[0]*(Subscript[μ, 0])[t]*(α[0]+(Subscript[γ, 0])[t])],
α->Function[t,(Subscript[α, 0])[t]-((Subscript[β, 0])[t])^2/(4*(α[0]+(Subscript[γ, 0])[t]))],
β->Function[t,-(β[0]*(Subscript[β, 0])[t])/(2*(α[0]+(Subscript[γ, 0])[t]))],
γ->Function[t,γ[0]-(β[0]^2)/(4*(α[0]+(Subscript[γ, 0])[t]))],
δ->Function[t,(Subscript[δ, 0])[t]-((Subscript[β, 0])[t]*(δ[0]+(Subscript[ϵ, 0])[t]))/(2*(α[0]+(Subscript[γ, 0])[t]))],
ϵ->Function[t,ϵ[0]-(β[0]*(δ[0]+(Subscript[ϵ, 0])[t]))/(2*(α[0]+(Subscript[γ, 0])[t]))],
κ->Function[t,κ[0]+(Subscript[κ, 0])[t]-((δ[0]+(Subscript[ϵ, 0])[t])^2)/(4*(α[0]+(Subscript[γ, 0])[t]))]};
TableForm[(#1[t]==#2[t]&@@@eqns224to230)]//TraditionalForm

Verify that Equations (2.23) - (2.30) are a solution of the Schrödinger equation (2.1).

Timing[Together[
schroedinger/.ψ->AB/.eqns224to230/.(riccati/.x:(α|β|γ|δ|ϵ|κ|μ)->(Subscript[x, 0]))
]]

Plug the solutions (2.24) - (2.30) into the Riccati-type system, i.e. (Subscript[c, 0])=0, and use the fact that (Subscript[α, 0]),...,(Subscript[κ, 0]) are the fundamental solution of the Riccati-type system.

Together[(#1-#2&@@@Take[riccati,6])/.eqns224to230/.(riccati/.x:(α|β|γ|δ|ϵ|κ|μ)->(Subscript[x, 0]))]

2.2.2 Cylindrical Symmetry, Proof of Lemma 1

(* The nonlinear parabolic equation (2.39) in 2D. *)
eqn239=-I*D[A[x,y,s],s]-a[s]*(D[A[x,y,s],x,x]+D[A[x,y,s],y,y])+b[s]*(x^2+y^2)*A[x,y,s]-I*c[s]*(x*D[A[x,y,s],x]+y*D[A[x,y,s],y])-2I*d[s]*A[x,y,s]-(x*(Subscript[f, 1])[s]+y*(Subscript[f, 2])[s])*A[x,y,s]+I*((Subscript[g, 1])[s]*D[A[x,y,s],x]+(Subscript[g, 2])[s]*D[A[x,y,s],y])+h[s]*Abs[A[x,y,s]]^p*A[x,y,s];

(* Substitute (2.38) into Equation (2.39). *)
expr=Together[(eqn239/.
A->Function[{x,y,s},Exp[I*(α[s]*(x^2+y^2)+(Subscript[δ, 1])[s]*x+(Subscript[δ, 2])[s]*y+(Subscript[κ, 1])[s]+(Subscript[κ, 2])[s])]*χ[β[s]*x+(Subscript[ϵ, 1])[s],β[s]*y+(Subscript[ϵ, 2])[s],γ[s]]/μ[s]]/.
h[s]->(Subscript[h, 0])*a[s]*β[s]^2*μ[s]^p)/Exp[I*(α[s]*(x^2+y^2)+(Subscript[δ, 1])[s]*x+(Subscript[δ, 2])[s]*y+(Subscript[κ, 1])[s]+(Subscript[κ, 2])[s])]]

(* Short-hand notation: omit the arguments of χ and the dependencies on s. *)
expr=expr/.χ[__]->χ/.Derivative[a__][χ][__]->Derivative[a][χ]/.a1_[s]->a1

(* Since all functions in question are real-valued, we can simplify the Im[...] part: *)
expr=expr/.Im[_]->0

(* From (2.47) we may conclude that μ[s] is always positive. *)
expr=PowerExpand[Numerator[Together[expr/.Abs[χ/μ]->Abs[χ]/μ]]]

(* Now Equation (2.40), this is the one we want to end up with. *)
eqn240=-I*Derivative[0,0,1][χ]+Derivative[2,0,0][χ]+Derivative[0,2,0][χ]-(Subscript[c, 0])*((β*x+(Subscript[ϵ, 1]))^2+(β*y+(Subscript[ϵ, 2]))^2)*χ-(Subscript[h, 0])*Abs[χ]^p*χ

(* Equate our expression with (2.40) such that |χ|^p disappears. *)
cond=Numerator[Together[eqn240+expr/Coefficient[expr,(Subscript[h, 0])*χ*Abs[χ]^p]]]

(* Coefficient comparison w.r.t. the partial derivatives of χ, plus w.r.t. x and y. *)
cond=Union[DeleteCases[Flatten[CoefficientList[cond/.Derivative[a1_,a2_,a3_][χ]->d1^a1*d2^a2*d3^a3/.χ->1,{d1,d2,d3,x,y}]],0]]

(* Separate real and imaginary parts (all occurring quantities are real). *)
cond=DeleteCases[Flatten[CoefficientList[cond/.Complex[a1_,a2_]:>a1+i*a2,i]],0]

(* Remove trivial factors. The expressions correspond to (2.41)-(2.47). *)
TableForm[cond=If[Head[#]===Times,Select[#,Not[FreeQ[#,Derivative]]&],#]&/@Factor[cond]]

(* Assuming that (Subscript[δ, 1])=(Subscript[δ, 2])=δ etc. we get precisely (2.46). *) 
Factor[cond[[-2]]/.{(Subscript[δ, _])->δ,(Subscript[κ, _])->κ,(Subscript[ϵ, _])->ϵ,(Subscript[g, _])->g}]

3.1 Airy Beams

(* This is Equation (3.1). *)
eqn31=I*D[ψ[x,t],t]+D[ψ[x,t],x,x]

(* Plug in the ansatz (3.2). *)
Together[eqn31/.ψ->Function[{x,t},Exp[I*(x-2t^2/3)*t]*F[x-t^2]]]

(* This reduces to the Airy equation (3.3). *)
Expand[(%/Exp[I*(x-2t^2/3)*t])/.x->z+t^2]

Explicit action of the Schrödinger group (3.5)

(* This is the Schrödinger group (3.5). *)
args={
ξ->ε[0]+((β[0]*x-2*β[0]*δ[0]*t)/(1+4*α[0]*t)),
τ->(β[0]^2*t)/(1+4*α[0]*t)-γ[0]};
SG=Function@@{{x,t},Sqrt[β[0]/(1+4α[0]*t)]*
Exp[I*((α[0]*x^2+δ[0]*x-δ[0]^2*t)/(1+4*α[0]*t)+κ[0])]*χ[ξ,τ]/.args};
TraditionalForm[SG[x,t]]

(* Plug (3.5) into (3.1) and divide out some content. *)
Together[(eqn31/.ψ->SG)/(SG[x,t]/.χ[__]->1)]

(* Write this more nicely. We get the second equation of (3.4). *)
%*((1+4 t*α[0])/β[0])^2/.(Reverse/@args)

Multi-parameter Airy beams (3.6)

MPAB[x_,t_]:=(1+4*α*t)^(-1/2)*Exp[I*((α*x^2+δ*x-δ^2*t)/(1+4*α*t))]*Exp[I*(ε+((β*x-2*β*δ*t)/(1+4*α*t))-(2*(β^4*t^2)/(3*(1+4*α*t)^2)))*(((β^2*t)/(1+4*α*t)))]*AiryAi[ε+((β*x-2*β*δ*t)/(1+4*α*t))-((β^4*t^2)/((1+4*α*t)^2))];
TraditionalForm[MPAB[x,t]]

Timing[FullSimplify[eqn31/.ψ->MPAB]]

Now we derive the Schrödinger equation (3.1) from the definition (3.6), using the HolonomicFunctions package.
At the same time, this proves that (3.6) solves (3.1).

(* This computes a system of PDEs w.r.t. t and x. *)
Timing[ann=Annihilator[MPAB[x,t],{Der[t],Der[x]}]]

(* This finds a linear combination of the above equations such that only B_t and B_xx occur. *)
FindRelation[ann,Support->{Der[t],Der[x]^2}]

Graphical example of “self-accelerating” Airy mode in Equation (3.6)

(* Special case of the Airy mode in Equation (3.6) when α=-1/4, β=1, δ=ε=0. *)
PowerExpand[ComplexExpand[Abs[MPAB[x,t]/.{t->s,α->-1/4,β->1,δ->0,ε->0}]]]

(* Figure 1. Levels of intensity of the Airy mode under consideration. *)ContourPlot[%^2,{x,-10,5},{s,-3,.5},PlotPoints->50,Axes->True,AxesLabel->{x,s}]

Quasi-diffraction-free finite energy Airy beams

(* Let us consider another special case of the Airy beam in Equation (3.6). *)
FEAB0=AiryAi[z]/Sqrt[γ-t]*Exp[w]/.{z->g^(1/3)*(-g*(γ*t)^2/(γ-t)^2+(γ*(x-v*t)/(γ-t))),w->I*((x-v*t/2)*(v*γ))/(2*(γ-t))+I*g*(x-v*t)*(γ^2*t)/(γ-t)^2}

FEAB=FEAB0*Exp[((I (3 (t^2) (x^2)-6 t (x^2) γ+3 (x^2) (γ^2)+8 (g^2) (t^3) (γ^3)))/(12 (((t-γ))^3)))]

(* Let us check that it satisfies the equation. *)
FullSimplify[I*D[FEAB,t]+D[FEAB,x,x]]

(* Compute the equation, from the expression, using HolonomicFunctions. *)
FindRelation[Annihilator[FEAB,{Der[t],Der[x]}],Support->{Der[t],Der[x]^2}]

(* Let us choose the following values of the parameters. *)
FEAB=FEAB/.{v->-I/10,g->1,γ->1};

(* Plotting at t=0. *)
Plot[Abs[FEAB/.t->0]^2,{x,-15,3},AxesLabel->{x,Abs[ψ]^2}]

(* Plotting at t=0.5. *)
Plot[Abs[FEAB/.t->1/2]^2,{x,-5,2},PlotStyle->Blue,AxesLabel->{x,Abs[ψ]^2}]

(* Showing plots at t=0 and t=0.5 together. *)
Show[%,%%,PlotRange->{{-5,2},{0,1/2}}]

(* "Rogue wave" animation *)
Animate[
Plot[Abs[FEAB/.t->tt]^2,{x,-10,6},PlotRange->{0,0.57},AxesLabel->{x,Abs[ψ]^2}],{tt,0,.95}]

3.2 Oscillating and Breathing Hermite-Gaussian Beams

(* This is Equation (3.7). *)
eqn37=2*I*D[A[x,t],t]+D[A[x,t],x,x]-(x^2)*A[x,t]

Standard harmonic oscillator solutions

SHO[x_,t_]:=(((E^(-(x^2)/2)) HermiteH[n,x])/(((2^n) n! (Sqrt[π]))^(1/2)))*Exp[-I*(n+1/2)*t];
TraditionalForm[SHO[x,t]]

TraditionalForm[Table[SHO[x,t],{n,0,2}]]

Timing[FullSimplify[eqn37/.A->SHO]]

“Missing” solutions (3.8) - (3.14)

MS[x_,t_]:=((HermiteH[n,β[t]*x+ε[t]])/((μ[t]*(2^n) n! (Sqrt[π]))^(1/2)))*Exp[-(β[t]*x+ε[t])^2/2]*Exp[I*(α[t]*x^2+δ[t]*x+κ[t])]*Exp[I*(2*n+1)*γ[t]];
TraditionalForm[MS[x,t]]

(* These are Equations (3.9) - (3.14). *)
eqns39to314={μ[t]->μ*Sqrt[β^4*Sin[t]^2+(2*α*Sin[t]+Cos[t])^2],
α[t]->(α*Cos[2*t]+1/4*Sin[2*t]*(β^4+4*α^2-1))/(β^4*Sin[t]^2+(2*α*Sin[t]+Cos[t])^2),
β[t]->β/Sqrt[(β^4*Sin[t]^2+(2*α*Sin[t]+Cos[t])^2)],
γ[t]->γ-(1/2)*ArcTan[β^2*Sin[t]/(2*α*Sin[t]+Cos[t])],
δ[t]->(δ*(2*α*Sin[t]+Cos[t])+ε*β^3*Sin[t])/(β^4*(Sin[t]^2)+(2*α*Sin[t]+Cos[t])^2),
ε[t]->(ε*(2*α*Sin[t]+Cos[t])-β*δ*Sin[t])/Sqrt[β^4*Sin[t]^2+(2*α*Sin[t]+Cos[t])^2],
κ[t]->κ+Sin[t]^2*(ε*β^2*(α*ε-β*δ)-α*δ^2)/(β^4*(Sin[t])^2+(2*α*Sin[t]+Cos[t])^2)+(1/4)*Sin[2*t]*(ε^2*β^2-δ^2)/(β^4*Sin[t]^2+(2*α*Sin[t]+Cos[t])^2)};
TraditionalForm[TableForm[eqns39to314]]

Timing[FullSimplify[eqn37/.A->Function@@{{x,t},MS[x,t]/.eqns39to314}]]

(* With a little effort one can reduce the simplification time drastically. *)
Timing[
(* The exponential part. *)
partExp=Exp[-(β[t]*x+ε[t])^2/2]*Exp[I*(α[t]*x^2+δ[t]*x+κ[t])]*Exp[I*((2*n+1)*γ[t])];
(* The HermiteH part (we replace HermiteH by H in order to avoid automatic evaluation). *)
partH=H[n,β[t]*x+ε[t]];
(* Rewrite Sin[2t] and Cos[2t]. *)
{partExp,partH}={partExp,partH}/.eqns39to314/.{Sin[2t]->2Sin[t]Cos[t],Cos[2t]->Cos[t]^2-Sin[t]^2};
(* Plug the expression into (3.1). *)
test=eqn37/.A->Function@@{{x,t},partExp*partH/Sqrt[μ[t]*2^n*n!*Sqrt[Pi]]/.eqns39to314};
(* Now we have to simplify the 2nd derivative of HermiteH by hand. *)
test=test/.Derivative[0,2][H][n_,x_]:>2*x*Derivative[0,1][H][n,x]-2*n*H[n,x];
(* We divide by the exponential part and clear denominators. *)
test=Numerator[Together[test/partExp]];
(* Consider the coefficients of H and H' separately and simplify. *)
Simplify[Coefficient[test,{partH,Derivative[0,1][H]@@partH}]]
]

Deriving equation (3.7) using the substitution t -> i Log[z]

(* We split (3.8) into two factors. *)
{expr1,expr2}={(μ[t]*2^n*n!*Sqrt[Pi])^(-1/2)*HermiteH[n,β[t]*x+ε[t]]*Exp[-(β[t]*x+ε[t])^2/2]*Exp[I*(α[t]*x^2+δ[t]*x+κ[t])],Exp[I*(2*n+1)*γ[t]]};
MS[x,t]===expr1*expr2

(* Now the substitutions; the last one (t -> I*Log[z]) has the advantage that Sin[t] and Cos[t] turn into rational functions in z. *)
{expr1,expr2}={expr1,expr2}/.eqns39to314/.
{Sin[2t]->2Sin[t]Cos[t],Cos[2t]->Cos[t]^2-Sin[t]^2}/.
t->I*Log[z]/.(a:(Sin|Cos)[_]:>FunctionExpand[a]);
expr2

(* Remove constant factor that confuses Annihilator here; it doesn't change the result. *)
expr2=expr2/.γ->0

Timing[ann=Annihilator[expr1*expr2,{Der[x],Der[z]}]]

(* This is indeed very large. *)
ByteCount/@ann

(* Find an operator with the desired support in the computed left ideal. *)
Timing[rel=FindRelation[ann,Support->{Der[x]^2,Der[z],1}]]

(* This is equation (3.7). *)
ApplyOreOperator[First[rel],A[x,I*Log[z]]]/.I*Log[z]->s

3.3 Hermite-Gaussian Beams

eqn315=2*I*D[B[x,s],s]+D[B[x,s],x,x]

Deriving equation (3.7) from (3.15) and (3.16)

(* Plug the substitution (3.16) into equation (3.15): *)
transEqn=eqn315/.B->Function[{x,s},Exp[I*s*x^2/2/(1+s^2)]/(1+s^2)^(1/4)*A[x/Sqrt[1+s^2],ArcTan[s]]];

(* Divide by the exponential factor and clear denominators. *)
transEqn=Numerator[Together[transEqn/Exp[I*s*x^2/2/(1+s^2)]]]

(* Now substitute the variables x and s in order to get (3.7). *)
Together[(transEqn/.{x->x*Sqrt[1+s^2],ArcTan[s]->s})/(1+s^2)]

Simplify[%==eqn37/.t->s]

Hermite-Gaussian modes (3.17) for the homogeneous parabolic equation

(* If gamma is set to  0, then this is Equation (3.17). *)
HGP[x_,s_]:=(1/((π^(1/4)) (Sqrt[(2^n) n!])))(E^(((I (-s (δ^2)+2 x (x α+δ))-(((x β+ε))^2)-2 s ε (-β δ+α ε))/(2 (1+2 s α+I s (β^2))))+I (1+2 n) (γ-(1/2) ArcTan[((s (β^2))/(1+2 s α))]))) (Sqrt[(β/(Sqrt[(((1+2 s α))^2)+(s^2) (β^4)]))]) HermiteH[n,((x β-s β δ+ε+2 s α ε)/(Sqrt[(((1+2 s α))^2)+(s^2) (β^4)]))];
TraditionalForm[HGP[x,s]]

Timing[FullSimplify[eqn315/.B->HGP]]

(* This computes a system of PDEs for chi w.r.t. s and x. *)
Timing[ann=Annihilator[HGP[x,s],{Der[s],Der[x]}]]

(* This finds a linear combination of the above equations such that only chi_s and chi_xx occur. *)
FindRelation[ann,Support->{Der[s],Der[x]^2}]

Graphical example of 1D self-focusing Gaussian mode in Equation (3.17)

(* This is intensity, or the modulus squared of the Gaussian mode for n=0 in Equation (3.17). *)
Clear[ABS2HGP0];
SetDelayed@@{ABS2HGP0[x_,s_,α_,β_,δ_,ε_],
Simplify[ComplexExpand[Abs[HGP[x,s]]^2/.{n->0,γ->0}]]};
TraditionalForm[ABS2HGP0[x,s,α,β,δ,ε]]

(* Compute the location of the maximum. *)
{max}=Solve[Thread[(D[ABS2HGP0[x,s,α,β,δ,ε],#]&/@{x,s})==0]/.Exp[a_]->1,{x,s}]

(* The maximum is given by *)
Simplify[ABS2HGP0[x,s,α,β,δ,ε]/.max,Element[β,Reals]]

(* Choosing special values of the parameters as follows: α=-2, β=1, δ=ε=0. *)
ABS2HGP0[x,s,-2,1,0,0]

max/.{α->-2, β->1, δ->0,ε->0}

(* The value of the maximum for the given parameters. *)
%%/.%

N[%]

(* Graph of intensity for the Gaussian mode under consideration when α=-2, β=1, δ=ε=0. *)Plot3D[ABS2HGP0[x,s,-2,1,0,0],{x,-1,1},{s,-.105,.575},PlotPoints->20,PlotRange->{0,2.33},AxesLabel->{x,s,"Intensity"}]

(* The corresponding intensity levels. *)ContourPlot[ABS2HGP0[x,s,-2,1,0,0],{x,-1,1},{s,-.105,.575},Axes->True,AxesLabel->{x,s}]

Elegant Hermite-Gaussian modes (3.19)

We derive Equation (3.19) from (3.17) by means of the substitution given in (3.18).

(* Perform the substitution of (3.18) into (3.17). *)
expr=HGP[x,s]/.γ->0/.
((1+2 α s)^2+β^4 s^2)^e_->((1+2α s+I*β^2 s)/Exp[I*ArcTan[β^2 s/(1+2α s)]])^(2e)/.α->I*β^2/2

(* Express the arctan by logarithms for further simplifications. *)
expr=Simplify[expr/.ArcTan[x_]->I/2*(Log[1-I*x]-Log[1+I*x])]

(* Now simplify the square-roots. *)
(* CAVEAT: PowerExpand simplifies Sqrt[x^2] to x. *)
expr=PowerExpand[expr/.Power[x_,a_]:>Power[Factor[x],a]]

(* Enforce that all Sqrt's have argument 1+2isβ^2. *)
expr=Simplify[expr/.Power[x_/;FreeQ[Together[x/(1+2*I*β^2*s)],s],a_]:>Power[Together[x/(1+2*I*β^2*s)],a]*Power[1+2*I*β^2*s,a]]

(* This is Equation (3.19) as given in the paper. *)
EH[x_,s_]:=Sqrt[β/(2^n*n!(1+2 I s β^2)^(n+1)Sqrt[Pi])]*HermiteH[n,((β (x-s δ+I s β ε)+ε)/(Sqrt[1+2 I s (β^2)]))]*Exp[-((2 (β^2)(x^2)-I δ(2x-s δ)+2 β ε(x-s δ)+(1+I s (β^2)) (ε^2))/(2(1+2I s (β^2))))];
TraditionalForm[EH[x,s]]

(* Verify that the above derivation agrees with (3.19). *)
Simplify[PowerExpand[expr/EH[x,s]]]

Timing[FullSimplify[eqn315/.B->EH]]

(* This computes a system of PDEs for EH w.r.t. s and x. *)
Timing[ann=Annihilator[EH[x,s],{Der[s],Der[x]}]]

(* This finds a linear combination of the above equations such that only EH_s and EH_xx occur. *)
FindRelation[ann,Support->{Der[s],Der[x]^2}]

Graphical example of change in elegant Gaussian mode direction of propagation

(* This is a renormalized version of Equation (3.20). *)ABS2EG[x_,s_,δ_,ε_]:=(Sqrt[(1/(1+4 (s^2) (β^4)))])(E^(-((2 ((β^2) (((x-s δ))^2)+β (x-s δ) ε+(1+3 (s^2) (β^4)) (ε^2)))/(1+4 (s^2) (β^4)))))/.β->1/Sqrt[2];
TraditionalForm[ABS2EG[x,s,δ,ε]]

ABS2EG[x,s,0,0]

(* Graph of intensity for the given Gaussian mode from Equation (3.20). *)Plot3D[ABS2EG[x,s,0,0],{x,-2,2},{s,-2.3,2.3},PlotPoints->20,PlotRange->{0,1},AxesLabel->{x,s,"Intensity"}]

(* Levels of intensity from the previous graph. *)ContourPlot[ABS2EG[x,s,0,0],{x,-2,2},{s,-2.3,2.3},Axes->True,AxesLabel->{x,s}]

(* Another choice of parameters. *)
FullSimplify[ABS2EG[x,s,1/2,1/2]]

(* Example of a slanted beam propagation. *)Plot3D[ABS2EG[x,s,1/2,1/2],{x,-4,3.5},{s,-2.1,2.1},PlotPoints->50,PlotRange->{0,1},AxesLabel->{x,s,"Intensity"}]

(* Figure 2. Levels of intensity from the previous graph. *)ContourPlot[ABS2EG[x,s,1/2,1/2],{x,-4,3.5},{s,-2.1,2.1},Axes->True,AxesLabel->{x,s}]

Graphical example of self-focusing

(* This is a verification of the expansion transformation of the Schrödinger group. *)
args={
X->(x/(1+m*s)),
S->s/(1+m*s)};
ET=Function@@{{x,s},Sqrt[1/(1+m*s)]*
Exp[I*((m*x^2)/2/(1+m*s))]*C[X,S]/.args};
TraditionalForm[ET[x,s]]

(* Plug the expansion into (3.15) and divide out some content. *)
Together[(eqn315/.B->ET)/(ET[x,s]/.C[__]->1)]

(* Write this more nicely. We get the same equation (3.15). *)
%*(1+m*s)^2/.(Reverse/@args)/.{X->x,S->s}

(* This is the transformed elegant Gaussian package from the paper up to a constant. *)
EGT[x_,s_]:=Exp[-((2 (x^2) (β^2)+2 (1+m s) x (-I δ+β ε)+(1+m s) ((ε^2)+s (I (δ^2)-2 β δ ε+(m+I (β^2)) (ε^2))))/(2 (1+m s) (1+m s+2 I s (β^2))))] (Sqrt[(β/(1+m s+2 I s (β^2)))])*Exp[I*((m*x^2)/2/(1+m*s))];
TraditionalForm[EGT[x,s]]

Timing[FullSimplify[eqn315/.B->EGT]]

(* This is intensity, or the modulus squared of the Gaussian package under consideration. *)
Clear[ABS2EGT];
SetDelayed@@{ABS2EGT[x_,s_,δ_,ε_,m_],ComplexExpand[Sqrt[2]*Abs[EGT[x,s]]^2]/.β->1/Sqrt[2]/.Exp[a_]:>Exp[Collect[#,ε,Factor]&/@Together[a]]};
TraditionalForm[ABS2EGT[x,s,δ,ε,m]]

{ABS2EGT[x,s,0,0,-2],ABS2EGT[0,1/2,0,0,-2]}

(* Graph of intensity for the given Gaussian mode. *)Plot3D[ABS2EGT[x,s,0,0,-2],{x,-2,2},{s,-1,1.5},PlotPoints->20,PlotRange->{0,2.25},AxesLabel->{x,s,"Intensity"}]

(* Levels of intensity from the previous graph. *)ContourPlot[ABS2EGT[x,s,0,0,-2],{x,-2,2},{s,-1,1.5},Axes->True,AxesLabel->{x,s}]

3.4 Breathing Spiral Laser Beams (Laguerre-Gaussian Modes)

Isotropic planar harmonic oscillator in a perpendicular magnetic field

(* This is Equation (3.22). *)
eqn322=I*D[Ψ[X,Y,T],T]+D[Ψ[X,Y,T],X,X]+D[Ψ[X,Y,T],Y,Y]-(X^2+Y^2)*Ψ[X,Y,T]+I*ω*(Y*D[Ψ[X,Y,T],X]-X*D[Ψ[X,Y,T],Y])

(* This is Equation (3.23) in Cartesian coordinates. *)Psi[X_,Y_,T_]:=Exp[-(X^2+Y^2)/2]* LaguerreL[n,m,X^2+Y^2]*Exp[-I*(4*n+2*(m+1)-m*ω)*T]*(X+I*Y)^m;
TraditionalForm[Psi[X,Y,T]]

TableForm[Table[TraditionalForm[Psi[X,Y,T]],{n,0,2}]]

Timing[FullSimplify[eqn322/.Ψ->Psi]]

Isotropic planar harmonic oscillator

FullSimplify[
Psi[X,Y,T] /.{X->Cos[ω*t]*x+Sin[ω*t]*y,Y->-Sin[ω*t]*x+Cos[ω*t]*y,T->t},Element[m,Integers]]

(* This function does not any more depend on omega. *)
Psi1[x_,y_,t_]:=Evaluate[%];
TraditionalForm[Psi1[x,y,t]]

Timing[FullSimplify[eqn322/.ω->0/.Ψ->Psi1]]

Planar oscillator

(* This is Equation (3.24). *)
eqn324=2*I*D[A[x,y,s],s]+D[A[x,y,s],x,x]+D[A[x,y,s],y,y]-(x^2+y^2)*A[x,y,s]

Psi2[x_,y_,s_]:=(E^((1/2) I (-2 (1+m+2 n) s+I ((x^2)+(y^2))))) (((x+I y))^m) LaguerreL[n,m,(x^2)+(y^2)];
TraditionalForm[Psi2[x,y,s]]

Timing[FullSimplify[eqn324/.A->Psi2]]

(* Compute a (Groebner) basis for the ideal of operators that annihilate Psi2. *)
ann=Annihilator[Psi2[x,y,s],{Der[s],Der[x],Der[y]}]

(* This proves that Psi2 satisfies equation (3.24). *)
OreReduce[2*I*Der[s]+Der[x]^2+Der[y]^2-x^2-y^2,ann]

(* This tries to FIND equation (3.24) automatically. *)
FindRelation[ann,Support->{1,Der[s],Der[x]^2,Der[y]^2}]

(* Equation (3.24) corresponds to a linear combination of the above two operators. *)
(2I)**%[[1]]+%[[2]]

(* Here is another ad hoc trick to FIND (3.24) by observing that it is free of m and n. *)
FindRelation[ann,Eliminate->{m,n}]

Simplify[ApplyOreOperator[First[%],A[x,y,s]]==eqn324]

3.5 Laguerre-Gaussian Beams

(* This is Equation (3.26). *)
eqn326=2*I*D[B[x,y,s],s]+D[B[x,y,s],x,x]+D[B[x,y,s],y,y]

Transformation (3.27)

(* Plug the substitution (3.27) into equation (3.26): *)
transEqn=eqn326/.B->Function[{x,y,s},Exp[I*s*(x^2+y^2)/2/(1+s^2)]/(1+s^2)^(1/2)*A[x/Sqrt[1+s^2],y/Sqrt[1+s^2],ArcTan[s]]];

(* Divide by the exponential factor and clear denominators. *)
transEqn=Numerator[Together[transEqn/Exp[I*s*(x^2+y^2)/2/(1+s^2)]]]

(* Now substitute the variables x, y and s in order to get (3.24). *)
Together[(transEqn/.{x->x*Sqrt[1+s^2],y->y*Sqrt[1+s^2],ArcTan[s]->s})/(1+s^2)]

Simplify[%==eqn324/.t->s]

Corresponding solution of equation (3.26) obtained by transformation (3.27)

(* This is one of the Laguerre-Gaussian solution of the parabolic equation (3.26). *)
LG[x_,y_,s_]:=(E^(-(((x^2)+(y^2))/(2+2 (s^2)))-I (1+m+2 n) ArcTan[s])) (((((x+I y)/(Sqrt[1+(s^2)]))))^m) LaguerreL[n,m,(((x^2)+(y^2))/(1+(s^2)))]*(1+s^2)^(-1/2)*Exp[I*s*(x^2+y^2)/(2*(1+s^2))];
TraditionalForm[LG[x,y,s]]

Timing[FullSimplify[eqn326/.B->LG]]

(* This computes a system of PDEs for LG w.r.t. s, x and y. *)
Timing[ann=Annihilator[LG[x,y,s],{Der[s],Der[x],Der[y]}]]

(* This finds a linear combination of the above operators such that only LG_s, LG_xx, LG_yy occur. *)
(* The result is exactly Equation (3.26). *)
FindRelation[ann,Support->{Der[s],Der[x]^2,Der[y]^2}]

Multi-parameter Laguerre-Gaussian solution (3.28) of equation (3.26)

(* This is a similar verification that the multi-parameter solution (3.28) safisfies the parabolic equation (3.26). *)
MLG[x_,y_,s_]:=(1/(1+2 s α+I*s (β^2)))(E^(-((((x^2)+(y^2)) (-2 I α+(β^2))+2 x (-I δ1+β ε1)+2 y (-I δ2+β ε2)+I s (((δ1)^2)+((δ2)^2))-2 s β (δ1 ε1+δ2 ε2)+(1+2 s α) (((ε1)^2)+((ε2)^2)))/(2 (1+2 s α+I s (β^2))))-I (m+2 n) ArcTan[((s (β^2))/(1+2 s α))]))  (((((β (x+I y-s (δ1+I δ2))+(1+2 s α)(ε1+I ε2))/(Sqrt[(((1+2 s α))^2)+(s^2) (β^4)]))))^m)  LaguerreL[n,m,(((((x β-s β δ1+ε1+2 s α ε1))^2)+(((y β-s β δ2+ε2+2 s α ε2))^2))/((((1+2 s α))^2)+(s^2) (β^4)))];
TraditionalForm[MLG[x,y,s]]

Timing[FullSimplify[eqn326/.B->MLG]]

(* This is a derivation of the parabolic equation (3.26) for the multi-parameter solution (3.28). *)
Timing[ann=Annihilator[MLG[x,y,s],{Der[s],Der[x],Der[y]}]]

ByteCount[ann]

(* As before, we find an appropriate linear combination of the above operators. *)
(* Again, the result is Equation (3.26). *)
Timing[FindRelation[ann,Support->{Der[s],Der[x]^2,Der[y]^2}]]

Graphical example of 2D self-focusing Gaussian mode in Equation (3.28)

(* This is intensity, or the modulus squared of the Gaussian mode n=m=0 in Equation (3.28). *)ABS2MG[x_,y_,s_,α_,β_,δ1_,δ2_,ε1_,ε2_]:=(1/((((1+2 s α))^2)+(s^2) (β^4)))*Exp[-((x^2) (β^2)+(y^2) (β^2)+(s^2) (β^2) ((δ1)^2)-2 s y (β^2) δ2+(s^2) (β^2) ((δ2)^2)-2 s β δ1 ε1-4 (s^2) α β δ1 ε1+((ε1)^2)+4 s α ((ε1)^2)+4 (s^2) (α^2) ((ε1)^2)+2 x β (-s β δ1+ε1+2 s α ε1)+2 (1+2 s α) β (y-s δ2) ε2+(((ε2+2 s α ε2))^2))/((((1+2 s α))^2)+(s^2) (β^4))]
TraditionalForm[ABS2MG[x,y,s,α,β,δ1,δ2,ε1,ε2]]

(* Checking the location of the maximum. *){D[ABS2MG[x,y,s,α,β,δ1,δ2,ε1,ε2],x],D[ABS2MG[x,y,s,α,β,δ1,δ2,ε1,ε2],y],D[ABS2MG[x,y,s,α,β,δ1,δ2,ε1,ε2],s]};
Simplify[%];
%/.{x->-((2 α δ1+(β^3) ε1)/(4 (α^2)+(β^4))),y->-((2 α δ2+(β^3) ε2)/(4 (α^2)+(β^4))),s->-((2 α)/(4 (α^2)+(β^4)))};
Simplify[%]

(* The maximum is given by *)
ABS2MG[x,y,s,α,β,δ1,δ2,ε1,ε2]/.{x->-((2 α δ1+(β^3) ε1)/(4 (α^2)+(β^4))),y->-((2 α δ2+(β^3) ε2)/(4 (α^2)+(β^4))),s->-((2 α)/(4 (α^2)+(β^4)))};
Simplify[%]

(* Choosing special values of the parameters as follows α=-2, β=1, δ1=δ2=ε1=ε2=0. *){ABS2MG[x,y,s,-2,1,0,0,0,0],ABS2MG[0,0,0,-2,1,0,0,0,0],ABS2MG[0,0,4/17,-2,1,0,0,0,0]}

(* Level surface where intensity drops by e. *)
ContourPlot3D[ABS2MG[x,y,s,-2,1,0,0,0,0]-(1/E)==0,{x,-1,1},{y,-1,1},{s,-.25,.75},PlotPoints->5,AxesLabel->{x,y,s}]

Multi-parameter “elegant” Laguerre-Gaussian modes (3.29)

ELG[x_,y_,s_]:=
(E^(-((2 ((x^2)+(y^2)) (β^2)+I s (((δ1)^2)+((δ2)^2))+2 x (-I δ1+β ε1)+2 y (-I δ2+β ε2)-2 s β (δ1 ε1+δ2 ε2)+(1+I s (β^2)) (((ε1)^2)+((ε2)^2)))/(2 (1+2 I s (β^2)))))) (((1+2 I s (β^2)))^(-1-m-n)) (((β (x+I y-s (δ1+I δ2))+(1+I s (β^2)) (ε1+I ε2)))^m) LaguerreL[n,m,(((((x β-s β δ1+ε1+I s (β^2) ε1))^2)+(((y β-s β δ2+ε2+I s (β^2) ε2))^2))/(1+2 I s (β^2)))];
TraditionalForm[ELG[x,y,s]]

Timing[FullSimplify[eqn326/.B->ELG]]

(* This is a derivation of the parabolic equation (3.26) for the multi-parameter solution (3.29). *)
Timing[ann=Annihilator[ELG[x,y,s],{Der[s],Der[x],Der[y]}];]

ByteCount[ann]

(* This finds a linear combination of the above equations such that only ELG_s, ELG_xx and ELG_yy occur. Once again, the result is Equation (3.26). *)
Timing[FindRelation[ann,Support->{Der[s],Der[x]^2,Der[y]^2}]]

Multi-parameter “diffraction free” Laguerre modes (3.30)

DL[x_,y_,s_]:=
(E^((1/2) (-I s (((δ1)^2)+((δ2)^2))-2 x (-I δ1+β ε1)-2 y (-I δ2+β ε2)+2 s β (δ1 ε1+δ2 ε2)-(1-I s (β^2)) (((ε1)^2)+((ε2)^2))))) (((β (x+I y-s (δ1+I δ2))+(1-I s (β^2)) (ε1+I ε2)))^m) *(((1-2 I s (β^2)))^n) LaguerreL[n,m,(((((x β-s β δ1+ε1-I s (β^2) ε1))^2)+(((y β-s β δ2+ε2-I s (β^2) ε2))^2))/(1-2 I s (β^2)))];
TraditionalForm[DL[x,y,s]]

Timing[FullSimplify[eqn326/.B->DL]]

(* This is a derivation of the parabolic equation (3.26) for the multi-parameter solution (3.30). *)
Timing[ann=Annihilator[DL[x,y,s],{Der[s],Der[x],Der[y]}];]

ByteCount[ann]

(* This finds a linear combination of the above equations such that only DL_s, DL_xx and DL_yy occur. Once again, the result is Equation (3.26). *)
Timing[FindRelation[ann,Support->{Der[s],Der[x]^2,Der[y]^2}]]

3.6 Bessel-Gaussian Beams

Simple Bessel-Gaussian solution of equation (3.26)

BG[x_,y_,s_]:=(((E^((((((1-I s))^2) t)/(1+(s^2)))+((I s ((x^2)+(y^2)))/(2 (1+(s^2))))-(((x^2)+(y^2))/(2+2 (s^2)))-I (1+m) ArcTan[s])) (((((x+I y)/(Sqrt[1+(s^2)]))))^m) Hypergeometric0F1[1+m,-(((((1-I s))^2) t ((x^2)+(y^2)))/(((1+(s^2)))^2))])/(Sqrt[1+(s^2)]));
TraditionalForm[BG[x,y,s]]

Timing[FullSimplify[eqn326/.B->BG]]

Derive equation (3.26) from the simple Bessel-Gaussian solution, using the HolonomicFunctions package. As a by-product, this also proves that BG satisfies equation (3.26).

Timing[ann=Annihilator[BG[x,y,s],{Der[s],Der[x],Der[y]}]]

(* This finds a linear combination of the above equations such that only BG_s, BG_xx, BG_yy occur. *)
FindRelation[ann,Support->{Der[s],Der[x]^2,Der[y]^2}]

Multi-parameter Bessel-Gaussian beams (3.32)

This solution is obtained with the help of (3.31).

MBG[x_,y_,s_]:=(1/(1+2 s α+I*s*(β^2)))(E^(-((((x^2)+(y^2)) (-2 I α+(β^2))+I s (((δ1)^2)+((δ2)^2))+2 x (-I δ1+β ε1)+2 y (-I δ2+β ε2)-2 s β (δ1 ε1+δ2 ε2)+(1+2 s α) (((ε1)^2)+((ε2)^2)))/(2 (1+2 s α+I s (β^2)))))) (((((β (x+I y-s (δ1+I δ2))+(1+2 s α) (ε1+I ε2))/(1+2 s α+I*s*(β^2)))))^m) Hypergeometric0F1[m+1,-(((((x β-s β δ1+ε1+2 s α ε1))^2)+(((y β-s β δ2+ε2+2 s α ε2))^2))/(((1+2 s α +I*s*(β^2)))^2))*t]*Exp[t*((1+2 s α -I*s*(β^2))/(1+2 s α+I*s* (β^2)))];
TraditionalForm[MBG[x,y,s]]

Timing[FullSimplify[eqn326/.B->MBG]]

Again, derive equation (3.26), but now from the solution (3.32), using the HolonomicFunctions package.

Timing[ann=Annihilator[MBG[x,y,s],{Der[s],Der[x],Der[y]}]]

Timing[FindRelation[ann,Support->{Der[s],Der[x]^2,Der[y]^2}]]

Multi-parameter elegant Bessel-Gaussian beams (3.33)

EBG[x_,y_,s_]:=
(1/(1+2 I s (β^2)))(E^((t/(1+2 I s (β^2)))-((2 ((x^2)+(y^2)) (β^2)+I s (((δ1)^2)+((δ2)^2))+2 x (-I δ1+β ε1)+2 y (-I δ2+β ε2)-2 s β (δ1 ε1+δ2 ε2)+(1+I s (β^2)) (((ε1)^2)+((ε2)^2)))/(2 (1+2 I s (β^2)))))) (((((β (x+I y-s (δ1+I δ2))+(1+I s (β^2)) (ε1+I ε2))/(1+2 I s (β^2)))))^m) Hypergeometric0F1[1+m,-((t ((((x β-s β δ1+ε1+I s (β^2) ε1))^2)+(((y β-s β δ2+ε2+I s (β^2) ε2))^2)))/(((1+2 I s (β^2)))^2))];
TraditionalForm[EBG[x,y,s]]

Timing[FullSimplify[eqn326/.B->EBG]]

Again, derive equation (3.26), but now from the solution (3.33), using the HolonomicFunctions package.

Timing[ann=Annihilator[EBG[x,y,s],{Der[s],Der[x],Der[y]}];]

ByteCount[ann]

Timing[FindRelation[ann,Support->{Der[s],Der[x]^2,Der[y]^2}]]

Multi-parameter “diffraction free” Bessel beams (3.34)

DB[x_,y_,s_]:=
(E^(t (1-2 I s (β^2))+(1/2) (-I s (((δ1)^2)+((δ2)^2))-2 x (-I δ1+β ε1)-2 y (-I δ2+β ε2)+2 s β (δ1 ε1+δ2 ε2)-(1-I s (β^2)) (((ε1)^2)+((ε2)^2))))) (((β (x+I y-s (δ1+I δ2))+(1-I s (β^2)) (ε1+I ε2)))^m) Hypergeometric0F1[1+m,-t ((((x β-s β δ1+ε1-I s (β^2) ε1))^2)+(((y β-s β δ2+ε2-I s (β^2) ε2))^2))];
TraditionalForm[DB[x,y,s]]

Timing[FullSimplify[eqn326/.B->DB]]

Again, derive equation (3.26), but now from the solution (3.34), using the HolonomicFunctions package.

Timing[ann=Annihilator[DB[x,y,s],{Der[s],Der[x],Der[y]}];]

Timing[FindRelation[ann,Support->{Der[s],Der[x]^2,Der[y]^2}]]

Graphical examples of “diffraction free” Bessel beams (3.34)

(* This is intensity, or the modulus squared of the "diffraction-free" Bessel mode *)
(* for m=0 and t=-1 in Equation (3.34). *)
ABS2DB[x_,y_,s_,α_,β_,δ1_,δ2_,ε1_,ε2_]:=
Exp[-2-2 x β ε1-((ε1)^2)-ε2 (2 y β+ε2)+2 s β (δ1 ε1+δ2 ε2)]*Abs[Hypergeometric0F1[1,((((x β-s β δ1+ε1-I s (β^2) ε1))^2)+(((y β-s β δ2+ε2-I s (β^2) ε2))^2))]]^2;
TraditionalForm[ABS2DB[x,y,s,α,β,δ1,δ2,ε1,ε2]]

{ABS2DB[x,y,s,α,β,δ1,δ2,ε1,ε2],ABS2DB[x,y,s,0,1/2,0,0,0,0],ABS2DB[.1,.1,0.1,0,1/2,0,0,0,0]}

(* Sample diffraction-free level surface. *)
ContourPlot3D[ABS2DB[x,y,s,0,1/2,0,0,0,0]==.15,{x,-.51,.51},{y,-.51,.51},{s,-.1,.1},PlotPoints->3,AxesLabel->{x,y,s}]

{ABS2DB[x,y,s,0,1/2,.5,.5,0,0],ABS2DB[.25,.25,0,0,1/2,.5,.5,0,0]}

(* Sample "slanted" diffraction-free level surface. *)ContourPlot3D[ABS2DB[x,y,s,0,1/2,.5,.5,0,0]==.15,{x,-.51,.51},{y,-.51,.51},{s,-.1,.1},PlotPoints->3,AxesLabel->{x,y,s}]

3.7 Spiral Beams

Transformation (3.35)

(* This is Equation (3.36). *)
eqn336=2*I*D[C[x,y,s],s]+D[C[x,y,s],x,x]+D[C[x,y,s],y,y]-(x^2+y^2)*C[x,y,s]-2*I*ω*(x*D[C[x,y,s],y]-y*D[C[x,y,s],x])

(* Plug the substitution (3.35) into equation (3.26): *)
transEqn=eqn326/.B->Function[{x,y,s},Exp[I*s*(x^2+y^2)/2/(1+s^2)]/(1+s^2)^(1/2)*C[(Cos[ω*ArcTan[s]]*x+Sin[ω*ArcTan[s]]*y)/Sqrt[1+s^2],(-Sin[ω*ArcTan[s]]*x+Cos[ω*ArcTan[s]]*y)/Sqrt[1+s^2],ArcTan[s]]];

(* Divide by the exponential factor and clear denominators. *)
transEqn=FullSimplify[Numerator[Together[transEqn/Exp[I*s*(x^2+y^2)/2/(1+s^2)]]]]

Simplify[transEqn/(1+s^2)/.
{x->Sqrt[1+s^2]*(x*Cos[ArcTan[s]*ω]-y*Sin[ArcTan[s]*ω]),
y->Sqrt[1+s^2]*(x*Sin[ArcTan[s]*ω]+y*Cos[ArcTan[s]*ω])}/.
ArcTan[s]->s]

Simplify[%==eqn336]

Solution (3.37)

(* This is solution (3.37) of Equation (3.26) in the laboratory frame of reference. *)Psi0[x_,y_,s_]:=(((E^((1/2) I ((((x^2)+(y^2))/(-I+s))-2 (1+m+2 n) ArcTan[s]))) (((((x+I y)/(Sqrt[1+(s^2)]))))^m)*LaguerreL[n,m,(((x^2)+(y^2))/(1+(s^2)))])/(Sqrt[1+(s^2)]));
TraditionalForm[Psi0[x,y,s]]

Timing[FullSimplify[eqn326/.B->Psi0]]

Identity (3.42)

(* This is Equation (3.41) in Cartesian coordinates. *)
(* sgn is either +1 or -1. *)
PsiLG0[X_,Y_,T_]:=Exp[-I*k*T]*(X+I sgn Y)^m*Exp[-(X^2+Y^2)/2]* LaguerreL[n,m,X^2+Y^2];
TraditionalForm[PsiLG0[X,Y,T]]

(* This is Equation (3.42). *)
eqn342=eqn336-2*(k+sgn m ω-m-2 n-1)*C[x,y,s]/.{x->X,y->Y,s->T}

(* Plug (3.41) into Equation (3.42). *)
Timing[FullSimplify[eqn342/.C->PsiLG0/.{{sgn->-1},{sgn->1}}]]

Now derive Equation (3.42) with the HolonomicFunctions package. We do both cases (sgn=1 and sgn=-1) separately.

ann=Annihilator[PsiLG0[X,Y,T]/.sgn->1,{Der[X],Der[Y],Der[T]}]

(* Equation (3.42) in operator notation. *)
ToOrePolynomial[eqn342,C[X,Y,T]]/.sgn->1

(* Reduction shows that the above operator is in the ideal ann. *)
OreReduce[%,ann/.sgn->1]

(* Find Equation (3.42) in the annihilator ideal. *)
First[FindRelation[ann,Support->{1,Der[X]^2,Der[Y]^2}]]+2I*(ann[[1]]+ω ann[[2]])

(* Compare this derived operator with the equation given in the paper. *)
Simplify[ApplyOreOperator[%,C[X,Y,T]]==eqn342/.sgn->1]

ann=Annihilator[PsiLG0[X,Y,T]/.sgn->-1,{Der[X],Der[Y],Der[T]}]

(* Equation (3.42) in operator notation. *)
ToOrePolynomial[eqn342,C[X,Y,T]]/.sgn->-1

(* Reduction shows that the above operator is in the ideal ann. *)
OreReduce[%,ann/.sgn->-1]

(* Find Equation (3.42) in the annihilator ideal. *)
First[FindRelation[ann,Support->{1,Der[X]^2,Der[Y]^2}]]+2I*(ann[[1]]+ω ann[[2]])

(* Compare this derived operator with the equation given in the paper. *)
Simplify[ApplyOreOperator[%,C[X,Y,T]]==eqn342/.sgn->-1]

3.8 “Smart” Lens Design

Gaussian beam and its energy density inside the lens

(* This is the Gaussian package from Equation (3.8) for n=0. *)
Clear[MS0];
SetDelayed@@{MS0[x_,t_],Simplify[MS[x,t]/.eqns39to314/.{n->0,μ->1/β,γ->0,κ->0}]};
TraditionalForm[MS0[x,s]]

(* We verify Equation (3.7). *)
Timing[FullSimplify[eqn37/.A->MS0]]

(* Intensity of the beam inside the lens. *)
Clear[ABS2MS0];
SetDelayed@@{ABS2MS0[x_,s_,α_,β_,δ_,ε_],FullSimplify[ComplexExpand[Abs[MS0[x,s]]^2],Element[{x,s,α,β,δ,ε},Reals]]};
TraditionalForm[ABS2MS0[x,s,α,β,δ,ε]]

ABS2MS0[0,0,α,β,δ,ε]

(* Looking for value of parameters: plot inside a long lens when α=2, β=1, and 0<s<Pi=l. *)Plot3D[ABS2MS0[x,s,2,1,0,0],{x,-1.2,1.2},{s,0,Pi},PlotPoints->40,PlotRange->{0,2.4},AxesLabel->{x,s,"Intensity"}]

(* Plot of intensity inside of a shorter lens when α=2, β=1, and 0<s<2.75=l. *)Plot3D[ABS2MS0[x,s,2,1,0,0],{x,-1.2,1.2},{s,0,2.75},PlotRange->{0,2.4},AxesLabel->{x,s,"Intensity"}]

Gaussian beam and its energy density outside the lens

(* This is the Gaussian package for n=0 in Equation (3.17). *)
Clear[HGP0];
SetDelayed@@{HGP0[x_,s_],PowerExpand[HGP[x,s]/.{n->0,γ->0}]};
TraditionalForm[HGP0[x,s]]

(* We verify Equation (3.15). *)
Timing[FullSimplify[eqn315/.B->HGP0]]

Clear[ABS2HGP0];
SetDelayed@@{ABS2HGP0[x_,s_,α_,β_,δ_,ε_],
Simplify[ComplexExpand[FullSimplify[Abs[HGP0[x,s]]^2,Element[{x,s,α,β,δ,ε},Reals]]]]};
TraditionalForm[ABS2HGP0[x,s,α,β,δ,ε]]

ABS2HGP0[0,0,α,β,δ,ε]

(* Plot of intensity outside the lens when α=2, β=1, and s<0. *)Plot3D[ABS2HGP0[x,s,2,1,0,0],{x,-1.2,1.2},{s,-.5,0},PlotPoints->20,PlotRange->{0,2.4},AxesLabel->{x,s,"Intensity"}]

(* We verify the continuity at s=0. *)
{HGP0[x,0],MS0[x,0]}
FullSimplify[HGP0[x,0]-MS0[x,0],β>0]

test={ABS2MS0[0,s,2,1,0,0],ABS2HGP0[0,s,α,β,0,0]};
Simplify[Last[Solve[Simplify[Equal@@@{test^2,D[test,s]}/.s->11/4],{α,β}]]]

{α1,β1}=Last/@%;
N[{α1,β1},15]

(* Intensity outside the lens when α=α1, β=β1, s>2.75=l. *)
FullSimplify[ABS2HGP0[x,s,N[α1],N[β1],0,0]]

(* We verify (symbolically!) the continuity of intensity at s=2.75=l. *)FullSimplify[ABS2MS0[x,11/4,2,1,0,0]-ABS2HGP0[x,11/4,α1,β1,0,0]]

(* Plot of intensity outside the lens when s>2.75=l. *)Plot3D[ABS2HGP0[x,s,α1,β1,0,0],{x,-1.2,1.2},{s,2.75,3.25},PlotRange->{0,2.4},AxesLabel->{x,s,"Intensity"}]

Example: Gaussian beam propagation through the lens

(* Figure 3. Propagation of the beam through the lens: Combine the last three graphs together. *)Plot3D[Which[s<0,ABS2HGP0[x,s,2,1,0,0],s<2.75,ABS2MS0[x,s,2,1,0,0],True,ABS2HGP0[x,s,α1,β1,0,0]],{x,-1.2,1.2},{s,-0.5,3.25},PlotPoints->50,PlotRange->{0,2.4},AxesLabel->{x,s,"Intensity"}]

4 Gröbner Basis Methods

(* Polynomial ideal generated by the equations of the Ermakov-type system. *)
idealErm=Most[#1-#2&@@@ermakov]/.x_[t]->x/.(Subscript[c, 0])->1

(* Eliminate a, b, c, ... *)
Timing[GroebnerBasis[idealErm,{α,β,δ,ϵ,α',β',γ',δ',ϵ',κ'},{a,b,c,d,f,g}]]

(* Try to eliminate the derivatives. There is no such relation. *)
Timing[GroebnerBasis[idealErm,{α,β,δ,ϵ,a,b,c,d,f,g},{α',β',γ',δ',ϵ',κ'}]]