Probing the vacuum fluctuations of a scalar ghost-free theory 
Jens Boos, Valeri P. Frolov, and Andrei Zelnikov. Last edited by Jens Boos, Jan 14, 2019.

(* import CustomTicks package which can be downloaded here: http://github.com/mark-caprio/CustomTicks *)
<<CustomTicks`
SetDirectory[NotebookDirectory[]];

Plot of the local renormalized vacuum polarization

PhiLocal[ϖ_,x_,λ_,m_]:=λ/(4π)(2ϖ Sin[2ϖ Abs[x]]-λ Cos[2ϖ x])/(Sqrt[ϖ^2+m^2](4ϖ^2+λ^2));
PhiSquaredLocal[x_,λ_,m_]:=NIntegrate[PhiLocal[ϖ,x,λ,m],{ϖ,0,10}];

λTMP:=0.5;
p01[size_] := Plot[{PhiSquaredLocal[x,λTMP,0.05],PhiSquaredLocal[x,λTMP,0.1],PhiSquaredLocal[x,λTMP,0.5]},{x,0,5},
PlotRange->{-0.2,0},
PlotStyle->{{Thickness[0.0025]},{Thickness[0.0025], Dashing[0.03]},{Thickness[0.0025], Dashing[0.01]}},
FrameStyle->{{Gray,AbsoluteThickness[2(size/800)]},{Gray,AbsoluteThickness[2(size/800)]},
			{Gray,AbsoluteThickness[2(size/800)]},{Gray,AbsoluteThickness[2(size/800)]}},
FrameTicksStyle->Directive[FontSize->25(size/800)],Frame->True,
FrameTicks->{{LinTicks[-0.2,0,TickLabelStep->2,MajorTickLength->0.01,MinorTickLength->0.005],None},
			{LinTicks[0,5,TickLabelStep->1,MajorTickLength->0.01,MinorTickLength->0.005],None}},
GridLines->Automatic,
ImageSize->size];
p01[500]
Export["plot-phisquared-loc.png",p01[2000]];

Plot of the vacuum polarization in the GF1 case
(logarithmically divergent term subtracted)

bbGF1[ϖ_,ℓ_,x_]:=-Im[Exp[I ϖ x]Erf[x/(2ℓ)+I ϖ ℓ]];
ccGF1[ϖ_,ℓ_,λ_]:=2ϖ/λ-Erfi[ϖ ℓ];
PhiGF1ren[ϖ_,ℓ_,x_,λ_,m_]:=1/(4π Sqrt[ϖ^2+m^2])((bbGF1[ϖ,ℓ,x]^2-Cos[ϖ x]^2-2Cos[ϖ x]bbGF1[ϖ,ℓ,x]ccGF1[ϖ,ℓ,λ])/(1+ccGF1[ϖ,ℓ,λ]^2)-(Exp[-x^2/(2ℓ^2)]-2Exp[-x^2/(4ℓ^2)]Cos[ϖ x]));
PhiSquaredGF1ren[x_,ℓ_,λ_,m_]:=NIntegrate[PhiGF1ren[ϖ,ℓ,x,λ,m],{ϖ,0,7},AccuracyGoal->5]-1/(2π)BesselK[0,m x]Exp[-x^2/(4ℓ^2)];

λTMP:=0.5;mTMP:=0.01;ℓTMP:=1;
p02[size_] := ListPlot[{
Table[{x,PhiSquaredLocal[x,λTMP,mTMP]},{x,0,5,0.05}],
Table[{x,PhiSquaredGF1ren[x,ℓTMP,λTMP,mTMP]},{x,0,5,0.05}]
},PlotRange->{-0.6,-0.1},Joined->True,PlotMarkers->None,	PlotStyle->{{Thickness[0.0025]},{Thickness[0.0025],Dashing[0.01]},{Thickness[0.0025], Dashing[0.03]},{Thickness[0.0025], DotDashed}},
FrameStyle->{{Gray,AbsoluteThickness[2(size/800)]},{Gray,AbsoluteThickness[2(size/800)]},
			{Gray,AbsoluteThickness[2(size/800)]},{Gray,AbsoluteThickness[2(size/800)]}},
FrameTicksStyle->Directive[FontSize->25(size/800)],Frame->True,
FrameTicks->{{LinTicks[-0.6,-0.1,TickLabelStep->2,TickLabelStart->1, MajorTickLength->0.01,MinorTickLength->0.005],None},
			{LinTicks[0,5,TickLabelStep->1,MajorTickLength->0.01,MinorTickLength->0.005],None}},
GridLines->Automatic,
ImageSize->size];
p02[500]
Export["plot-phisquared-loc-gf1.png",p02[2000]];

Plot of (Subscript[f, ϖ]) for different values of ϖℓ

Fomega[ξ_,ϖ_,ℓ_]:=(1-Exp[-ϖ^4ℓ^4(1-ξ^2)^2])/(1-ξ^2);
p03[size_] := Plot[{5Fomega[ξ,0.5,1],Fomega[ξ,2,1]},{ξ,0,10},MaxRecursion->10,
PlotRange->{-3,3},		PlotStyle->{{Thickness[0.0025]},{Thickness[0.0025], Dashing[0.005]},{Thickness[0.0025], Dashing[0.03]},{Gray,Thickness[0.0025], Dashing[0.01]}},
FrameStyle->{{Gray,AbsoluteThickness[2(size/800)]},{Gray,AbsoluteThickness[2(size/800)]},
			{Gray,AbsoluteThickness[2(size/800)]},{Gray,AbsoluteThickness[2(size/800)]}},
FrameTicksStyle->Directive[FontSize->25(size/800)],Frame->True,
FrameTicks->{{LinTicks[-3.2,3.2,TickLabelStep->2,MajorTickLength->0.01,MinorTickLength->0.005],None},
			{LinTicks[0,10,TickLabelStep->1,MajorTickLength->0.01,MinorTickLength->0.005],None}},
GridLines->Automatic,
ImageSize->size];
p03[500]
Export["plot-fomega.png",p03[2000]];

Plot of the vacuum polarization in the GF2 case

(* auxiliary quantities *)
b[ϖ_,ℓ_]:=(ϖ ℓ)^2;
xtilde[ϖ_,x_]:=ϖ x;
ξp[a2_,ϖ_,ℓ_]:=Sqrt[1+2a2/b[ϖ,ℓ]];
ξm[a2_,ϖ_,ℓ_]:=Sqrt[1-2a2/b[ϖ,ℓ]];
c1[a2_,ϖ_,ℓ_]:=Exp[-4a2^2-b[ϖ,ℓ]^2(1-ξp[a2,ϖ,ℓ])];
c3[a2_,ϖ_,ℓ_]:=Exp[-4a2^2+b[ϖ,ℓ]^2(1-ξm[a2,ϖ,ℓ])];
mm[a2_,ϖ_,ℓ_]:=1/2(1-Exp[-4a2^2])(1/2-b[ϖ,ℓ]^2/a2^2);
nn[a2_,ϖ_,ℓ_]:=1/2(1-Exp[-4a2^2])(b[ϖ,ℓ]^2/a2^2-1);

(* special exponential integral taking care of proper branch cut *)
Ei[z_]:=Piecewise[{
{ExpIntegralEi[z],Re[z]≥0},
{ExpIntegralEi[z]+I π,Re[z]<0}
}];

(* approximation of the main integrand and its approximation *)
F[ξ_,ϖ_,ℓ_,x_]:=Cos[xtilde[ϖ,x]ξ](1-Exp[-b[ϖ,ℓ]^2(1-ξ^2)^2])/(1-ξ^2);
Fapprox[ξ_,a2_,ϖ_,ℓ_,x_]:=Cos[xtilde[ϖ,x]ξ]Piecewise[{
{0,b[ϖ,ℓ]==0},
{-1/(1+ξ^2),b[ϖ,ℓ]≤3a2},
{Piecewise[{
{(1-c1[a2,ϖ,ℓ]Exp[-b[ϖ,ℓ]^2(1-ξ)])/(1-ξ^2),ξ<ξm[a2,ϖ,ℓ]},
{mm[a2,ϖ,ℓ]ξ+nn[a2,ϖ,ℓ],ξ≥ξm[a2,ϖ,ℓ]&&ξ≤ξp[a2,ϖ,ℓ]},
{(1-c3[a2,ϖ,ℓ]Exp[-b[ϖ,ℓ]^2(ξ-1)])/(1-ξ^2),ξ>ξp[a2,ϖ,ℓ]}
}],b[ϖ,ℓ]>3a2}
}];
Difference[ξ_,a2_,ϖ_,ℓ_,x_]:=F[ξ,ϖ,ℓ,x]-Fapprox[ξ,a2,ϖ,ℓ,x];

(* Analytic expressions *)
AI1[a2_,ϖ_,ℓ_,x_]:=-π/2Exp[-Abs[xtilde[ϖ,x]]];
AI2[a2_,ϖ_,ℓ_,x_]:=1/2Cos[xtilde[ϖ,x]](CosIntegral[xtilde[ϖ,x](1+ξm[a2,ϖ,ℓ])]-CosIntegral[xtilde[ϖ,x](1-ξm[a2,ϖ,ℓ])]-CosIntegral[xtilde[ϖ,x](1+ξp[a2,ϖ,ℓ])]+CosIntegral[xtilde[ϖ,x](ξp[a2,ϖ,ℓ]-1)])+1/xtilde[ϖ,x]^2(mm[a2,ϖ,ℓ](Cos[xtilde[ϖ,x]ξp[a2,ϖ,ℓ]]-Cos[xtilde[ϖ,x]ξm[a2,ϖ,ℓ]])+(mm[a2,ϖ,ℓ]ξp[a2,ϖ,ℓ]+nn[a2,ϖ,ℓ])xtilde[ϖ,x]Sin[xtilde[ϖ,x]ξp[a2,ϖ,ℓ]]-(mm[a2,ϖ,ℓ]ξm[a2,ϖ,ℓ]+nn[a2,ϖ,ℓ])xtilde[ϖ,x]Sin[xtilde[ϖ,x]ξm[a2,ϖ,ℓ]])+1/2Sin[xtilde[ϖ,x]](π+SinIntegral[xtilde[ϖ,x](1+ξm[a2,ϖ,ℓ])]-SinIntegral[xtilde[ϖ,x](1-ξm[a2,ϖ,ℓ])]-SinIntegral[xtilde[ϖ,x](1+ξp[a2,ϖ,ℓ])]+SinIntegral[xtilde[ϖ,x](1-ξp[a2,ϖ,ℓ])])+c1[a2,ϖ,ℓ]/2Re[Exp[I xtilde[ϖ,x]]ExpIntegralEi[(-b[ϖ,ℓ]^2-I xtilde[ϖ,x])(1-ξm[a2,ϖ,ℓ])]-Exp[-I xtilde[ϖ,x]-2b[ϖ,ℓ]^2]ExpIntegralEi[(b[ϖ,ℓ]^2+I xtilde[ϖ,x])(1+ξm[a2,ϖ,ℓ])]]-c1[a2,ϖ,ℓ]/2Re[Exp[I xtilde[ϖ,x]]ExpIntegralEi[-b[ϖ,ℓ]^2-I xtilde[ϖ,x]]-Exp[-I xtilde[ϖ,x]-2b[ϖ,ℓ]^2]ExpIntegralEi[b[ϖ,ℓ]^2+I xtilde[ϖ,x]]]-c3[a2,ϖ,ℓ]/2Re[Exp[-I xtilde[ϖ,x]]Ei[(b[ϖ,ℓ]^2+I xtilde[ϖ,x])(1-ξp[a2,ϖ,ℓ])]-Exp[-I xtilde[ϖ,x]+2b[ϖ,ℓ]^2]Ei[-(b[ϖ,ℓ]^2+I xtilde[ϖ,x])(1+ξp[a2,ϖ,ℓ])]];
AI20[a2_,ϖ_,ℓ_,x_]:=1/2mm[a2,ϖ,ℓ](ξp[a2,ϖ,ℓ]^2-ξm[a2,ϖ,ℓ]^2)+nn[a2,ϖ,ℓ](ξp[a2,ϖ,ℓ]-ξm[a2,ϖ,ℓ])-1/2Log[(1-ξm[a2,ϖ,ℓ])/(1+ξm[a2,ϖ,ℓ])]+c1[a2,ϖ,ℓ]/2(ExpIntegralEi[b[ϖ,ℓ]^2(ξm[a2,ϖ,ℓ]-1)]-ExpIntegralEi[-b[ϖ,ℓ]^2]+Exp[-2b[ϖ,ℓ]^2]ExpIntegralEi[b[ϖ,ℓ]^2]-Exp[-2b[ϖ,ℓ]^2]ExpIntegralEi[b[ϖ,ℓ]^2(1+ξm[a2,ϖ,ℓ])])+1/2Log[(ξp[a2,ϖ,ℓ]-1)/(ξp[a2,ϖ,ℓ]+1)]-c3[a2,ϖ,ℓ]/2(ExpIntegralEi[b[ϖ,ℓ]^2(1-ξp[a2,ϖ,ℓ])]-Exp[2b[ϖ,ℓ]^2]ExpIntegralEi[-b[ϖ,ℓ]^2(1+ξp[a2,ϖ,ℓ])]);
AI[a2_,ϖ_,ℓ_,x_]:=Piecewise[{
{0,b[ϖ,ℓ]==0},
{1/π AI1[a2,ϖ,ℓ,x],b[ϖ,ℓ]≤3a2},
{1/π Piecewise[{{AI20[a2,ϖ,ℓ,x],x==0},{AI2[a2,ϖ,ℓ,x],x≠0}}],b[ϖ,ℓ]>3a2}
}];

(* numerical integral of the remainder *)
NI[a2_,ϖ_?NumericQ,ℓ_?NumericQ,x_]:=NI[a2,ϖ,ℓ,x]=Piecewise[{
{NIntegrate[1/π Difference[ξ,a2,ϖ,ℓ,x],{ξ,0,20/(1+ϖ ℓ)},AccuracyGoal->5],b[ϖ,ℓ]≤3a2},
{NIntegrate[1/π Difference[ξ,a2,ϖ,ℓ,x],{ξ,0.9(1-1/(ϖ ℓ)),1.1(1+3/(ϖ ℓ)^2)},AccuracyGoal->5],b[ϖ,ℓ]>3a2}
}];

(* overall expression (Subscript[g, ω])(x) (we call it 𝒩 here for historical reasons in our notation ) *)
𝒩[a2_,ϖ_,ℓ_,x_]:=AI[a2,ϖ,ℓ,x]+NI[a2,ϖ,ℓ,x];

(* analytic result for (Subscript[g, ω])(0) *)
𝒩0[ϖ_,ℓ_]:=Sqrt[b[ϖ,ℓ]]/π(Sqrt[2]π b[ϖ,ℓ]/(6Gamma[3/4])HypergeometricPFQ[{3/4,5/4},{3/2,7/4},-b[ϖ,ℓ]^2]-Gamma[3/4]HypergeometricPFQ[{1/4,3/4},{1/2,5/4},-b[ϖ,ℓ]^2] );

(* vacuum polarization, finally *)
psA[ϖ_,x_]:=Cos[ϖ x];
psB[a2_,ϖ_?NumericQ,ℓ_?NumericQ,x_]:=psB[a2,ϖ,ℓ,x]=2𝒩[a2,ϖ,ℓ,x]-Sin[ϖ x];
psC[ϖ_,ℓ_,λ_]:=2𝒩0[ϖ,ℓ]+2ϖ/λ;
PhiSquaredGF2[a2_,x_,ℓ_,λ_,m_]:=NIntegrate[1/(4π Sqrt[ϖ^2+m^2])(psB[a2,ϖ,ℓ,x]^2-psA[ϖ,x]^2-2psA[ϖ,x]psB[a2,ϖ,ℓ,x]psC[ϖ,ℓ,λ])/(1+psC[ϖ,ℓ,λ]^2),{ϖ,0,5},AccuracyGoal-> 5]

a2TMP:=0.5604532215;ℓTMP:=1;mTMP:=0.01;
p04[size_] := ListPlot[{
Table[{x,PhiSquaredLocal[x,0.5,mTMP]},{x,0,10.1,0.251}],
Table[{x,PhiSquaredLocal[x,2,mTMP]},{x,0,10.1,0.251}],
Table[{x,PhiSquaredGF2[a2TMP,x,ℓTMP,0.5,mTMP]},{x,0,10.1,0.251}],
Table[{x,PhiSquaredGF2[a2TMP,x,ℓTMP,2,mTMP]},{x,0,10.1,0.251}]
},PlotRange->{-0.45,-0.1},Joined->True,PlotMarkers->None,	PlotStyle->{{Thickness[0.0025]},{Thickness[0.0025], Dashing[0.03]},{Thickness[0.0025], Dashing[0.01]},{Thickness[0.0025], DotDashed}},
FrameStyle->{{Gray,AbsoluteThickness[2(size/800)]},{Gray,AbsoluteThickness[2(size/800)]},
			{Gray,AbsoluteThickness[2(size/800)]},{Gray,AbsoluteThickness[2(size/800)]}},
FrameTicksStyle->Directive[FontSize->25(size/800)],Frame->True,
FrameTicks->{{LinTicks[-0.45,-0.1,TickLabelStep->1,MajorTickLength->0.01,MinorTickLength->0.005],None},
			{LinTicks[0,10,TickLabelStep->1,MajorTickLength->0.01,MinorTickLength->0.005],None}},
GridLines->Automatic,
ImageSize->size];
p04[500]
Export["plot-phisquared-loc-gf2.png",p04[2000]];

p05[size_] := ListPlot[{
Table[{x,PhiSquaredLocal[x,0.5,mTMP]},{x,0,1,0.0251}],
Table[{x,PhiSquaredLocal[x,2,mTMP]},{x,0,1,0.0251}],
Table[{x,PhiSquaredGF2[a2TMP,x,ℓTMP,0.5,mTMP]},{x,0,1,0.0251}],
Table[{x,PhiSquaredGF2[a2TMP,x,ℓTMP,2,mTMP]},{x,0,1,0.0251}]
},PlotRange->{-0.24,-0.51},Joined->True,PlotMarkers->None,	PlotStyle->{{Thickness[0.0025]},{Thickness[0.0025], Dashing[0.03]},{Thickness[0.0025], Dashing[0.01]},{Thickness[0.0025], DotDashed}},
FrameStyle->{{Gray,AbsoluteThickness[2(size/800)]},{Gray,AbsoluteThickness[2(size/800)]},
			{Gray,AbsoluteThickness[2(size/800)]},{Gray,AbsoluteThickness[2(size/800)]}},
FrameTicksStyle->Directive[FontSize->25(size/800)],Frame->True,
FrameTicks->{{LinTicks[-0.5,-0.25,TickLabelStep->1,MajorTickLength->0.01,MinorTickLength->0.005],None},
			{LinTicks[0,1,TickLabelStep->1,MajorTickLength->0.01,MinorTickLength->0.005],None}},
GridLines->Automatic,
ImageSize->size];
p05[500]
Export["plot-phisquared-loc-gf2-detail.png",p05[2000]];

Plot of the difference of <(φ^2)(0)> in the GF2 and local case

p06[size_] := ListPlot[{
Table[{λ,PhiSquaredGF2[a2TMP,0,ℓTMP,λ,0.01]-PhiSquaredLocal[0,λ,0.01]},{λ,0.02,5,0.05321}],
Table[{λ,PhiSquaredGF2[a2TMP,0,ℓTMP,λ,0.1]-PhiSquaredLocal[0,λ,0.1]},{λ,0.02,5,0.05321}],
Table[{λ,PhiSquaredGF2[a2TMP,0,ℓTMP,λ,0.5]-PhiSquaredLocal[0,λ,0.5]},{λ,0.02,5,0.05321}],
Table[{λ,PhiSquaredGF2[a2TMP,0,ℓTMP,λ,1]-PhiSquaredLocal[0,λ,1]},{λ,0.02,5,0.05321}]
},PlotRange->{0,0.035},Joined->True,PlotMarkers->None,	PlotStyle->{{Thickness[0.0025]},{Thickness[0.0025], Dashing[0.03]},{Thickness[0.0025], Dashing[0.01]},{Thickness[0.0025], DotDashed}},
FrameStyle->{{Gray,AbsoluteThickness[2(size/800)]},{Gray,AbsoluteThickness[2(size/800)]},
			{Gray,AbsoluteThickness[2(size/800)]},{Gray,AbsoluteThickness[2(size/800)]}},
FrameTicksStyle->Directive[FontSize->25(size/800)],Frame->True,
FrameTicks->{{LinTicks[0,0.035,TickLabelStep->1,MajorTickLength->0.01,MinorTickLength->0.005],None},
			{LinTicks[0,5,TickLabelStep->1,MajorTickLength->0.01,MinorTickLength->0.005],None}},
GridLines->Automatic,
ImageSize->size];
p06[500]
Export["plot-phisquared0-comparison.png",p06[2000]];

Asymptotics for analytic result for <(φ^2)(0)> in the GF2 case

p07[size_] := Plot[{𝒩0[ϖ,ℓTMP],-Gamma[3/4]Sqrt[b[ϖ,ℓTMP]]/π,-1/(4Sqrt[π]b[ϖ,ℓTMP])},{ϖ,0,3},MaxRecursion->10,
PlotRange->{-0.25,0.05},		PlotStyle->{{Thickness[0.0025]},{Thickness[0.0025], Dashing[0.03]},{Thickness[0.0025], Dashing[0.01]},{Thickness[0.0025], Dashing[0.01]}},
FrameStyle->{{Gray,AbsoluteThickness[2(size/800)]},{Gray,AbsoluteThickness[2(size/800)]},
			{Gray,AbsoluteThickness[2(size/800)]},{Gray,AbsoluteThickness[2(size/800)]}},
FrameTicksStyle->Directive[FontSize->25(size/800)],Frame->True,
FrameTicks->{{LinTicks[-0.25,0.5,TickLabelStep->1,MajorTickLength->0.01,MinorTickLength->0.005],None},
			{LinTicks[0,3,TickLabelStep->2,MajorTickLength->0.01,MinorTickLength->0.005],None}},
GridLines->Automatic,
ImageSize->size];
p07[500]
Export["plot-n-x=0-exact.png",p07[2000]];

Asymptotics for approximate result for <(φ^2)(x)> in the GF2 case

a2TMP:=0.5604532215;ℓTMP:=1;
p08[size_] := ListPlot[{
Table[{ϖ,2𝒩[a2TMP,ϖ,ℓTMP,0.0]-Sin[0.0ϖ]},{ϖ,0.1,15,0.1251}],
Table[{ϖ,2𝒩[a2TMP,ϖ,ℓTMP,0.5]-Sin[0.5ϖ]},{ϖ,0.1,15,0.1251}],
Table[{ϖ,2𝒩[a2TMP,ϖ,ℓTMP,1.0]-Sin[1.0ϖ]},{ϖ,0.1,15,0.1251}],
Table[{ϖ,2𝒩[a2TMP,ϖ,ℓTMP,2.0]-Sin[2.0ϖ]},{ϖ,0.1,15,0.1251}]
},PlotRange->{-0.9,0.5},Joined->True,PlotMarkers->None,	PlotStyle->{{Thickness[0.0025]},{Thickness[0.0025], Dashing[0.03]},{Thickness[0.0025], Dashing[0.01]},{Thickness[0.0025], DotDashed}},
FrameStyle->{{Gray,AbsoluteThickness[2(size/800)]},{Gray,AbsoluteThickness[2(size/800)]},
			{Gray,AbsoluteThickness[2(size/800)]},{Gray,AbsoluteThickness[2(size/800)]}},
FrameTicksStyle->Directive[FontSize->25(size/800)],Frame->True,
FrameTicks->{{LinTicks[-0.8,0.5,TickLabelStep->2,MajorTickLength->0.01,MinorTickLength->0.005],None},
			{LinTicks[0,15,TickLabelStep->1,MajorTickLength->0.01,MinorTickLength->0.005],None}},
GridLines->Automatic,
ImageSize->size];
p08[500]
Export["plot-n-numerical.png",p08[2000]];

Visualization of numerical integration procedure in the GF2 case

(* b < b0: visualization of the difference between the original function and the improved approximation *)
a2TMP:=0.5604532215;
FF[ξ_,b_]:=(1-Exp[-b^2(1-ξ^2)^2])/(1-ξ^2);
Tail[ξ_,x_]:=-Cos[ξ x]/(1+ξ^2);
p09a[bTMP_,a2TMP_,xTMP_,size_] := Plot[{Cos[ξ xTMP]FF[ξ,bTMP],Cos[ξ xTMP]FF[ξ,bTMP]-Tail[ξ,xTMP],Tail[ξ,xTMP]},{ξ,0,10},MaxRecursion->10,
PlotRange->{-0.11,0.21},		PlotStyle->{{Thickness[0.0025]},{Thickness[0.0025]},{Thickness[0.0025], Dashing[0.005]},{Gray,Thickness[0.0025], Dashing[0.03]}},
FrameStyle->{{Gray,AbsoluteThickness[2(size/800)]},{Gray,AbsoluteThickness[2(size/800)]},
			{Gray,AbsoluteThickness[2(size/800)]},{Gray,AbsoluteThickness[2(size/800)]}},
FrameTicksStyle->Directive[FontSize->25(size/800)],Frame->True,
FrameTicks->{{LinTicks[-0.1,0.2,TickLabelStep->2,MajorTickLength->0.01,MinorTickLength->0.005],None},
			{LinTicks[0,10,TickLabelStep->1,MajorTickLength->0.01,MinorTickLength->0.005],None}},
GridLines->Automatic,
Filling->{2->{Axis,Orange}},
ImageSize->size];
p09b[bTMP_,a2TMP_,xTMP_,size_] := Plot[{Cos[ξ xTMP]FF[ξ,bTMP],Cos[ξ xTMP]FF[ξ,bTMP]-Tail[ξ,xTMP],Tail[ξ,xTMP]},{ξ,0,10},MaxRecursion->10,
PlotRange->All,		PlotStyle->{{Thickness[0.0025]},{Thickness[0.0025]},{Thickness[0.0025], Dashing[0.005]},{Thickness[0.0025]}},
FrameStyle->{{Gray,AbsoluteThickness[2(size/800)]},{Gray,AbsoluteThickness[2(size/800)]},
			{Gray,AbsoluteThickness[2(size/800)]},{Gray,AbsoluteThickness[2(size/800)]}},
FrameTicksStyle->Directive[FontSize->25(size/800)],Frame->True,
FrameTicks->{{LinTicks[-1.1,1.1,TickLabelStep->1,MajorTickLength->0.01,MinorTickLength->0.005],None},
			{LinTicks[0,10,TickLabelStep->1,MajorTickLength->0.01,MinorTickLength->0.005],None}},
GridLines->Automatic,
Filling->{2->{Axis,Orange}},
ImageSize->size];
Export["plot-b=0.1-x=0.png",p09a[0.1,a2TMP,0,2000]];
Export["plot-b=0.1-x=10.png",p09b[0.1,a2TMP,10,2000]];
p09a[0.1,a2TMP,0,500]
p09b[0.1,a2TMP,10,500]

(* b > b0: visualization of the difference between the original function and the improved approximation *)
a2TMP:=0.5604532215;
ξξm[b_,a2_]:=Sqrt[1-2a2/b];
ξξp[b_,a2_]:=Sqrt[1+2a2/b];
FFlinear[ξ_,b_,a2_]:=1/2(1-Exp[-4a2^2])((1/2-b^2/a2^2)ξ+b^2/a2^2-1);
cm[b_,a2_]:=Exp[-4a2^2+b^2(1-Sqrt[1-2a2/b])];
cp[b_,a2_]:=Exp[-4a2^2-b^2(1-Sqrt[1+2a2/b])];
FFapprox[ξ_,b_,a2_]:=Piecewise[{
{F[ξ,b,a2]+1/(1+ξ^2),b≤5 a2},
{Piecewise[{
{(1-cm[b,a2]Exp[-b^2(1-ξ)])/(1-ξ^2),ξ<ξξm[b,a2]},
{FFlinear[ξ,b,a2],ξ≥ξξm[b,a2]&&ξ≤ξξp[b,a2]},
{(1-cp[b,a2]Exp[-b^2(ξ-1)])/(1-ξ^2),ξ>ξξp[b,a2]}
}],b>5 a2}
}];

p10[bTMP_,a2TMP_,xTMP_,size_] := Plot[{ Cos[ξ xTMP]FF[ξ,bTMP],Cos[ξ xTMP](FF[ξ,bTMP]-FFapprox[ξ,bTMP,a2TMP]),Cos[ξ xTMP]FFapprox[ξ,bTMP,a2TMP]},{ξ,0,3},MaxRecursion->10,
PlotRange->{-2.2,2.2},		PlotStyle->{{Thickness[0.0025]},{Thickness[0.0025]},{Thickness[0.0025], Dashing[0.005]},{Thickness[0.0025]}},
FrameStyle->{{Gray,AbsoluteThickness[2(size/800)]},{Gray,AbsoluteThickness[2(size/800)]},
			{Gray,AbsoluteThickness[2(size/800)]},{Gray,AbsoluteThickness[2(size/800)]}},
FrameTicksStyle->Directive[FontSize->25(size/800)],Frame->True,
FrameTicks->{{LinTicks[-2,2,TickLabelStep->1,MajorTickLength->0.01,MinorTickLength->0.005],None},
			{LinTicks[0,3,TickLabelStep->1,MajorTickLength->0.01,MinorTickLength->0.005],None}},
GridLines->Automatic,
Filling->{2->{Axis,Orange}},
ImageSize->size];
p10[3,a2TMP,0,500]
p10[3,a2TMP,10,500]
Export["plot-b=3-x=0.png",p10[3,a2TMP,0,2000]];
Export["plot-b=3-x=10.png",p10[3,a2TMP,10,2000]];