Stokes diagram for the overdense barrier

by W. Garrett Mitchener

<< NumericalMath`

Q[z_] = -(1 - z) (1 + z)

Plot[Q[z], {z, -3, 3}]

Q[1 + ϵ] // Expand

Q[-1 + ϵ] // Expand

sl1NSol = 
      With[{ϵ = 1/100}, 
        NDSolve[{z'[t] == 
              I/Sqrt[Q[z[t]]], 
z[0] == 
              1 + Exp[2  π I/6] ϵ}, z, {t, 0, 10}]]

ComplexToPoint[w_] = {Re[w], Im[w]}

z[1] /. sl1NSol

ParametricPlot[
      Evaluate[ComplexToPoint[z[t]] /. sl1NSol], 
{t, 0, 
        10}]

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

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

poles = {}

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

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

ToMultiplicities[{}, soFar_] := soFar

ToMultiplicities[{x_, xs___}, soFar_] := 
      ToMultiplicities[{xs}, 
With[
{search = 
              Position[soFar, ({y_, c_} /; AlmostEqual[y, x])]}, 
          If[
Length[search] == 
              0, 
Append[
              soFar, {x, 1}], 
ReplacePart[soFar, 
              Part[soFar, 
                  search[[1, 
                    1]]] /. {y_, c_} -> {y, 
                    c + 1}, 
search[[1, 
                1]]]
]
]
]

Position[{{1, 2}, {2, 3}, {3, 1}}, ({y_, c_} /; 
          AlmostEqual[y, 2])]

Part[{{1, 2}, {2, 3}, {3, 1}}, 2]

ToMultiplicities[{1, 1, 2, 3, 4}, {}]

ToMultiplicities[xs_] := ToMultiplicities[xs, {}]

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\(\@ 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 turining 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_}] := 
      If[k == -2, {}, 
        With[
{targetSquared = 
              Unit[((-1)/NLimit)[Q[z]/(((z - a))^k), z -> a]]}, 
deltaZ /. 
            NSolve[(deltaZ^((k + 2))) == targetSquared, 
              deltaZ]]]

StokesStartingOffsets[{-1, 1}]

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.

So, take this example.  This one needs the -(Sqrt[□]).

I/(deltaZ (Sqrt[ Q])[-1 + deltaZ]) /. 
      deltaZ -> 
        0.01 StokesStartingOffsets[{-1, 
                1}][[1]]

tStart = 0;  tEnd = 10;

StokesLine[{a_, k_}, 
        dirZ_, ϵ_] := 
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]

slEx = 
      StokesLine[{-1, 
          1}, StokesStartingOffsets[{-1, 
              1}][[3]], 
        0.01]

slEx[0.2]

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

ParametricPlot[ComplexToPoint[-1 + slEx[t]], {t, 0, 10}, 
      PlotRange -> {{-2, 2}, {-2, 
            2}}, 
AspectRatio -> 
        1, 
AxesOrigin -> {-2, -2}]

This gives all Stokes lines for a turning point:

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

allSLEx = AllStokesLines[First[zerosWMult], 0.01]

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

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

PlotStokesLines[tp_, ϵ_] := 
      Show[MakeStokesLines[tp, ϵ], 
        DisplayFunction -> $DisplayFunction]

Show[
      PlotStokesLines[First[zerosWMult], 
        0.01], 
PlotRange -> {{-2, 2}, {-2, 
            2}}]

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_}] := 
      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_, ϵ_] := 
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]

alEx = 
      AntiStokesLine[{-1, 
          1}, AntiStokesStartingOffsets[{-1, 
              1}][[1]], 
        0.01]

ParametricPlot[ComplexToPoint[-1 + alEx[t]], {t, 0, 10}, 
      PlotRange -> {{-2, 2}, {-2, 
            2}}, 
AspectRatio -> 
        1, 
AxesOrigin -> {-2, -2}]

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

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

PlotAntiStokesLines[tp_, ϵ_] := 
      Show[MakeAntiStokesLines[tp, ϵ], 
        DisplayFunction -> $DisplayFunction]

Show[
      PlotAntiStokesLines[First[zerosWMult], 
        0.01], 
PlotRange -> {{-2, 2}, {-2, 
            2}}, 
AspectRatio -> 1]

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

PlotStokesDiagram[tps_, ϵ_] := 
      Show[MakeStokesDiagram[
          tps, ϵ], 
DisplayFunction -> $DisplayFunction]

spPic = 
      Show[PlotStokesDiagram[Join[zerosWMult, polesWMult], 
          0.01], 
PlotRange -> {{-2, 2}, {-2, 
              2}}, 
AspectRatio -> 1]

tpPic = 
      ListPlot[ComplexToPoint /@ (Join[zeros, poles]), 
        PlotStyle -> PointSize[0.025]]

Show[{spPic, 
        tpPic}, 
AxesOrigin -> {-2, -2}]