Clear["Global`*"]

Non-Adiabatic theories
comparison

Parameters to set:

λ=788*(10^(-9));(* in m *)
Ip0 = 24.59/27.21;(* He: 24.59, Ne: 21.56, H: 13.5984 in au*) ;

αN =1.38; (* He: 1.38, Ne: 2.67, H: 0 *)
αI = 9/32; (* He: 9/32, Ne: 1.33, H: 0 *) 

Nt = 300;

Parameters automatically:

c=2.998*(10^8);
ω = 2π*c/λ;
TU = 2.418*(10^(-17)); (* 1 au = ... s *)
ωAU = ω*TU;

Fieldx[t_,F0_,ϵ_]:=F0*(1/Sqrt[ϵ^2+1]*Cos[ωAU*t]);
Fieldy[t_,F0_,ϵ_]:=F0*(ϵ/Sqrt[ϵ^2+1]*Sin[ωAU*t]);
F[t_,F0_,ϵ_]:=Sqrt[Fieldx[t,F0,ϵ]^2 + Fieldy[t,F0,ϵ]^2];

Ip[t_,F0_, ϵ_]:=Ip0+(αN-αI)*F[t,F0,ϵ]^2/2;
γ[t_,F0_, ϵ_]:=ωAU*Sqrt[2*Ip[0,F0,ϵ]]/F[t,F0,ϵ];

PPT

s0[t_,F0_, ϵ_]:=s/.FindRoot[ArcTanh[(Sqrt[(((s^2)+((γ[t,F0,ϵ])^2))/(1+((γ[t,F0,ϵ])^2)))])]==(ϵ/(ϵ-s))((Sqrt[(((s^2)+((γ[t,F0,ϵ])^2))/(1+((γ[t,F0,ϵ])^2)))])),{s,ϵ/2}]

cx[t_,F0_, ϵ_]:= ((s0[0,F0,ϵ]*(1-(ϵ^2)))/((ϵ-s0[0,F0,ϵ])*(1-ϵ*s0[0,F0,ϵ])))(Sqrt[((((s0[0,F0,ϵ])^2)+((γ[0,F0,ϵ])^2))/(1+((γ[0,F0,ϵ])^2)))]);
cy[t_,F0_, ϵ_]:= (ϵ/(ϵ-s0[0,F0,ϵ]))(Sqrt[((((s0[0,F0,ϵ])^2)+((γ[0,F0,ϵ])^2))/(1+((γ[0,F0,ϵ])^2)))])+(((s0[0,F0,ϵ])^2)/(1-ϵ s0[0,F0,ϵ]))(Sqrt[((1+((γ[0,F0,ϵ])^2))/(((s0[0,F0,ϵ])^2)+((γ[0,F0,ϵ])^2)))]);
cz[t_,F0_, ϵ_]:= (ϵ/(ϵ-s0[0,F0,ϵ]))(Sqrt[((((s0[0,F0,ϵ])^2)+((γ[0,F0,ϵ])^2))/(1+((γ[0,F0,ϵ])^2)))]);

g[t_,F0_, ϵ_]:= (1+((1+(ϵ^2))/(2*((γ[t,F0,ϵ])^2))))*ArcSinh[(Sqrt[((((γ[t,F0,ϵ])^2)+((s0[t,F0,ϵ])^2))/(1-((s0[t,F0,ϵ])^2)))])]-((1+(ϵ^2)-2*ϵ*s0[t,F0,ϵ])/(2*((γ[t,F0,ϵ])^2)*(1-((s0[t,F0,ϵ])^2))))(Sqrt[(1+((γ[t,F0,ϵ])^2))(((s0[t,F0,ϵ])^2)+((γ[t,F0,ϵ])^2))])

p0[F0_, ϵ_]:=((Sqrt[2*Ip[0,F0,ϵ]])/(γ[0,F0, ϵ]))*(ϵ-s0[0,F0,ϵ])*(Sqrt[((1+((γ[0,F0,ϵ])^2))/(1-((s0[0,F0,ϵ])^2)))]);

PPT Ionisation Probability

PPTProbab [t_,F0_,ϵ_, px_, py_, pz_]:=Exp[-((2*Ip[0,F0,ϵ])/(ωAU))*g[t,F0,ϵ]]*Exp[-((cx[t,F0,ϵ]*(px^2))/(ωAU))-((cz[t,F0,ϵ]*(pz^2))/(ωAU))]*(Exp[-((cy[t,F0,ϵ]*(((py-p0[F0,ϵ]))^2))/(ωAU))]+Exp[-((cy[t,F0,ϵ]*(((py+p0[F0,ϵ]))^2))/(ωAU))]);
PPTProbabNorm[t_,F0_,ϵ_, px_, py_, pz_] := PPTProbab[t,F0,ϵ, px, py, pz]/PPTProbab[0,F0,ϵ, 0, p0[F0,ϵ],0];

Bondar: linear polarisation only!

α[t_,F0_,px_,po_]:= ((γ[t,F0,0])/2)*(√((((px^2)+(po^2))/(2*Ip[0,F0,0]))+((2*px)/(γ[t,F0,0]*(Sqrt[2*Ip[0,F0,0]])))+(1/((γ[t,F0,0])^2))+1)+√((((px^2)+(po^2))/(2*Ip[0,F0,0]))-((2*px)/(γ[t,F0,0]*(Sqrt[2*Ip[0,F0,0]])))+(1/((γ[t,F0,0])^2))+1));
β[t_,F0_,px_,po_]:= ((γ[t,F0,0])/2)*(√((((px^2)+(po^2))/(2*Ip[0,F0,0]))+((2*px)/(γ[t,F0,0]*(Sqrt[2*Ip[0,F0,0]])))+(1/((γ[t,F0,0])^2))+1)-√((((px^2)+(po^2))/(2*Ip[0,F0,0]))-((2*px)/(γ[t,F0,0]*(Sqrt[2*Ip[0,F0,0]])))+(1/((γ[t,F0,0])^2))+1));

fBondar[t_,F0_,px_,po_]:=(1+(1/(2*((γ[t,F0,0])^2)))+(((px^2)+(po^2))/(2*Ip[0,F0,0])))*ArcCosh[α[t,F0,px,po]]-(Sqrt[((α[t,F0,px,po])^2)-1])*(((β[t,F0,px,po]*px)/(γ[t,F0,0]))*(Sqrt[(2/(Ip[0,F0,0]))])+((α[t,F0,px,po]*(1-2*((β[t,F0,px,po])^2)))/(2*((γ[t,F0,0])^2))));

BondarProbab[t_,F0_,ϵ_, px_, py_, pz_]:=Exp[-((2*Ip[0,F0,0])/(ωAU))fBondar[t,F0,px,(Sqrt[(py^2)+(pz^2)])]];
BondarProbabNorm[t_,F0_,ϵ_, px_, py_, pz_]:= BondarProbab[t,F0,0, px, py, pz]/BondarProbab[0,F0,0, 0, 0,0];

Plots PPT vs Bondar

Steps = 150;

Time = 0;
Fconst = 0.06;
eps = 0.0000001;
PX = 0;
pxSpan = Array[#&, Steps, {-2.25,2.25}];
pySpan = Array[#&, Steps, {-2.25,2.25}];
PY = 0;
PZ = 0;
pzSpan = Array[#&, Steps, {-2.25,2.25}];
Ppy= Transpose[{pySpan,PPTProbabNorm[Time,Fconst,eps,PX ,pySpan,PZ]}];
Bpy = Transpose[{pySpan,BondarProbabNorm[Time,Fconst,eps,PX ,pySpan,PZ]}];
ListLinePlot[{Ppy,Bpy}, PlotRange->All, Frame->True, FrameLabel->{"py / au", "Ionisation Probability"},
PlotStyle->{Blue,{Red,Dashed}},
PlotLegends->{"PPT","Bondar"}]

Ppz= Transpose[{pySpan,PPTProbabNorm[Time,Fconst,eps,PX ,p0[Fconst,eps],pzSpan]}];
Bpz = Transpose[{pxSpan,BondarProbabNorm[Time,Fconst,eps,PX ,PY,pzSpan]}];ListLinePlot[{Ppz,Bpz}, PlotRange->All, Frame->True, FrameLabel->{"pz / au", "Ionisation Probability"},
PlotStyle->{Blue,{Red,Dashed}},
PlotLegends->{"PPT", "Bondar"}]

Ppx = Transpose[{pySpan,PPTProbabNorm[Time,Fconst,eps,pxSpan ,p0[Fconst,eps],PZ]}];
Bpx = Transpose[{pxSpan,BondarProbabNorm[Time,Fconst,eps,pxSpan ,0,PZ]}];ListLinePlot[{Ppx,Bpx}, PlotRange->All, Frame->True, FrameLabel->{"px / au", "Ionisation Probability"},
PlotStyle->{Blue,{Red,Dashed}},
PlotLegends->{"PPT", "Bondar"}]

TSpan =  Array[((#)/(ωAU))&, Steps, {((-π)/3),(π/3)}];
PX = 0;
PZ = 0;

γ[0,Fconst,eps]

PT = {{TSpan[[1]],PPTProbabNorm[TSpan[[1]],Fconst,eps,PX ,p0[Fconst,eps],PZ]}};
For[i=2,i<Steps+1,i++,
AppendTo[PT,
{TSpan[[i]],PPTProbabNorm[TSpan[[i]],Fconst,eps,PX ,p0[Fconst,eps],PZ]}
]
]

BT = {{TSpan[[1]],BondarProbabNorm[TSpan[[1]],Fconst,eps,PX ,0,PZ]}};
For[i=2,i<Steps+1,i++,
AppendTo[BT,
{TSpan[[i]],BondarProbabNorm[TSpan[[i]],Fconst,eps,PX ,0,PZ]}
]
]

KT=Transpose[{TSpan,Exp[((- 2(((2*Ip[0,Fconst,eps]))^(3/2)))/(3*F[TSpan,Fconst,eps]))]/Exp[((- 2(((2*Ip[0,Fconst,eps]))^(3/2)))/(3*F[0,Fconst,eps]))]}];
ListLinePlot[{PT,BT,KT}, PlotRange->All, Frame->True, 
PlotStyle->{Blue,{Red,Dashed},Black},FrameLabel->{"time / au", "Ionisation Probability"},
PlotLegends->{"PPT", "Bondar","Keldysh"}]
(*  *)

epsDensity = 0.6;
Pdensity = Transpose[{ConstantArray[pxSpan[[1]],Steps],pySpan,PPTProbabNorm[Time,Fconst,epsDensity,pxSpan[[1]] ,pySpan,PZ]}];
For[i=2,i<Steps+1,i++,
Pdensity=Join[Pdensity,
Transpose[{ConstantArray[pxSpan[[i]],Steps],pySpan,PPTProbabNorm[Time,Fconst,epsDensity,pxSpan[[i]] ,pySpan,PZ]}]
]
]
Pdensity[[All,3]]=Pdensity[[All,3]]*0.7+0.3; (* to get nice colours in the plot *)

ListDensityPlot[Pdensity,ColorFunction->(Hue[-#]&), ColorFunctionScaling->False, 
PlotRange->{{-2.25,2.25},{-2.25,2.25},All},
Axes->True,
AxesLabel->{"px / au", "py / au"}]