Stokes diagram for Lanier Watkin's Regge pole problem

by W. Garrett Mitchener

Needs["NumericalCalculus`"]

This function is from LW’s original problem.  I have no idea what it represents.

Q[x_]=With[
{n=5,
aa=.02,
a=.02,
bb=.04,
b=.04,
Z=1,
k=2,
kk=0.8359,
l=-6.01806039558861+1.84413640038481I},
With[
{A=a*Z^(1/3),
B=b*Z^(2/3),
CC=-2*Z},Together[k^2-1*(CC/((x*(A*x+1)*(B*x+1)^2)))-((l*(l+1))/( (x^2)))]]]

zeros=z/.Solve[Q[z]==0,z]

poles=z/.Solve[Denominator[Q[z]]==0,z]

AlmostZero[z_]:=Abs[z]<(10^(-8))

AlmostEqual[x_,y_]:=AlmostZero[x-y]

ToMultiplicities takes a list of numbers {(Subscript[x, 1]),(Subscript[x, 2]),...} and produces a list {{(Subscript[r, 1]),(Subscript[n, 1])},...} representing the fact that (Subscript[r, k]) occurs (Subscript[n, k]) times in the original list of (Subscript[x, j])’s.  Since we’re dealing with approximate numbers, I can’t use ==.  Instead, I’m using AlmostEqual defined above in terms of a distance threshold.

ToMultiplicities[xs_]:=Map[
Function[nearlyEqualItems,
{Mean[nearlyEqualItems],Length[nearlyEqualItems]}],
GatherBy[xs,AlmostEqual]]

zerosWMult=ToMultiplicities[zeros]

polesWMult=ToMultiplicities[poles]/.{x_,k_}->{x,-k}

Unit[z_]=(z/(Abs[z]))

WKB approximation is:

ψ''+Q[z]ψ==0
ψ=(1/((Q[z])^(1/4)))(e^(±(∫_a^z)(Sqrt[Q[t]])dt))

Stokes lines are oscillatory:

(Sqrt[Q[z]])dz=imaginary

Anti-Stokes lines are growing or decaying:

(Sqrt[Q[z]])dz=real

Take a turning point a that's of order k.  We need to start the Stokes and anti-Stokes lines near it.

Q[a+w]=A (w^k)+...
Think: dz=ϵ (e^(I ω t))
(Sqrt[Q[a+dz]])dz=(Sqrt[A])(dz^(k/2))dz

Stokes line:

(dz^((k+2)/2))~same phase as~(I/(Sqrt[A]))=target

StokesStartingOffsets[{a_,k_}]:=Module[{delataZ},
If[k==-2,{},
With[
{targetSquared=Unit[((-1)/(NLimit[((Q[z])/(((z-a))^k)),z->a]))]},
deltaZ/.NSolve[(deltaZ^((k+2)))==targetSquared,deltaZ]]]]

StokesStartingOffsets[First[zerosWMult]]

To trace a Stokes line: Take a ϵ step in the dz direction, then integrate this ODE:

(dz/dt)=Unit[(I/(Sqrt[Q[a+z[t]]]))]

But to get the right sign on that (Sqrt[□]), we need for I/(Sqrt[Q]) to be in approximately the same direction as z.  In other words, (I/(Sqrt[Q]))/z should be approximately phase 0.

tStart=0; tEnd=10;

This function takes a turning point a that is a root or pole of degree k, a direction dirZ, and an offset size ϵ.  It computes the Stokes lines starting at an offset ϵ from the turning point in the dirZ direction as an InterpolatingFunction.

StokesLine[{a_,k_},dirZ_,ϵ_]:=
Module[{z},
With[
{sign=If[Re[I/(ϵ dirZ(Sqrt[Q[a+ϵ dirZ]]))]<0,-1,1]},
z/.NDSolve[{z'[t]==Unit[(I/(sign(Sqrt[Q[a+z[t]]])))],z[0]==ϵ dirZ},z,{t,tStart,tEnd}]//First]]

ComplexToPoint[z_]={Re[z],Im[z]}

This gives all Stokes lines for a turning point, tp={a,k} as for StokesLine:

AllStokesLines[tp_,ϵ_]:=Map[StokesLine[tp,#,ϵ]&,StokesStartingOffsets[tp]]

allSLEx=AllStokesLines[First[zerosWMult],0.01]

Interpolating functions are defined on a closed interval [(Subscript[t, min]),(Subscript[t, max])] (or a cross product of closed intervals).  The following utility function just reads the endpoints of the domain out of an InterpolatingFunction object:

IFDomain[InterpolatingFunction[{{tMin_,tMax_}},___]]:={tMin,tMax}

(I don’t absolutely need that because I use the same (Subscript[t, min]) and (Subscript[t, max]) all the way through this notebook, but I didn’t know if I’d be able to do that when I was writing the worksheet.)

This makes a list of pictures with the given turning point and a Stokes line emanating from it.  (I could also do this with a single call to ParametricPlot if I was sure that I could use the same time interval for each one, but again, I wasn’t sure if that would be possible, so I set this up to draw one Stokes line at a time.)

MakeStokesLines[tp_,ϵ_]:=
Map[
Function[sl,
With[
{domain=IFDomain[sl]},
ParametricPlot[ComplexToPoint[First[tp]+sl[t]],{t,domain[[1]],domain[[2]]},
PlotStyle->{Hue[0.4]}
]
]
],
AllStokesLines[tp,ϵ]]

PlotStokesLines[tp_,ϵ_,options___]:=Show[MakeStokesLines[tp,ϵ],options]

PlotStokesLines[First[zerosWMult],0.01,
PlotRange->{-5,5},
AspectRatio->3/4]

Now for anti-Stokes lines, things are almost the same except we need:

(dz^((k+2)/2))~same phase as~(1/(Sqrt[A]))=target

AntiStokesStartingOffsets[{a_,k_}]:=Module[{deltaZ},
If[k==-2,{},With[
{targetSquared=Unit[(1/(NLimit[((Q[z])/(((z-a))^k)),z->a]))]},
deltaZ/.Solve[(deltaZ^((k+2)))==targetSquared,deltaZ]]]
]

AntiStokesLine[{a_,k_},dirZ_,ϵ_]:=
Module[{z},
With[
{sign=If[Re[1/(ϵ dirZ(Sqrt[Q[a+ϵ dirZ]]))]<0,-1,1]},z/.NDSolve[{z'[t]==Unit[(1/(sign(Sqrt[Q[a+z[t]]])))],z[0]==ϵ dirZ},z,{t,tStart,tEnd}]//First]
]

AllAntiStokesLines[tp_,ϵ_]:=Map[AntiStokesLine[tp,#,ϵ]&,AntiStokesStartingOffsets[tp]]

MakeAntiStokesLines[tp_,ϵ_]:=
Map[
Function[sl,
With[
{domain=IFDomain[sl]},
ParametricPlot[ComplexToPoint[First[tp]+sl[t]],{t,domain[[1]],domain[[2]]},
PlotStyle->{Hue[0]}
]
]
],
AllAntiStokesLines[tp,ϵ]]

PlotAntiStokesLines[tp_,ϵ_,options___]:=Show[MakeAntiStokesLines[tp,ϵ],options]

PlotAntiStokesLines[First[zerosWMult],0.01,
AspectRatio->3/4]

MakeStokesDiagram[tps_,ϵ_]:=Join@@Map[{MakeStokesLines[#,ϵ],MakeAntiStokesLines[#,ϵ]}&,
tps]

PlotStokesDiagram[tps_,ϵ_,options___]:=Show[MakeStokesDiagram[tps,ϵ],
options]

spPic=PlotStokesDiagram[Join[zerosWMult,polesWMult],0.0001,
AspectRatio->3/4,
PlotRange->{-5,5}]

tpPic=ListPlot[
{ComplexToPoint/@zeros,ComplexToPoint/@poles},
PlotMarkers->{Automatic,Medium}]

Show[{spPic,tpPic},
Frame->True,Axes->False]