Script for the paper “Uniform enclosures for the phase and zeros of Bessel functions and their derivatives”

https://michaellevitin.net/bessels.html

Nikolay Filonov, Michael Levitin, Iosif Polterovich, and David A. Sher

Behind the scenes

Directories

curdir=SetDirectory[NotebookDirectory[]];
SaveDir="./";

Colors and symbols

mycolors=ColorData[97,"ColorList"][[1;;8]]
mytexcolors="\\definecolor{mycolors1}{rgb}{0.368417, 0.506779, 0.709798}\\definecolor{mycolors2}{rgb}{0.880722, 0.611041, 0.142051}\\definecolor{mycolors3}{rgb}{0.560181, 0.691569, 0.194885}";
{filledcircle,filledsquare,filleddiamond,filleduptriangle,filleddowntriangle,circle,square,diamond,uptriangle,downtriangle}=Graphics`PlotMarkers[][[All,1]]

MaTeX and lines

Needs["MaTeX`"]
SetOptions[MaTeX,"BasePreamble"->{"\\usepackage{amsmath}","\\usepackage{fourier}","\\usepackage{ebgaramond}","\\usepackage{bm}",
"\\usepackage{xcolor}\\definecolor{darkgreen}{rgb}{0.00, 0.67, 0.00}"(*,"\\usepackage{accents}",
"\\DeclareMathAccent{\\wtilde}{\\mathord}{largesymbols}{\"65}",
"\\DeclareRobustCommand{\\utilde}[1]{\\underaccent{\\wtilde}{#1}}}"*)},FontSize->12,Magnification->1];

Thk=AbsoluteThickness[1.8]; 
VeryThk=AbsoluteThickness[3.]; 
Thn=AbsoluteThickness[0.6]; 
res= 600; size = 4.00; sizesmall=2.50; sizewide=5.00;
hwcoeff=1.2;
Imgsize ={UpTo[Round[72 size]], UpTo[Round[72 hwcoeff size]]};
Imgsizesmall ={UpTo[Round[72 sizesmall]], UpTo[Round[72  hwcoeff sizesmall]]};
Imgsizewide ={UpTo[Round[72 sizewide]], UpTo[Round[72  hwcoeff sizewide]]};
lightgr=GrayLevel[.9];darkgr=GrayLevel[.6]; 
mydashing0=Charting`ResolvePlotTheme["Monochrome",Plot][[7]][[2]][[5]][[2]];
mydashing=Table[Directive[mydashing0[[j,2]], mydashing0[[j,4]]],{j,1,8}];
mydsh=mydashing[[2]]; mydot=mydashing[[3]];
PrTicks[xs_, ts_]:=Thread[{xs, MaTeX[ts]}]

mods for [Ho17]

Mods[a_,b_]:=a-b Round[a/b];

findallRoots from https://mathematica.stackexchange.com/questions/16439/find-all-roots-of-an-interpolating-function-solution-to-a-differential-equation/16444#16444

Clear[findAllRoots]
SyntaxInformation[findAllRoots]={"LocalVariables"->{"Plot",{2,2}},"ArgumentsPattern"->{_,_,OptionsPattern[]}};
SetAttributes[findAllRoots,HoldAll];

Options[findAllRoots]=Join[{"ShowPlot"->False,PlotRange->All},FilterRules[Options[Plot],Except[PlotRange]]];

findAllRoots[fn_,{l_,lmin_,lmax_},opts:OptionsPattern[]]:=Module[{pl,p,x,localFunction,brackets},localFunction=ReleaseHold[Hold[fn]/. HoldPattern[l]:>x];
If[lmin≠lmax,pl=Plot[localFunction,{x,lmin,lmax},Evaluate@FilterRules[Join[{opts},Options[findAllRoots]],Options[Plot]]];
p=Cases[pl,Line[{x__}]:>x,Infinity];
If[OptionValue["ShowPlot"],Print[Show[pl,PlotLabel->"Finding roots for this function",ImageSize->200,BaseStyle->{FontSize->8}]]],p={}];
brackets=Map[First,Select[(*This Split trick pretends that two points on the curve are "equal" if the function values have _opposite _ sign.Pairs of such sign-changes form the brackets for the subsequent FindRoot*)Split[p,Sign[Last[#2]]==-Sign[Last[#1]]&],Length[#1]==2&],{2}];
x/. Apply[FindRoot[localFunction==0,{x,##1}]&,brackets,{1}]/. x->{}];

cleanCountourPlot  from https://mathematica.stackexchange.com/questions/3190/saner-alternative-to-contourplot-fill/3279#3279

cleanContourPlot[cp_Graphics]:=Module[{points,groups,regions,lines},groups=Cases[cp,{style__,g_GraphicsGroup}:>{{style},g},Infinity];
points=First@Cases[cp,GraphicsComplex[pts_,___]:>pts,Infinity];
regions=Table[Module[{group,style,polys,edges,cover,graph},{style,group}=g;
polys=Join@@Cases[group,Polygon[pt_,___]:>pt,Infinity];
edges=Join@@(Partition[#,2,1,1]&/@polys);
cover=Cases[Tally[Sort/@edges],{e_,1}:>e];
graph=Graph[UndirectedEdge@@@cover];
{Sequence@@style,FilledCurve[List/@Line/@First/@Map[First,FindEulerianCycle/@(Subgraph[graph,#]&)/@ConnectedComponents[graph],{3}]]}],{g,groups}];
lines=Cases[cp,_Tooltip,Infinity];
Graphics[GraphicsComplex[points,{regions,lines}],Sequence@@Options[cp]]]

cleanContourPlot[Legended[cp_Graphics,rest___]]:=Legended[cleanContourPlot[cp],rest]

§1. Introduction and main results

§1.1. Setup I: Bessel functions

Asymptotics of

θasympt = x-Pi/4(2ν+1)+(4ν^2-1)/(8x)+(4ν^2-1)(4ν^2-25)/(384x^3);

Figure

xtib={5,10}; ytib={-1,-0.5,0,0.5,1};
figCnutau=GraphicsColumn[{GraphicsRow[Table[Plot[Cos[Pi τ]BesselJ[ν,x]+Sin[Pi τ]BesselY[ν,x]/.τ->{0,1/4,1/2,3/4}//Evaluate, {x,0,10+ν+ν^(1/3)},PlotStyle->mycolors[[4;;7]], AxesLabel->{MaTeX["x"], None},Ticks->{PrTicks[xtib,xtib], PrTicks[xtib,xtib]},Epilog->{Inset[Framed[MaTeX["\\nu="<>ToString[ν]],RoundingRadius->3],{Right, Bottom},Scaled[{1.5,-0.25}]]}], {ν,{0,3}}]],
LineLegend[mycolors[[4;;7]], MaTeX[{"J_\\nu(x)", "\\mathcal{C}_{\\nu,\\frac{1}{4}}(x)","Y_\\nu(x)", "\\mathcal{C}_{\\nu,\\frac{3}{4}}(x)"}],LegendLayout->"Row"]}]

Export[SaveDir<>"figCnutau.pdf", figCnutau]

Computing  following [Ho17]

θ[ν_,x_]:= If[Abs[x]>=Abs[ν],Sqrt[x^2-ν^2] -ν  ArcCos[ν/x]-Pi/4,-Pi/4]+Mods[ArcTan[BesselJ[ν,x],BesselY[ν,x]]-If[Abs[x]>=Abs[ν],Sqrt[x^2-ν^2] -ν  ArcCos[ν/x]-Pi/4,-Pi/4], 2Pi];

Figure

yts1=Table[yt, {yt, -Pi/2, 3 Pi,Pi/2}]; 
xts1=Table[x,{x,2,10,2}];
yti1=PrTicks[yts1, {"-\\pi/2","0", "\\pi/2", "\\pi",  "3\\pi/2", "2\\pi",  "5\\pi/2", "3\\pi"}];
xti1=PrTicks[xts1,xts1];

plotθ=Plot[Table[θ[ν,x],{ν,{0,2,5}}]//Evaluate,{x,0,10}, PlotStyle->mycolors,Ticks->{xti1,yti1}, AxesLabel->MaTeX[{"x","\\theta_\\nu(x)"}],Epilog->{{Thin, Black,Dashed, Table[Line[{{0,yt},{10,yt}}], {yt, 0, 3 Pi,Pi/2}]}, {mycolors[[1]], Dashed,Table[Line[{{BesselJZero[0,k],0}, {BesselJZero[0,k], -Pi/2+Pi k}}],{k,1,3}],Table[Line[{{BesselYZero[0,k],0}, {BesselYZero[0,k], -Pi+Pi k}}],{k,1,3}]},
{mycolors[[1]],Table[Inset[filledcircle,{BesselJZero[0,k],0},Scaled[{1/2,1/2}]], {k,1,3}]},
{mycolors[[1]],Table[Inset[circle,{BesselYZero[0,k],0}], {k,1,3}]},
{mycolors[[2]], Dashed,Table[Line[{{BesselJZero[2,k],0}, {BesselJZero[2,k], -Pi/2+Pi k}}],{k,1,2}],Table[Line[{{BesselYZero[2,k],0}, {BesselYZero[2,k], -Pi+Pi k}}],{k,1,3}]},
{mycolors[[2]],Table[Inset[filledcircle,{BesselJZero[2,k],0}], {k,1,2}]},
{mycolors[[2]],Table[Inset[circle,{BesselYZero[2,k],0}], {k,1,3}]},
{mycolors[[3]], Dashed,Table[Line[{{BesselJZero[5,k],0}, {BesselJZero[5,k], -Pi/2+Pi k}}],{k,1,1}],Table[Line[{{BesselYZero[5,k],0}, {BesselYZero[5,k], -Pi+Pi k}}],{k,1,2}]},
{mycolors[[3]],Table[Inset[filledcircle,{BesselJZero[5,k],0}], {k,1,1}]},
{mycolors[[3]],Table[Inset[circle,{BesselYZero[5,k],0}], {k,1,2}]}
}, PlotLegends->MaTeX[{"\\nu=0","\\nu=2","\\nu=5"}]]

Export[SaveDir<>"plotthetabessels.pdf", plotθ]

Figure

figtaustar=Plot[1/Pi θ[ν,ν]+1/2,{ν,0,3},PlotRange->{0.0000,0.18},PlotStyle->{Thick,mycolors[[4]]}, Epilog->{mycolors[[4]],Thick, Line[{{0,0},{0,0.042732809954193596}}]}, (* some fiddling near ν==0*) AxesLabel->MaTeX[{"\\nu", "\\tau^*_\\nu"}],
Ticks->{{Range[3], MaTeX[Range[3]]}//Transpose, {0.05 Range[3], MaTeX[0.05 Range[3]]}//Transpose}]

Export[SaveDir<>"figtaustar.pdf", figtaustar]

§1.2. Setup II: derivatives of Bessel functions

Derivatives of Bessel functions

DBesselJ[ν_,x_]:=(1/2) (BesselJ[-1+ν,x]-BesselJ[1+ν,x]);
DBesselY[ν_,x_]:=(1/2) (BesselY[-1+ν,x]-BesselY[1+ν,x]);

Asymptotics of

ϕasympt = x-Pi/4(2ν-1)+(4ν^2+3)/(8x)+(16ν^4+184ν^2-63)/(384x^3);

Computing  following [Ho17]

ϕ[ν_,x_]:= If[Abs[x]>=Abs[ν],Sqrt[x^2-ν^2] -ν  ArcCos[ν/x]+Pi/4,Pi/4]+Mods[ArcTan[DBesselJ[ν,x],DBesselY[ν,x]]-If[Abs[x]>=Abs[ν],Sqrt[x^2-ν^2] -ν  ArcCos[ν/x]+Pi/4,Pi/4], 2Pi];

Figure

xtib={3,6,9,12}; ytib={-1,-0.5,0,0.5,1};
figCnutauprime=GraphicsColumn[{GraphicsRow[Table[Plot[Cos[Pi τ]DBesselJ[ν,x]+Sin[Pi τ]DBesselY[ν,x]/.τ->{0,1/4,1/2,9/10}//Evaluate, {x,0,10+ν+ν^(1/3)},PlotStyle->mycolors[[4;;7]], AxesLabel->{MaTeX["x"], None},Ticks->{PrTicks[xtib, xtib], PrTicks[ytib, ytib]},Epilog->{Inset[Framed[MaTeX["\\nu="<>ToString[ν]],RoundingRadius->3],{Right,Top},Scaled[{1.5,1.5}]]}], {ν,{0,3}}]],
LineLegend[mycolors[[4;;7]], MaTeX[{"J'_\\nu(x)", "\\mathcal{C}'_{\\nu,\\frac{1}{4}}(x)","Y'_\\nu(x)", "\\mathcal{C}'_{\\nu,\\frac{9}{10}}(x)"}],LegendLayout->"Row"]}]

Export[SaveDir<>"figCnutauprime.pdf", figCnutauprime]

Figure

plotϕjprimez=Table[findAllRoots[DBesselJ[ν,x],{x,0,10}], {ν,{0,2,5}}]; 
PrependTo[plotϕjprimez[[1]],0];
plotϕjprimez
plotϕyprimez=Table[findAllRoots[DBesselY[ν,x],{x,0,10}], {ν,{0,2,5}}]

plotϕ=Plot[Table[ϕ[ν,x],{ν,{0,2,5}}]//Evaluate,{x,0,10}, PlotStyle->mycolors,Ticks->{xti1,yti1}, AxesLabel->MaTeX[{"x","\\varphi_\\nu(x)"}],AxesOrigin->{0,0},Epilog->{
{Thin, Black,Dashed, Table[Line[{{0,yt},{10,yt}}], {yt, 0, 3 Pi,Pi/2}]}, 
Table[
{mycolors[[j]], 
Dashed,
Table[Line[{{plotϕjprimez[[j,k]],0}, {plotϕjprimez[[j,k]], -Pi/2+Pi k}}],{k,1,Length[plotϕjprimez[[j]]]}],
Table[Line[{{plotϕyprimez[[j,k]],0}, {plotϕyprimez[[j,k]], Pi k}}],{k,1,Length[plotϕyprimez[[j]]]}]
},{j,1,3}],
Table[
{mycolors[[j]],
Table[Inset[filledcircle,{plotϕjprimez[[j,k]],0}], {k,1,Length[plotϕjprimez[[j]]]}],
Table[Inset[circle,{plotϕyprimez[[j,k]],0}], {k,1,Length[plotϕyprimez[[j]]]}]
}, {j,1,3}]
}, PlotLegends->MaTeX[{"\\nu=0","\\nu=2","\\nu=5"}]]

Export[SaveDir<>"figphi.pdf", plotϕ]

§1.4. Definitions and properties of the auxiliary functions I

,  its derivative, and asymptotics

θUp[ν_,x_]:=Sqrt[x^2-ν^2]-ν ArcCos[ν/x]-Pi/4;  (* (Subscript[(OverTilde[θ]), ν])(x) *)
DθUp=Simplify[D[θUp[ν,x],x],x>ν>=0] (* (Subscript[(OverTilde[θ]), ν]^('))(x) *)

θUpasympt=Series[θUp[ν,x], {x,Infinity,2}]//Simplify

(Underscript[θ, ~])(x), its derivative, and asymptotics

θDown[ν_,x_]:=Sqrt[x^2-ν^2]-ν ArcCos[ν/x]-Pi/4-(3x^2+2ν^2)/(24(x^2-ν^2)^(3/2));  (* (Subscript[(Underscript[θ, ~]), ν])(x) *)
DθDown=Simplify[D[θDown[ν,x],x],x>ν>=0]; (* (Subscript[(Underscript[θ, ~]), ν]^('))(x) *)
DθDown//TraditionalForm
Simplify[DθDown/.x->Sqrt[ν^2+χ^2],χ>0]

θDownasympt=Series[θDown[ν,x], {x,Infinity,4}]//Simplify

§1.5. Main results I: bounding the phase and zeros of Bessel functions

Inverse functions of bounds with precision

θUpInv[ν_,y_, prec_:MachinePrecision]:=x/.FindRoot[θUp[ν,x]==y,{x,Max[ν,y+1/4Pi (2ν+1)],prec}, WorkingPrecision->prec,Method->"AffineCovariantNewton"];
(* θDownInv[ν_,y_, prec_:MachinePrecision]:=x/.FindRoot[θDown[ν,x]==y,{x,Max[ν,y+1/4Pi (2ν+1)],prec}, WorkingPrecision->prec]//Re; *)
θDownInv[ν_,y_, prec_:MachinePrecision]:=x/.FindRoot[θDown[ν,x]==y,{x,Max[ν,y+1/4Pi (2ν+1)],prec}, WorkingPrecision->prec,Method->"AffineCovariantNewton"];

Figure  bounds

νs={0,5,50,200}; xmaxs= {3,10,58,214}; 
xt={{1,2,3},{6,8,10},{54,58},{207,214}};
xticks=Table[{xt[[j]], MaTeX[xt[[j]], Magnification->0.8]}//Transpose,{j,1,4}];
yticks={{-1/2,0,1/4, 1/2}, MaTeX[{"-\\frac{1}{2}", "0","\\frac{1}{4}", "\\frac{1}{2}"}, Magnification->0.8]}//Transpose;
wp=16;
jzeros=N[BesselJZero[νs,1],wp]
jzerosdown=Table[θUpInv[νs[[j]],Pi/2, wp],{j,1,4}]
jzerosup=Table[θDownInv[νs[[j]],Pi/2, wp],{j,1,4}]
yzeros=N[BesselYZero[νs,1],wp]
yzerosdown=Table[θUpInv[νs[[j]],0, wp],{j,1,4}]
yzerosup=Table[θDownInv[νs[[j]],0, wp],{j,1,4}]
jyzeros= Table[NSolveValues[{(BesselJ[ν,x]-BesselY[ν,x]/.ν->νs[[j]])==0,νs[[j]]<x<xmaxs[[j]]},x, WorkingPrecision->wp][[1]],{j,1,4}]
jyzerosdown=Table[θUpInv[νs[[j]],Pi/4, wp],{j,1,4}]
jyzerosup=Table[θDownInv[νs[[j]],Pi/4, wp],{j,1,4}]
yzeros2=N[BesselYZero[νs,2],wp]
yzerosdown2=Table[θUpInv[νs[[j]],Pi, wp],{j,1,4}]
yzerosup2=Table[θDownInv[νs[[j]],Pi, wp],{j,1,4}]

xmaxs2= {4.1,11,61,217};
xt={{2,4},{6,8,10},{54,58},{208,216}};
xticks=Table[{xt[[j]], MaTeX[xt[[j]], Magnification->0.8]}//Transpose,{j,1,4}];
yticks={yt={-0.5,0,0.25, 0.5,1}, MaTeX[yt, Magnification->0.8],Table[1,5], Table[Directive[Thin,Dashed],5]}//Transpose;
 figcomparisonJ2=GraphicsGrid[ArrayReshape[Table[Plot[1/Pi{θ[ν,x],θUp[ν,x],Max[-Pi/2,θDown[ν,x] ]}/.ν->νs[[j]]//Evaluate, {x,νs[[j]],xmaxs2[[j]]}, PlotLabel->MaTeX["\\nu="<>ToString[νs[[j]]]],Epilog->{Black, Thin,
Line[{{jzeros[[j]],0},{jzeros[[j]],0.75}}],
Line[{{yzeros[[j]],0},{yzeros[[j]],-0.6}}],
Line[{{yzeros2[[j]],0},{yzeros2[[j]],1.15}}],
Line[{{jyzeros[[j]],0.25},{jyzeros[[j]],-0.25}}],
Inset[MaTeX[ToString[DecimalForm[jzerosdown[[j]],6]]<>"<j_{\\nu,1}\\approx"<>ToString[DecimalForm[jzeros[[j]],6]]<>"<"<>ToString[DecimalForm[jzerosup[[j]],6]], Magnification->0.8],{jzeros[[j]],0.75},Scaled[{0.8,0}]],
Inset[MaTeX[ToString[DecimalForm[yzerosdown[[j]],6]]<>"<y_{\\nu,1}\\approx"<>ToString[DecimalForm[yzeros[[j]],6]]<>"<"<>ToString[DecimalForm[yzerosup[[j]],6]], Magnification->0.8],{yzeros[[j]],-0.6},Scaled[{0,1}]],
Inset[MaTeX[ToString[DecimalForm[yzerosdown2[[j]],6]]<>"<y_{\\nu,2}\\approx"<>ToString[DecimalForm[yzeros2[[j]],6]]<>"<"<>ToString[DecimalForm[yzerosup2[[j]],6]], Magnification->0.8],{yzeros2[[j]],1.1},Scaled[{1,0}]],
Inset[MaTeX[ToString[DecimalForm[jyzerosdown[[j]],6]]<>"<c_{\\nu,\\frac{3}{4},1}\\approx"<>ToString[DecimalForm[jyzeros[[j]],6]]<>"<"<>ToString[DecimalForm[jyzerosup[[j]],6]], Magnification->0.8],{jyzeros[[j]],-0.25},Scaled[{0.2,1}]]
}, PlotRange->{-0.75,1.25}, Ticks->{xticks[[j]],yticks}],{j,1,4}], {2,2}],Frame->All]

Export["fig_comparison_J.pdf", figcomparisonJ2]

§1.6. Definitions and properties of the auxiliary functions II

(Subscript[(Underscript[ϕ, ~]), ν])(x), its derivative, and asymptotics

ϕDown[ν_,x_]:=Sqrt[x^2-ν^2]-ν ArcCos[ν/x]+Pi/4;  (* (Subscript[(Underscript[ϕ, ~]), ν])(x) *)
DϕDown=Simplify[D[ϕDown[ν,x],x],x>ν>=0] (* (Subscript[(Underscript[ϕ, ~]), ν]^('))(x) *)

ϕDownasympt=Series[ϕDown[ν,x], {x,Infinity,2}]//Simplify

(Subscript[(OverTilde[ϕ]), ν])(x), its derivative, and asymptotics

ϕUp[ν_,x_]:=Sqrt[x^2-ν^2]-ν ArcCos[ν/x]+Pi/4+(9x^2-2ν^2)/(24(x^2-ν^2)^(3/2));
DϕUp=Simplify[D[ϕUp[ν,x],x],x>ν>=0]

ϕUpasympt=Series[ϕUp[ν,x], {x,Infinity,3}]//Simplify

(Subscript[p, ν]) and

pν[ν_,x_]:=8 (x^6)-8 (ν^6)+4 (x^2) (ν^2) (-1+6 (ν^2))-3 (x^4) (1+8 (ν^2));

xstar[ν_]:=If[ν==0, Sqrt[3/8],Root[-8 (ν^6)+(-4 (ν^2)+24 (ν^4)) ((#1)^2)+(-3-24 (ν^2)) ((#1)^4)+8 ((#1)^6)&,2]];
zstar[ν_]:=ϕUp[ν,xstar[ν]]

Assuming[ ν>0,Series[xstar[ν], {ν,Infinity,1}]]
Assuming[ν>0,Series[ϕUp[ν,ν+(1/4) (7^(1/3)) (ν^(1/3))],{ν,Infinity,1}]]//FullSimplify

zstar[0]/Pi//N
zstarinf=Pi/4+Sqrt[7/18];
zstarinf/Pi//N

ztzt={ztz={0.5, 0.6}, MaTeX[ztz]}//Transpose;
PrependTo[ztzt,{N[zstarinf/Pi], MaTeX["0.4485"],1}];
AppendTo[ztzt,{N[zstar[0]/Pi], MaTeX["0.6398"]}];figzstar=Plot[zstar[ν]/Pi, {ν,0,10},PlotRange->{Full,Full}, AxesLabel->MaTeX[{"\\nu","\\frac{z^\\star_\\nu}{\\pi}"}],AxesOrigin->{0,zstarinf/Pi-0.02},Ticks->{{xtz={3,6,9}, MaTeX[xtz]}//Transpose, ztzt},AspectRatio->1/2]

Export[SaveDir<>"figzstar.pdf", figzstar]

§1.7. Main results II: bounding zeros of derivatives of Bessel functions

Figure “Forbidden region”

zstar[0]/Pi-1/2//N
FindRoot[zstar[ν]/Pi==1/2,{ν,1}]
figforbid=RegionPlot [τ<=zstar[ν]/Pi-1/2,{ν,0,1.4},{τ,0,0.15},Frame->Off,Axes->True, AxesLabel->MaTeX[{"\\nu","\\tau"}],AspectRatio->1/2,PlotStyle->Directive[Opacity[0.25],mycolors[[4]]],BoundaryStyle->mycolors[[4]],Ticks->{{xtr={1,1.19876},MaTeX[xtr]}//Transpose,{ytr={0.1,0.1398},MaTeX[ytr]}//Transpose}]

Export[SaveDir<>"figforbid.pdf", figforbid]

Inverse functions of bounds with precision

ϕUpInv[ν_,y_, prec_:MachinePrecision]:=x/.FindRoot[ϕUp[ν,x]==y,{x,Max[ν,y+1/4Pi (2ν+1)],prec}, WorkingPrecision->prec,Method->"AffineCovariantNewton"];
ϕDownInv[ν_,y_, prec_:MachinePrecision]:=x/.FindRoot[ϕDown[ν,x]==y,{x,Max[ν,y+1/4Pi (2ν+1)],prec}, WorkingPrecision->prec,Method->"AffineCovariantNewton"];
wp=8;

Figure  bounds

xmaxs2= {4.1,11,61,217};
jpzeros=Table[NSolveValues[{(DBesselJ[ν,x]/.ν->νs[[j]])==0,νs[[j]]<x<xmaxs2[[j]]},x, WorkingPrecision->wp],{j,1,4}];
PrependTo[jpzeros[[1]],0];
jpzeros
jpzerosdown=Table[If[j==1&&k==0,,ϕUpInv[νs[[j]],Pi/2+Pi k, wp]],{j,1,4}, {k,0,1}]
jpzerosup=Table[If[j==1&&k==0,,ϕDownInv[νs[[j]],Pi/2+Pi k, wp]],{j,1,4}, {k,0,1}]
ypzeros=Table[If[j>1,NSolveValues[{(DBesselY[ν,x]/.ν->νs[[j]])==0,νs[[j]]<x<xmaxs2[[j]]},x, WorkingPrecision->wp][[1]],N[BesselYZero[1,1],wp]],{j,1,4}]
ypzerosdown=Table[ϕUpInv[νs[[j]],Pi, wp],{j,1,4}]
ypzerosup=Table[ϕDownInv[νs[[j]],Pi, wp],{j,1,4}]
jypzeros=Table[NSolveValues[{(DBesselJ[ν,x]-DBesselY[ν,x]/.ν->νs[[j]])==0,Floor[ypzeros[[j]]]<x<Ceiling[jpzeros[[j]][[2]]]},x, WorkingPrecision->wp][[1]],{j,1,4}]
jypzerosdown=Table[NSolveValues[{ϕUp[νs[[j]],x]==5Pi/4,Floor[ypzeros[[j]]]<x<Ceiling[jpzeros[[j]][[2]]]},x, WorkingPrecision->wp][[1]],{j,1,4}]
jypzerosup=Table[NSolveValues[{ϕDown[νs[[j]],x]==5Pi/4,Floor[ypzeros[[j]]]<x<Ceiling[jpzeros[[j]][[2]]]},x, WorkingPrecision->wp][[1]],{j,1,4}]

yticks={yt={0,0.5, 1,1.25,1.5}, MaTeX[yt, Magnification->0.8],Table[1,5], Table[Directive[Thin,Dashed],5]}//Transpose;
xt={{2,4},{7,10},{54,58},{208,216}};
xticks=Table[{xt[[j]], MaTeX[xt[[j]], Magnification->0.8]}//Transpose,{j,1,4}];
Do[PrependTo[xticks[[j]],{xstar[νs[[j]]],MaTeX["x^\\star_\\nu", Magnification->0.8],0}],{j,1,4}]
figcomparisonJp=GraphicsGrid[ArrayReshape[Table[
Show[
Plot[1/Pi{ϕ[ν,x],ϕUp[ν,x],ϕDown[ν,x]}/.ν->νs[[j]]//Evaluate, {x,xstar[νs[[j]]],xmaxs2[[j]]}, PlotLabel->MaTeX["\\nu="<>ToString[νs[[j]]]], Epilog->{Black, Thin,
Line[{{jpzeros[[j]][[2]],0},{jpzeros[[j]][[2]],1.65}}],
Line[{{jypzeros[[j]],1.25},{jypzeros[[j]],-0.45}}],
Line[{{ypzeros[[j]],0},{ypzeros[[j]],1.3}}],
If[j>1,Line[{{jpzeros[[j]][[1]],0.5},{jpzeros[[j]][[1]],-0.25}}]],
Line[{{xstar[νs[[j]]],0},{xstar[νs[[j]]],zstar[νs[[j]]]/Pi}}],
Inset[MaTeX[ToString[DecimalForm[jpzerosdown[[j]][[2]],6]]<>"<j'_{\\nu,2}\\approx"<>ToString[DecimalForm[jpzeros[[j]][[2]],6]]<>"<"<>ToString[DecimalForm[jpzerosup[[j]][[2]],6]], Magnification->0.8],{Right,1.65},Scaled[{1,0}]],
If[j>1,Inset[MaTeX[ToString[DecimalForm[jpzerosdown[[j]][[1]],6]]<>"<j'_{\\nu,1}\\approx"<>ToString[DecimalForm[jpzeros[[j]][[1]],6]]<>"<"<>ToString[DecimalForm[jpzerosup[[j]][[1]],6]], Magnification->0.8],{νs[[j]],-0.25},Scaled[{0,1}]]],
Inset[MaTeX[ToString[DecimalForm[ypzerosdown[[j]],6]]<>"<y'_{\\nu,1}\\approx"<>ToString[DecimalForm[ypzeros[[j]],6]]<>"<"<>ToString[DecimalForm[ypzerosup[[j]],6]], Magnification->0.8],{ypzeros[[j]],1.3},Scaled[{0.8,0}]],
Inset[MaTeX[ToString[DecimalForm[jypzerosdown[[j]],6]]<>"<c'_{\\nu,\\frac{3}{4},1}\\approx"<>ToString[DecimalForm[jypzeros[[j]],6]]<>"<"<>ToString[DecimalForm[jypzerosup[[j]],6]], Magnification->0.8],{Right,-0.45},Scaled[{1,1}]]
}, PlotRange->{{νs[[j]],xmaxs2[[j]]},{-0.7,1.85}}, Ticks->{xticks[[j]],yticks}],
Plot[1/Pi{ϕ[ν,x],zstar[ν],ϕDown[ν,x]}/.ν->νs[[j]]//Evaluate, {x,νs[[j]],xstar[νs[[j]]]}, PlotRange->{Automatic,{-0.7,1.85}}]
],{j,1,4}], {2,2}],Frame->All]

Export[SaveDir<>"fig_comparison_Jp.pdf", figcomparisonJp]

§2. Proofs

§2.1. Liouville’s transform and phase functions

ℱ[f_,t_,x_]:=Cos[f-Pi t]/Sqrt[D[f,x]];
𝒱[f_,x_]:=D[f,x]^2+1/2 D[f, {x,3}]/D[f,x]-3/4(D[f, {x,2}]/D[f,x])^2//Simplify
CheckSchrö[𝒜_,𝒫_,x_]:=D[𝒜,x,x]+𝒫 𝒜 //Simplify; (* checks Schrödinger equation for function 𝒜 and potential 𝒫 *)

Lemma 2.1

CheckSchrö[ℱ[f[x],t,x], 𝒱[f[x],x],x]

Lemma 2.3

Simplify[D[ArcTan[𝒜2[x]/𝒜1[x]],x]]/.Wronskian[{𝒜1[x],𝒜2[x]},x]->𝒲 (*𝒜1[x] (𝒜2)'[x]-𝒜2[x] (𝒜1)'[x]->𝒲*)

Simplify[FullSimplify[ℱ[ArcTan[𝒜2[x]/𝒜1[x]],t,x]//TrigExpand, 𝒜1[x]>0 && 𝒲>0]/.Wronskian[{𝒜1[x],𝒜2[x]},x]->𝒲,((𝒜1[x])^2)+((𝒜2[x])^2)>0]

(* extra check *)
Simplify[𝒱[ArcTan[𝒜2[x]/𝒜1[x]],x]/. {𝒜1''[x]->-𝒫[x] 𝒜1[x],𝒜1'''[x]->-𝒫'[x] 𝒜1[x]-𝒫[x] 𝒜1'[x],𝒜2''[x]->-𝒫[x] 𝒜2[x],𝒜2'''[x]->-𝒫'[x] 𝒜2[x]-𝒫[x] 𝒜2'[x]}]

Lemma 2.4

CheckSchrö[Sqrt[x]BesselJ[ν,x],1-(ν^2-1/4)/x^2,x]//FullSimplify

𝒱θ[ν_,x_]:=1-(ν^2-1/4)/x^2;

CheckSchrö[x^(3/2)/Sqrt[x^2-ν^2]DBesselJ[ν,x],1-(ν^2-1/4)/x^2-(2ν^2+x^2)/(x^2-ν^2)^2,x]//FullSimplify

𝒱ϕ[ν_,x_]:=1-(ν^2-1/4)/x^2-(2ν^2+x^2)/(x^2-ν^2)^2;

§2.3. Proof of Theorem 1.4

Upper bound

Map[Collect[#,x]&,Simplify[𝒱[θUp[ν,x],x], x>ν>=0]]

𝒱θUp[ν_,x_]:=((4 (x^6)-12 (x^4) (ν^2)+(ν^4)-4 (ν^6)+6 (x^2) (ν^2) (-1+2 (ν^2)))/(4 ((((x^3)-x (ν^2)))^2)));

Condition ((Subscript[C, 2]))

Simplify[𝒱θ[ν,x]-𝒱θUp[ν,x], x>ν>=0]

Condition ((Subscript[C, 3]^(')))

θUpasympt-θasympt//Simplify

Lower bound

{q1,q2}=Simplify[𝒱[θDown[ν,x],x], x>ν>=0]//Together//NumeratorDenominator;

Collect[q1, x,Simplify]//TraditionalForm

𝒱θDown[ν_,x_]:=((x^16)+32 (x^18)+1920 (x^20)+2048 (x^22)+4096 (x^24)+16 (x^14) (ν^2)+704 (x^16) (ν^2)+18688 (x^18) (ν^2)-16384 (x^20) (ν^2)-49152 (x^22) (ν^2)+96 (x^12) (ν^4)-3568 (x^14) (ν^4)-117376 (x^16) (ν^4)+56320 (x^18) (ν^4)+270336 (x^20) (ν^4)+256 (x^10) (ν^6)+9360 (x^12) (ν^6)+209408 (x^14) (ν^6)-107520 (x^16) (ν^6)-901120 (x^18) (ν^6)+256 (x^8) (ν^8)-14960 (x^10) (ν^8)-80000 (x^12) (ν^8)+122880 (x^14) (ν^8)+2027520 (x^16) (ν^8)+13168 (x^8) (ν^10)-183040 (x^10) (ν^10)-86016 (x^12) (ν^10)-3244032 (x^14) (ν^10)-5504 (x^6) (ν^12)+244352 (x^8) (ν^12)+43008 (x^10) (ν^12)+3784704 (x^12) (ν^12)+768 (x^4) (ν^14)-103424 (x^6) (ν^14)-30720 (x^8) (ν^14)-3244032 (x^10) (ν^14)+4352 (x^4) (ν^16)+30720 (x^6) (ν^16)+2027520 (x^8) (ν^16)+5120 (x^2) (ν^18)-20480 (x^4) (ν^18)-901120 (x^6) (ν^18)+7168 (x^2) (ν^20)+270336 (x^4) (ν^20)-1024 (ν^22)-49152 (x^2) (ν^22)+4096 (ν^24))/(64 (x^2) ((((x^2)-(ν^2)))^5) ((((x^4)+8 (x^6)+4 (x^2) (ν^2)-24 (x^4) (ν^2)+24 (x^2) (ν^4)-8 (ν^6)))^2));

Condition  ((Subscript[C, 2]))

Collect[FullSimplify[(64 (χ^5) (((5 (ν^4)+6 (ν^2) χ+(χ^2)+8 (χ^3)))^2))(𝒱θDown[ν,Sqrt[ν^2+χ]]-𝒱θ[ν,Sqrt[ν^2+χ]])],χ, Simplify]//TraditionalForm

Condition ((Subscript[C, 3]^(')))

θasympt-θDownasympt//Simplify

§2.4 Proof of Lemma 1.8

DϕUpN=Collect[Simplify[pν[ν,Sqrt[ν^2+ξ]],ξ>0], ξ, Simplify]
DerϕUpN=D[DϕUpN,ξ]
Dχ0=SolveValues[DerϕUpN==0&&ξ>0,  ξ,Assumptions->ν>=0]//Apart

dnu=FullSimplify[((D[ϕUp[μ,x],μ]//Simplify))/.x-> κ μ, κ>1 && μ>=0]

μsn=Simplify[SolveValues[dnu==0 , μ,Assumptions->κ>1][[1]], κ>1 ]

κeq=Simplify[pν[μ, κ μ]/.μ->μsn, κ>1 ]

Simplify[SolveValues[(κeq/.ArcSec[κ]->A)==0,A],κ>1 ][[1]]-ArcSec[κ]

K[κ_]:=((2-25 (κ^2)+23 (κ^4))/(3 (κ^2) (Sqrt[-1+(κ^2)]) (4+3 (κ^2))))-ArcSec[κ];
Simplify[K'[κ], κ>1 ]

§. Proof of Theorem 1.9

Lower bound

𝒱ϕDown[ν_,x_]:=𝒱θUp[ν,x];

Condition ((Subscript[C, 2]))

Simplify[𝒱ϕDown[ν,x]-𝒱ϕ[ν,x], x>ν>=0]

Condition ((Subscript[C, 3]^(')))

ϕasympt-ϕDownasympt//Simplify

Upper bound

{q3,q4}=Simplify[𝒱[ϕUp[ν,x],x], x>ν>=0]//Together//NumeratorDenominator

Collect[q3, x,Simplify]//TraditionalForm

𝒱ϕUp[ν_,x_]:=(81 (x^16)-864 (x^18)-1152 (x^20)-6144 (x^22)+4096 (x^24)+432 (x^14) (ν^2)+1344 (x^16) (ν^2)-24320 (x^18) (ν^2)+40960 (x^20) (ν^2)-49152 (x^22) (ν^2)+864 (x^12) (ν^4)-9200 (x^14) (ν^4)+138624 (x^16) (ν^4)-91136 (x^18) (ν^4)+270336 (x^20) (ν^4)+768 (x^10) (ν^6)+26512 (x^12) (ν^6)-261632 (x^14) (ν^6)-9216 (x^16) (ν^6)-901120 (x^18) (ν^6)+256 (x^8) (ν^8)-25584 (x^10) (ν^8)+165760 (x^12) (ν^8)+466944 (x^14) (ν^8)+2027520 (x^16) (ν^8)+7664 (x^8) (ν^10)+91392 (x^10) (ν^10)-1118208 (x^12) (ν^10)-3244032 (x^14) (ν^10)-640 (x^6) (ν^12)-183680 (x^8) (ν^12)+1419264 (x^10) (ν^12)+3784704 (x^12) (ν^12)+768 (x^4) (ν^14)+80896 (x^6) (ν^14)-1112064 (x^8) (ν^14)-3244032 (x^10) (ν^14)-768 (x^4) (ν^16)+546816 (x^6) (ν^16)+2027520 (x^8) (ν^16)-5120 (x^2) (ν^18)-159744 (x^4) (ν^18)-901120 (x^6) (ν^18)+23552 (x^2) (ν^20)+270336 (x^4) (ν^20)-1024 (ν^22)-49152 (x^2) (ν^22)+4096 (ν^24))/(64 (x^2) ((((x^2)-(ν^2)))^5) (((-3 (x^4)+8 (x^6)-4 (x^2) (ν^2)-24 (x^4) (ν^2)+24 (x^2) (ν^4)-8 (ν^6)))^2));

Condition ((Subscript[C, 2]))

Simplify[4096 (x^2) ((((x^2)-(ν^2)))^10)D[ϕUp[ν,x],x]^2(𝒱ϕ[ν,x]-𝒱ϕUp[ν,x]), x>ν>=0]

δ𝒱ϕ[ν_,x_]:=4032 (x^18)+4096 (ν^18)+1024 (x^2) (ν^14) (-1+17 (ν^2))-3584 (x^4) (ν^12) (-1+41 (ν^2))+16 (x^16) (27+1208 (ν^2))-64 (x^8) (ν^6) (12-337 (ν^2)+2338 (ν^4))+64 (x^6) (ν^8) (-4-164 (ν^2)+4543 (ν^4))-16 (x^10) (ν^4) (54+1415 (ν^2)+12796 (ν^4))+16 (x^12) (ν^2) (-27+701 (ν^2)+20272 (ν^4))-(x^14) (81+2640 (ν^2)+158656 (ν^4));

(δ𝒱ϕ[0,x]x^(-14)//Simplify)/.x->Sqrt[3/8+ξ]//Simplify

δ𝒱ϕ1= Collect[ν^(-14) δ𝒱ϕ[ν,ν+ν ζ],ν, Simplify]

δineq1=δ𝒱ϕ1/.ν^2->0

Collect[(ν^(-4)) pν[ν,ν(ζ+1)]//Simplify, ν]
νineq=FullSimplify[Reduce[{(ν^(-4)) pν[ν,ν(ζ+1)]>=0,ζ>0,ν>0}, ν], ζ>0]
Map[FullSimplify[#^4, ζ>0]&, νineq]

Collect[δineq1/.ν->νineq[[2]], ζ, Simplify]//FullSimplify

Condition ((Subscript[C, 3]^(')))

ϕUpasympt-ϕasympt//Simplify

§3. Derivatives of ultraspherical Bessel functions

§3.1. Setup III

Definitions

D[x^(-η)BesselJ[ν,x],x]//FullSimplify

Uprime[ν_,η_, x_]:=(x^(-η)) (BesselJ[-1+ν,x]-(((η+ν) BesselJ[ν,x])/x))
Wprime[ν_,η_, x_]:=(x^(-η)) (BesselY[-1+ν,x]-(((η+ν) BesselY[ν,x])/x))

Figure  and

xtib={5,10}; ytib={-10,-5,0,5,10};
figUprimeWprime=GraphicsColumn[{GraphicsRow[Table[Plot[{x^2 Uprime[ν,1,x],x^11 Uprime[ν,10,x],x^2 Wprime[ν,1,x],x^11 Wprime[ν,10,x]} //Evaluate, {x,0,10+ν+ν^(1/3)},PlotStyle->mycolors[[4;;7]], AxesLabel->{MaTeX["x"], None},Ticks->{PrTicks[xtib, xtib], PrTicks[ytib, ytib]},Epilog->{Inset[Framed[MaTeX["\\nu="<>ToString[ν]],RoundingRadius->3],{Right, If[ν==0,Bottom,Top]},Scaled[{1.5,If[ν==0,-0.25,1.25]}]]}], {ν,{0,3}}]],
LineLegend[mycolors[[4;;7]], MaTeX[{"x^2 U'_{\\nu,1}(x)", "x^{11} U'_{\\nu,10}(x)","x^2 W'_{\\nu,1}(x)", "x^{11} W'_{\\nu,10}(x)"}],LegendLayout->"Row"]}]

Export[SaveDir<>"figUWprime.pdf", figUprimeWprime]

§3.2. Phase function of ultraspherical Bessel derivatives

Lemma 3.1

μ2[ν_,η_]:=ν^2-η^2;
𝒱ψ[ν_,η_,x_]:=1-(ν^2-1/4)/x^2+2(1-η)/(x^2-μ2[ν,η])-3x^2/(x^2-μ2[ν,η])^2;

CheckSchrö[x^(η+3/2)/Sqrt[x^2-μ2[ν,η]]Uprime[ν,η,x],𝒱ψ[ν,η,x], x]//FullSimplify
CheckSchrö[x^(η+3/2)/Sqrt[x^2-μ2[ν,η]]Wprime[ν,η,x],𝒱ψ[ν,η,x], x]//FullSimplify

Wronskian[{x^(η+3/2)/Sqrt[x^2-μ2[ν,η]]Uprime[ν,η,x],x^(η+3/2)/Sqrt[x^2-μ2[ν,η]]Wprime[ν,η,x]},x]

D[ArcTan[Wprime[ν,η,x]/Uprime[ν,η,x]], x]-2(x^2-μ2[ν,η])/(Pi x^(2η+3)(Uprime[ν,η,x]^2+Wprime[ν,η,x]^2))//FullSimplify

Lemma 3.2

seriesLsquared=Collect[Asymptotic[x^(2η)(Uprime[ν,η, x]^2+Wprime[ν,η, x]^2),{x,Infinity,5}]//Simplify,x]

ser=Series[2/Pi (x^2-ν^2+η^2)/x^3 / seriesLsquared,{x,Infinity,4}]//Normal;
seras=Collect[Integrate[ser,x]+C,x,Collect[Numerator[#],{ν,η},Simplify]/Denominator[#]&]

x^(-η)Asymptotic[Sqrt[seriesLsquared] Cos[seras], {x,Infinity,1}]
Asymptotic[Uprime[ν,η,x], {x,Infinity,1}]

ψasympt=seras/.C->-(1/4) π (-1+2 ν)

§3.3. Definitions and properties of the auxiliary functions III
§3.4. Main results III: bounding the phase and zeros of derivatives of ultraspherical Bessel functions

(Subscript[(Underscript[ψ, ~]), ν,η])(x) and  its derivative

ψDown[ν_,η_,x_]:=Sqrt[x^2-μ^2]+η/Sqrt[x^2-μ^2]+Pi/4 (η^2/μ+2(μ-ν)+1)-(μ+η^2/(2μ))ArcCos[μ/x];
Collect[Limit[ψDown[ν,η,x]/.μ->Sqrt[μ2[ν,η]],η->ν,Direction->"FromBelow",Assumptions->x>0],x]
ψDown1[ν_,η_,x_]:=If[ν>η,Sqrt[x^2-ν^2+η^2]+η/Sqrt[x^2-ν^2+η^2]+Pi/4 (η^2/Sqrt[ν^2-η^2]+2(Sqrt[ν^2-η^2]-ν)+1)-(Sqrt[ν^2-η^2]+η^2/(2Sqrt[ν^2-η^2]))ArcCos[Sqrt[ν^2-η^2]/x],
x+((ν (2+ν))/(2 x))+(1/4) (π-2 π ν)
];

DψDown=FullSimplify[D[ψDown[ν,η,x],x]//Together,x>μ>0]

(Subscript[x, μ,η]^(#))

FullSimplify[SolveValues[Numerator[DψDown]==0,x][[4]], η>0]

xhash[ν_,η_]:=(1/2) (Sqrt[η (2+η)+4 μ2[ν,η]+(Sqrt[η (η (((2+η))^2)+16 μ2[ν,η])])])
xhashμ[μ_,η_]:=(1/2) (Sqrt[η (2+η)+4 μ^2+(Sqrt[η (η (((2+η))^2)+16 μ^2)])])

(Subscript[𝒱, (Subscript[(Underscript[ψ, ~]), ν,η])])

𝒱ψDown=𝒱[ψDown[ν,η,x],x];
{q5,q6}=Simplify[𝒱ψDown//Together//NumeratorDenominator,x>μ>0];
(* q5=Collect[q5, x, Collect[#,η,Collect[#,μ,Simplify]&]&] *)
q5=Collect[q5, {x},Simplify];
q6=FullSimplify[Expand[q6],x>0];
q5//TraditionalForm
q6//TraditionalForm

𝒱ψ[Sqrt[μ^2+η^2],η,x]

(Subscript[r, μ,η])

rμη=Collect[Together[Simplify[16 x^4(x^2-μ^2)^6 DψDown^2(𝒱ψDown-𝒱ψ[Sqrt[μ^2+η^2],η,x]), x>μ>0]],{x,μ}, Simplify]

rμη//TraditionalForm

rμη/.μ->0//Simplify

(x^(&))

xamp[ν_,η_]:=NSolveValues[(12 (x^14)-72 (x^12) η+12 (x^10) (η^2)-20 (x^12) (η^2)-20 (x^10) (η^3)+16 (x^12) (η^3)+16 (x^8) (η^4)-45 (x^10) (η^4)+4 (x^12) (η^4)+32 (x^8) (η^5)-24 (x^10) (η^5)+24 (x^8) (η^6)-4 (x^10) (η^6)+8 (x^8) (η^7)+(x^8) (η^8)-44 (x^12) (μ^2)+40 (x^10) η (μ^2)+52 (x^8) (η^2) (μ^2)+48 (x^10) (η^2) (μ^2)+96 (x^8) (η^3) (μ^2)-80 (x^10) (η^3) (μ^2)+147 (x^8) (η^4) (μ^2)-24 (x^10) (η^4) (μ^2)-32 (x^6) (η^5) (μ^2)+96 (x^8) (η^5) (μ^2)-48 (x^6) (η^6) (μ^2)+20 (x^8) (η^6) (μ^2)-24 (x^6) (η^7) (μ^2)-4 (x^6) (η^8) (μ^2)+40 (x^10) (μ^4)+280 (x^8) η (μ^4)-80 (x^6) (η^2) (μ^4)+8 (x^8) (η^2) (μ^4)-116 (x^6) (η^3) (μ^4)+160 (x^8) (η^3) (μ^4)-171 (x^6) (η^4) (μ^4)+60 (x^8) (η^4) (μ^4)-144 (x^6) (η^5) (μ^4)+24 (x^4) (η^6) (μ^4)-40 (x^6) (η^6) (μ^4)+24 (x^4) (η^7) (μ^4)+6 (x^4) (η^8) (μ^4)+40 (x^8) (μ^6)-360 (x^6) η (μ^6)+16 (x^4) (η^2) (μ^6)-112 (x^6) (η^2) (μ^6)+24 (x^4) (η^3) (μ^6)-160 (x^6) (η^3) (μ^6)+81 (x^4) (η^4) (μ^6)-80 (x^6) (η^4) (μ^6)+96 (x^4) (η^5) (μ^6)+40 (x^4) (η^6) (μ^6)-8 (x^2) (η^7) (μ^6)-4 (x^2) (η^8) (μ^6)-100 (x^6) (μ^8)+80 (x^4) η (μ^8)+108 (x^4) (η^2) (μ^8)+16 (x^2) (η^3) (μ^8)+80 (x^4) (η^3) (μ^8)-12 (x^2) (η^4) (μ^8)+60 (x^4) (η^4) (μ^8)-24 (x^2) (η^5) (μ^8)-20 (x^2) (η^6) (μ^8)+(η^8) (μ^8)+68 (x^4) (μ^10)+32 (x^2) η (μ^10)-32 (x^2) (η^2) (μ^10)-16 (x^2) (η^3) (μ^10)-24 (x^2) (η^4) (μ^10)+4 (η^6) (μ^10)-16 (x^2) (μ^12)+4 (η^4) (μ^12)/.μ->Sqrt[μ2[ν,η]])==0,x,Reals]//Max
zamp[ν_,η_]:=ψDown1[ν,η,xamp[ν,η]];

Lemma 3.4

FullSimplify[((rμη/.{x->xhashμ[μ,η]}))]/.(Sqrt[η (η (((2+η))^2)+16 (μ^2))])->ρ

Asymptotics of (Subscript[(Underscript[ψ, ~]), ν,η])(x)

ψDownasympt=Series[ψDown[ν,η,x],{x,Infinity,3}]/.μ->Sqrt[ν^2-η^2]//Simplify

Condition  ((Subscript[C, 3]^(')))

ψasympt-ψDownasympt/.μ->Sqrt[μ2[ν,η]]//Simplify

Computing  following the ideas from [Ho17]. NB: may not always work for some values of the parameters beyond those in the paper

ψ[ν_,η_,x_]:=Module[{aaa},aaa=ArcTan[BesselJ[-1+ν,x]-(((η+ν) BesselJ[ν,x])/x),BesselY[-1+ν,x]-(((η+ν) BesselY[ν,x])/x)];aaa-2 π Round[(1/(2 π))(aaa-If[x≥xhash[ν,η],Max[ψDown1[ν,η, x],0],(π/2)])] 
];

Inverse function of the bound

ψDownInv[ν_,η_,y_, prec_:MachinePrecision]:=If[y>zamp[ν,η],SolveValues[ψDown1[ν,η,x]==y && x>xamp[ν,η],x, WorkingPrecision->prec][[1]],Null];

Figure  bounds

νs3={5,5,50,200};
ηs3={1,5,25,80};
xmaxs3= {11,11,61,217};
xamp3=Table[xamp[νs3[[j]],ηs3[[j]]],{j,1,4}]
zamp3=Table[zamp[νs3[[j]],ηs3[[j]]],{j,1,4}]
μs3=Sqrt[μ2[νs3,ηs3]]//N
wp=8;

yticks={yt={0.5, 1,1.5}, MaTeX[yt, Magnification->0.8],Table[1,3], Table[Directive[Thin,Dashed],3]}//Transpose;
xt={{6,10},{4,8},{48,56},{190,210}};
xticks=Table[{xt[[j]], MaTeX[xt[[j]], Magnification->0.8]}//Transpose,{j,1,4}];
Do[PrependTo[xticks[[j]],{xamp3[[j]],MaTeX["x^@_{\\nu,\\eta}", Magnification->0.8],0}],{j,1,4}]

upzeros=Table[NSolveValues[{BesselJ[-1+νs3[[j]],x]-(((ηs3[[j]]+νs3[[j]]) BesselJ[νs3[[j]],x])/x)==0,μs3[[j]]<x<xmaxs3[[j]]},x, WorkingPrecision->wp],{j,1,4}];
PrependTo[upzeros[[2]],μs3[[2]]];
upzeros
wpzeros=Table[NSolveValues[{BesselY[-1+νs3[[j]],x]-(((ηs3[[j]]+νs3[[j]]) BesselY[νs3[[j]],x])/x)==0,μs3[[j]]<x<xmaxs3[[j]]},x, WorkingPrecision->wp],{j,1,4}]
upzerosup=Table[ψDownInv[νs3[[j]],ηs3[[j]],Pi/2+Pi k],{j,1,4},{k,0,1}]
wpzerosup=Table[ψDownInv[νs3[[j]],ηs3[[j]],Pi],{j,1,4}]

figcomparisonUp=GraphicsGrid[ArrayReshape[Table[
Show[
Plot[1/Pi{ψ[νs3[[j]],ηs3[[j]],x],ψDown1[νs3[[j]],ηs3[[j]],x]}//Evaluate, {x,xamp3[[j]],xmaxs3[[j]]}, PlotLabel->MaTeX["\\nu="<>ToString[νs3[[j]]]<>", \\eta="<>ToString[ηs3[[j]]]], Epilog->{Black, Thin,
Line[{{upzeros[[j]][[2]],1.5},{upzeros[[j]][[2]],1.65}}],
Line[{{wpzeros[[j]][[1]],1.0},{wpzeros[[j]][[1]],1.25}}],
Line[{{upzeros[[j]][[1]],0.5},{upzeros[[j]][[1]],-0.5}}],
Line[{{xamp3[[j]],0},{xamp3[[j]],zamp3[[j]]/Pi}}],
Inset[MaTeX["u'_{\\nu,\\eta,2}\\approx"<>ToString[DecimalForm[upzeros[[j]][[2]],6]]<>"<"<>ToString[DecimalForm[upzerosup[[j]][[2]],6]], Magnification->0.8],{upzeros[[j]][[2]],1.65},Scaled[{1,0}]],
Inset[MaTeX["w'_{\\nu,\\eta,1}\\approx"<>ToString[DecimalForm[wpzeros[[j]][[1]],6]]<>"<"<>ToString[DecimalForm[wpzerosup[[j]],6]], Magnification->0.8],{wpzeros[[j]][[1]],1.25},Scaled[{1,0}]],
Inset[If[j!=2,
MaTeX["u'_{\\nu,\\eta,1}\\approx"<>ToString[DecimalForm[upzeros[[j]][[1]],6]]<>If[j<3,
 "\\text{ no bound}",
"<"<>ToString[DecimalForm[upzerosup[[j]][[1]],6]]], Magnification->0.8],
MaTeX["u'_{\\nu,\\eta,1}=0", Magnification->0.8]
],{upzeros[[j]][[1]],-0.5},Scaled[{0,1}]]
},PlotRange->{{μs3[[j]],xmaxs3[[j]]},{-0.85,1.85}}, Ticks->{xticks[[j]],yticks}], 
Plot[1/Pi{ψ[νs3[[j]],ηs3[[j]],x],ψDown1[νs3[[j]],ηs3[[j]],x]}//Evaluate, {x,μs3[[j]],xamp3[[j]]}, PlotStyle->{Automatic,Dashed}, PlotRange->{Automatic,{-0.7,1.85}}]
],{j,1,4}], {2,2}],Frame->All]

Export[SaveDir<>"figcomparisonUp.pdf",figcomparisonUp]

§4. Benchmarking and conclusions

Series coefficients

A1[ν_,β_]:=β; (*  *)
A2[ν_,β_]:=β-(4ν^2-1)/(8β); (*  *)
A3[ν_,β_]:=β-(4ν^2-1)/(8β)-4(4ν^2-1)(28ν^2-31)/(3(8β)^3); (*  *)
β[ν_, k_]:=Pi(ν/2+k-1/4); (*  *)
(* A4[ν_,β_]:=β-(4ν^2-1)/(8β)-4(4ν^2-1)(28ν^2-31)/(3(8β)^3)-32(4ν^2-1)(1328ν^43928ν^2+3779)/(15(8β)^5); *)

Series[A1[ν,β[ν,k]]+A2[ν,β[ν,k]]+A3[ν,β[ν,k]], {k,Infinity,4}]

Previous bounds

(Subscript[(OverTilde[q]), ν,k])

QWUp[ν_,k_,  prec_:MachinePrecision]:=Module[{az=AiryAiZero[k]},If[ν>0,N[ν-az (ν/2)^(1/3)+3/20 az^2(2/ν)^(1/3),wp],None]];

(Subscript[(Underscript[q, ~]), ν,k])

QWDown[ν_,k_, prec_:MachinePrecision]:=Module[{az=AiryAiZero[k]},
If[ν>0,N[ν-az (ν/2)^(1/3),wp],None]];

(Subscript[(OverTilde[ℓ]), ν,k])

ELUp[ν_,k_]:=If[ν<=1/2,
A2[ν,β[ν,k]],
A3[ν,β[ν,k]]
];

(Subscript[(Underscript[ℓ, ~]), ν,k])

ELDown[ν_,k_]:=Which[ν<=1/2,A3[ν,β[ν,k]], 
ν<Sqrt[31/28], A2[ν,β[ν,k]],
True, None]

(Subscript[(OverTilde[ℓ]), ν,k])

ELprimeUp[ν_,k_,  prec_:MachinePrecision]:=Module[{az=AiryAiZero[k], az1=If[k==1,0, AiryAiZero[k-1]], apz,z},
apz=Abs[z/.FindRoot[AiryAiPrime[z],{z,az,az1}, WorkingPrecision->wp]];
N[apz (((((8 (apz^(3/2)))/27)+(ν/2)))^(1/3))+ν+((9 (apz^2))/(10 (2^(2/3)) (((16 (apz^(3/2))+27 ν))^(1/3)))),wp]];

Appendix B

§B.1. Bounds for zeros of Bessel functions

Parameters

νsB1= {0,1/2, 1, 5,10,50,100, 1000, 10000,100000,500000};
ksB1={1,2,5,10,50,100, 1000, 10000,100000,500000};
wp=16;

Pretty LATEX formatting

MyMax[args_]:=Select[args,NumberQ]//Max;
MyMin[args_]:=Select[args,NumberQ]//Min;

MyDecForm[x_,l_]:=Module[{lint},If[NumberQ[x],lint=StringLength[ToString[Floor[x]]]; ToString[DecimalForm[x,{l,l-lint} ]]," "]];

endst="\\\\\n";
endeverystr="\\\\\\nopagebreak\n";

MakeLineJ[k_, bd_, l_, Myf_,showdiff_:False, dl_:1]:=Module[{str,minmax, b, stradd},
str= " & " <>ToString[k];
minmax=Myf[bd];
Do[
stradd= If[TrueQ[b==minmax],
 " & {\\color{red}"<> MyDecForm[b,l]<>"}" ,
 " & "<> MyDecForm[b,l]];
str=str<> stradd,{b,bd}];
If[showdiff, 
If[bd[[1]]!=minmax,str=str<> " & "<> ToString[TeXForm[ScientificForm[bd[[1]]/minmax-1,dl,NumberPoint->If[dl==1,"","." ]]]],str=str <>" & "]
 ];
str
];

Table 1 (remove semicolon to get LATEX code)

StringJoin[Table[StringJoin["\\multirow{"<>ToString[Length[ksB1]]<>"}{*}{$",
ToString[TeXForm[ν]],
"$}",
Table[StringJoin[MakeLineJ[k, {θDownInv[ν,Pi (k -1/2), wp],N[ELUp[ν,k],wp],QWUp[ν,k,wp]}, 15, MyMin, True,1],If[k==Last[ksB1],endst<>"\\midrule\n",endeverystr]],{k,ksB1 }]
],{ν,νsB1}]];

Table 2 (remove semicolon to get LATEX code)

StringJoin[Table[StringJoin["\\multirow{"<>ToString[Length[ksB1]]<>"}{*}{$",
ToString[TeXForm[ν]],
"$}",
Table[StringJoin[MakeLineJ[k, {θUpInv[ν,Pi (k -1/2), wp],N[ELDown[ν,k],wp],QWDown[ν,k,wp]}, 15, MyMax, True,1 ],If[k==Last[ksB1],endst<>"\\midrule\n",endeverystr]],{k,ksB1 }]
],{ν,νsB1}]];

Figure accuracy of Bessel zeros bounds

νsB1Logs= (Range[28]-4)/4;
ksB1Logs=(Range[25]-1)/4;
wp=16;
xyticks={0,2,4,6};
errs=Flatten[Table[{νlog,klog,Log[10,θDownInv[10^νlog,Pi (Floor[10^klog] -1/2), wp]/θUpInv[10^νlog,Pi (Floor[10^klog] -1/2), wp]-1]}, {νlog,νsB1Logs}, {klog,ksB1Logs}],1];
(* figerrj=cleanContourPlot[ListContourPlot[errs,Contours->Table[c,{c,-3,-14,-1}],PlotLegends->BarLegend[Automatic,All],
Frame ->{{True,False},{True,False}},FrameLabel->MaTeX[{"\\log_{10}\\nu","\\log_{10} k"}],FrameTicks->{{xyticks, MaTeX[xyticks]}//Transpose,{xyticks, MaTeX[xyticks]}//Transpose}]] *)
figerrj=Labeled[cleanContourPlot[ListContourPlot[errs,Contours->Table[c,{c,-3,-14,-1}],PlotLegends->BarLegend[Automatic,All],
Frame ->{{True,False},{True,False}},FrameTicks->{{xyticks, MaTeX[xyticks]}//Transpose,{xyticks, MaTeX[xyticks]}//Transpose}]],
MaTeX[{"\\log_{10}\\nu","\\log_{10} k"}],{Bottom,Left},RotateLabel->True,Spacings->{0,-0.7}]

Export[SaveDir<>"figerrj.pdf", figerrj]

§B.2. Bounds for zeros of derivatives of Bessel functions

Pretty LATEX formatting

MakeLineJprime[k_, outbd_,bd_, l_, Myf_,showdiff_:False, dl_:1]:=Module[{str,minmax, b, stradd},
str= " & " <>ToString[k]<> " & "<> MyDecForm[outbd,l];
minmax=Myf[bd];
Do[
stradd= If[TrueQ[b==minmax],
 " & {\\color{red}"<> MyDecForm[b,l]<>"}" ,
 " & "<> MyDecForm[b,l]];
str=str<> stradd,{b,bd}];
If[showdiff, 
If[bd[[1]]!=minmax,str=str<> " & "<> ToString[TeXForm[ScientificForm[bd[[1]]/minmax-1,dl,NumberPoint->If[dl==1,"","." ]]]],str=str <>" & "]
 ];
str
];

Table 3 (remove semicolon to get LATEX code)

νstar0=ν/.FindRoot[zstar[ν]/Pi==1/2,{ν,1}];

StringJoin[Table[StringJoin["\\multirow{"<>ToString[Length[ksB1]]<>"}{*}{$",
ToString[TeXForm[ν]],
"$}",
Table[StringJoin[MakeLineJprime[k,If[ν>νstar0|| k>1,ϕUpInv[ν,Pi (k -1/2), wp],None], {ϕDownInv[ν,Pi (k -1/2), wp],ELprimeUp[ν,k,wp]}, 15, MyMin, True, 1],If[k==Last[ksB1],endst<>"\\midrule\n",endeverystr]],{k,ksB1 }]
],{ν,νsB1}]];
(* remove semicolon to get LaTEX for Table 3*)

Figure accuracy of Bessel derivative zeros bounds

errsprime=Flatten[Table[If[10^νlog>νstar0||Floor[10^klog] >1,
{νlog,klog,Log[10,θDownInv[10^νlog,Pi (Floor[10^klog] -1/2), wp]/θUpInv[10^νlog,Pi (Floor[10^klog] -1/2), wp]-1]}], {νlog,νsB1Logs}, {klog,ksB1Logs}],1];
errsprime=Select[errsprime,!TrueQ[#==Null]&];
(* figerrjprime=cleanContourPlot[ListContourPlot[errsprime,Contours->Table[c,{c,-3,-14,-1}],PlotLegends->BarLegend[Automatic,All],
Frame ->{{True,False},{True,False}},FrameLabel->MaTeX[{"\\log_{10}\\nu","\\log_{10} k"}],FrameTicks->{{xyticks, MaTeX[xyticks]}//Transpose,{xyticks, MaTeX[xyticks]}//Transpose}]] *)
figerrjprime=Labeled[cleanContourPlot[ListContourPlot[errsprime,Contours->Table[c,{c,-3,-14,-1}],PlotLegends->BarLegend[Automatic,All],
Frame ->{{True,False},{True,False}},FrameTicks->{{xyticks, MaTeX[xyticks]}//Transpose,{xyticks, MaTeX[xyticks]}//Transpose}]],MaTeX[{"\\log_{10}\\nu","\\log_{10} k"}],{Bottom,Left},RotateLabel->True,Spacings->{0,-0.7}]

Export[SaveDir<>"figerrjprime.pdf", figerrjprime]