Aharonov-Bohm effect

Quantum Mechanics
K. Konishi, G. Paffuti

Solve the Schrödinger equation for the Aharonov-Bohm process,  with a line-shaped solenoid or a finite cylinder-shaped solenoid.

ClearAll["Global`*"];
Off[General::spell, General::spell1];

Introduction

In quantum mechanics the Schrödinger equation for a charged particle in an electromagnetic field is

i ℏ(∂/(∂t))ψ[t,x] = (1/(2m))((((ℏ/i) ∇ -  (e/c)A))^2) ψ[t,x] + e Φ ψ[t,x];

This equation is invariant under the gauge transformation

ψ-> Exp[i (e/(ℏc)) Λ]ψ;     A-> A + ∇Λ  ; Φ-> Φ - (1/c)((∂Λ)/(∂t))  ;

generalizing the invariance of the classical Maxwell equations. While in classical mechanics Newton's equations can be written in terms of the gauge invariant electric and magnetic fields:

E = - ∇Φ -(1/c) ((∂A)/(∂t)) ;   B = ∇⋀ A  ,

the use of gauge potentials A, Φ is unavoidable in quantum mechanics, see eq.(1).  This raises the question whether there exist observables effects not directly related to classical e.m. fields. For the sake of simplicity let us consider a static case, with Φ = 0, and A  independent of t. Eq.(3) implies E = 0. If in a given domain D we know A, B is given by its rotor. The nontrivial point is that while the value of B is fixed by the local properties of A, i.e. by its derivatives, the very existence of a globally well defined vector field A can give rise to new mathematical structures.  A vector can be integrated along lines and the Stokes theorem applied to a surface Σ with contour ∂Σ reads

(∮_(∂Σ)) A dx = (∫_Σ) B dσ .

Eq.(4) specialized to a small circle around an arbitrary point x implies eq.(3). The interesting case is the one in which the region Σ is not simply connected (i.e., it has a hole, roughly speaking). As an example consider the 2-dimensional case shown below

Show[Graphics[{{RGBColor[0.718227, 0.794781, 0.937957],Disk[{0,0},2]},{RGBColor[1,1,1],Disk[{0,0},0.9]},Circle[{0,0},2],Circle[{0,0},0.9],
{RGBColor[1, 0, 0],Disk[{0,0},0.5]},
Text[Style["Σ",FontFamily->"Times",FontSize->14],{1.5,0}],
Text[Style["(Subscript[C, 1])",FontFamily->"Times",FontSize->14],{2.02,0.8}],
Text[Style["\!\(\*SubscriptBox[\(C\), \(\(\\ \)\(2\)\)]\)",FontFamily->"Times",FontSize->14],{0.86,0.7}],
Arrow[{{2,0},{2,0.2}}],Arrow[{{0,0.9},{-0.2,0.9}}]},
AspectRatio->Automatic]]

(∮_(Subscript[C, 1])) ((Subscript[A, x])dx + (Subscript[A, y])dy)- (∮_(Subscript[C, 2])) ((Subscript[A, x])dx + (Subscript[A, y])dy)= (∫_Σ)(Subscript[B, z])dx dy

Suppose that B = 0. The previous relation implies invariance under the deformation of the contour and, if (Subscript[C, 2]) can be shrunk to zero the line integral vanishes. But if the region in which A is defined is not simply connected, e.g. if the red circle in the figure is excluded, we have a gauge invariant quantity, the integral around a closed circuit, which can be nontrivial and unrelated to a local magnetic field. The question poses itself: are quantum mechanical predictions affected by the value of this line integral?. The Ahronov-Bohm effect answers affirmatively this question: there are observable effects. The question has been settled beyond any doubt in relatively recent  experiments  (see the main text for the discussion),  ref.[1]

Description of the system

Let us consider an infinitely long and impenetrable cylindrical solenoid of radius R. We will use cylindrical coordinate (r,θ, z)with solenoid axis as z axis.  Φ(r) is the magnetic flux trough a circle of radius r . From eq.(4) the vector potential can be taken directed along θ and Φ(r) = 2π r (Subscript[A, θ]). In this gauge the vector potential vanishes for r -> ∞. Outside the solenoid the magnetic field vanishes, as is easily checked in cartesian coordinates

((Subscript[A, x]),(Subscript[A, y])) = ((Φ[R])/(2π r))(-Sin[θ],Cos[θ]) = ((Φ[R])/(2π (r^2)))(-y,x)

BfieldOutside = Simplify[(∂_x)((x/((x^2)+(y^2))))-(∂_y)(((-y)/((x^2)+(y^2))))]

A limiting case is R→0 with Φ(R) fixed, i.e. an infinitely thin solenoid. This is the original example studied by Aharonov and Bohm. We will study the scattering of a monochromatic wave impinging the solenoid from positive x direction. In the first part of this notebook we consider the limiting case R→0.

Boundary conditions and gauge invariance

In cylindrical coordinates the Schrödinger equation outside the solenoid reads

((∂^2)/(∂ (r^2)))ψ+ (1/r)(∂/(∂r))ψ + (1/(r^2))((((∂/(∂θ))- iα))^2)ψ + (k^2)ψ = 0.

k is the wave number, the energy of the incoming wave is E = (((ℏ k))^2)/2m. Using the relation 2πr (Subscript[A, θ])= Φ(R) we have introduced the parameter

α = (e/(ℏc))((Φ[R])/(2π))

The impenetrability of the solenoid amounts to the condition ψ(R) = 0. To compute the scattering solution we have to fix the asymptotic behavior of ψ. In general ψ can be decomposed as an incident and a scattered wave:

ψ = (Subscript[ψ, in])+ (Subscript[ψ, sc])

This decomposition is not necessarily true in the forward direction, θ = π in our case, where interference effects between transmitted and incident waves can be important. In usual cases (Subscript[ψ, in])would be written as (Subscript[ψ, in])= Exp[- i k x] corresponding to a plane wave incoming from the positive x direction. In the present case the incoming wave is distorted by the long range nature of the vector field. Physically we demand that for x very large and positive ψ  describes  a particle with wave number k. The gauge invariant current describing the flux of particles is

j = (1/(2m))((SuperStar[ψ])(p - (e/c)A)ψ - ψ(p - (e/c)A)(SuperStar[ψ]))=(ℏ/(2m i))((SuperStar[ψ])∇ψ - ψ∇(SuperStar[ψ]))- (e/mc)A|ψ((|)^2).

Locally

A = ((Φ[R])/(2π (r^2)))(-y,x)=((Φ[R])/(2π ))∇θ; (e/mc)A =(ℏ/m)α ∇θ .

The distorted wave must cancel the additional factor with respect the plane wave, i.e.

(Subscript[ψ, in])= Exp[-i k x + α θ].

The scattered wave must be a diverging cylindrical wave, i.e. asymptotically

ψ =Exp[-i k x + α θ] + ((f[θ])/( (r^(1/2))))Exp[i k r].

f is the scattering amplitude and the cross section for the scattering is

dσ = |f(θ)((|)^2)dθ.

Monodromy and gauge invariance

The solution of the Schrödinger equation must be a function, i.e. must be a single valued function (monodromy) crossing the negative x axis (θ= π). This means that the discontinuity in (Subscript[ψ, in]) along this line must be cancelled by (Subscript[ψ, sc]). Formally the α dependence on eq.(5) can be reabsorbed by a change of variables (gauge transformation)

ψ=(Subscript[ψ, α=0])Exp[i α θ].

This is not an admissible gauge transformation, it is not single valued and in effect the solution given below has not this form. Let us note that eq.(10) is formally of the form

ψ=(Subscript[ψ, α=0])Exp[i α∫ (∂_i)θ (dx^i)]=(Subscript[ψ, α=0])Exp[i (e/(ℏc))∫ (Subscript[A, i]) (dx^i)]

this formula is inconsistent because the vector field is only locally defined as (∂_i)θ, i.e. is not globally a gradient (a gauge transformation).

Analytical solution

The system is symmetric under rotation around z axis, the solution at fixed energy can be expanded in eigenfunctions of (Subscript[L, z])

ψ = (∑_(m=-∞)^∞)(e^(i m θ))(Subscript[f, m])[r].

At an even simpler level ψ must be periodic in θ (monodromy) and eq.(11) is just its Fourier expansion. Substituting in eq,(5) the radial functions (Subscript[f, m]) must satisfy

(((d^2)(Subscript[f, m]))/(d (r^2)))+ (1/r)((Subscript[df, m])/dr)ψ - (1/(r^2))(((m- α))^2)(Subscript[f, m]) + (k^2)(Subscript[f, m]) = 0.

The change of variables z = k r reduces this equation to the defining equation for Bessel functions and we have for the generic solution

ψ = (∑_(m=-∞)^∞)(e^(i m θ))( (Subscript[A, m])(Subscript[J, m-α])[kr]+ (Subscript[B, m])(Subscript[J, -(m-α)])[kr]).

The condition ψ(0) = 0 (we are discussing the limiting case R=0 for the cylinder) implies that only bessel functions of positive order can appear in the expansion. i.e.

ψ = (∑_(m=-∞)^∞)(e^(i m θ)) (Subscript[a, m])(Subscript[J, |m-α|])[kr]

The asymptotic condition  fixes the coefficients (Subscript[a, m]). The simplest way to proceed is to consider each angular momentum separately. If we study ψ at fixed m we must find an incoming wave corresponding to a plane wave and an outgoing cylindrical wave due to scattering. For a plane wave the following expansion holds:

Exp[-i k r Cos[θ]]=  (∑_(m=-∞)^∞)(((-i))^m) (Subscript[J, m])[kr] Exp[i m θ].

For Bessel functions asymptotically

(Subscript[J, ν])[kr]∼ (Sqrt[(2/(π kr))])Cos[ kr - (π/2)ν - (π/4)] = (Sqrt[(2/(π kr))]) (1/2)(Exp[i(kr - (π/2)ν - (π/4))] + Exp[-i(kr - (π/2)ν - (π/4))]).

The convergent part of the plane wave and ψ is respectively, using -i = Exp[-i π/2]

(Sqrt[(2/(π kr))]) (1/2) (((-i))^m) Exp[-i(kr - (π/2)m - (π/4))]=(Sqrt[(2/(π kr))]) (1/2)Exp[-i(kr - (π/4))];
(Sqrt[(2/(π kr))]) (1/2)(Subscript[a, m ])Exp[-i(kr - (π/2)|m-α| - (π/4))];

it follows  (Subscript[a, m]) = (((-i))^(|m-α|))and finally

ψ = (∑_(m=-∞)^∞)(e^(i m θ)) (((-i))^(|m-α|))(Subscript[J, |m-α|])[kr].

Let us note that α dependence is contained in the argument of the Bessel functions, not in θ-dependent part, as in the improper gauge transformation.

Asymptotic expansion

The nontrivial task is now to show that eq.(15) asymptotically behaves as requested in eq.(8). Let us sketch the proof, leaving the details to ref.[2]. Substitution of asymptotic expansion (14) in (8) is not permitted. In effect (14) is valid only when the argument of the Bessel function is greater than its order, kr > ν, i.e. it is not a uniform expansion. We have to sum over ν and the contribution of higher terms in the sum is not negligible, i.e. the limit r→∞ does not commute with the infinite summation.
The way out is to introduce an integral representation for the Bessel functions:

(Subscript[J, ν])[z] = (1/(2π)) (∫_C) Exp[i (ν t - z Sin[t])] dt .

The complex path C must close to infinity. Writing t = u + i v, Sin[t] = Sin[u]Cos[v] + i Cos[u]Sinh[v], the convergence is obtained, for z > 0, where Cos[u]Sinh[v]<0, i.e. in the regions indicated in the following figure

pl1=Plot[If[Abs[x]>0.8π,100, 1/(x-0.8π)^2*1/(x+0.8π)^2],{x,-1.2 π,1.2 π},PlotRange-> {{-1.2 π,1.2 π},{-5,5}},
PlotStyle-> {GrayLevel[0],Thickness[.01]},Ticks-> {{-π,((-π)/2),(π/2),π},{}},AxesLabel-> {"Re[t]","Im[t]"},
(* AspectRatio-> Automatic, *)
Epilog-> {Text["C",{0.9π,3}],Dashing[{0.015,0.015}],Line[{{-(π/2),-5},{-(π/2),5}}],Line[{{(π/2),-5},{(π/2),5}}]}
]

Substitution of (16) into (15) gives two geometric series ( one for each sign of m-α) which can be easily summed producing

ψ = (1/(2π)) ∫dt Exp[-i kr Sin[t]] (((Exp[i( [α]θ +(t-π/2)(α-[α]))])/(1-Exp[i(t-π/2 -θ)]))+((Exp[i( ([α]+1)θ +(t-π/2)([α]-α+1))])/(1-Exp[i(t-π/2 +θ)])))

[α] is the integer part of α. For large values of r the integral can be approximated by the saddle point method. Saddles points are given by the stationary phase condition d/dt Sin[t] = Cos[t] = 0, i.e. t = +π/2 and t=-π/2. These saddle points must be crossed along vertical lines, so the contour must be deformed in the lower complex plane. The integrand has poles for t = π/2 + θ (first term) and t = π/2 - θ (second term). The poles contribution must be subtracted and this contribution gives rise to the incident wave

(Subscript[ψ, pole])= Exp[-i k rCos[θ] + i α θ];

The saddle point contribution gives

(Subscript[ψ, SP])= ((Sin[π α])/((((2π i kr))^(1/2))Cos[(θ/2)]))(((-1))^([α]))Exp[i([α]+(1/2))θ] Exp[i k r];

i.e. the scattered wave. Putting all terms together the asymptotic expansion of ψ is

ψ -> Exp[-i k rCos[θ] + i α θ] + ((Sin[π α])/((((2π i kr))^(1/2))Cos[(θ/2)]))(((-1))^([α]))Exp[i([α]+(1/2))θ] Exp[i k r];  0 < θ < π.

As it is evident the expansion is not valid for θ = π (the poles collapse to the saddle points). We quote the result for this particular case

ψ[r,θ=π] -> Exp[i(kr+ [α]π)]Cos[π(α - [α])].

This particular case shows that ψ vanishes at θ=π for half-integer values of α. This property follows also from the relations which follows from the expression (15):

ψ[r,θ;-α] = ψ[r,-θ;α]; ψ[r,θ;α+N]=Exp[iN θ]ψ[r,θ;α].

Finally let us mention that for α = 1/2 the series can be resummed in closed form:

ψ[r,θ;α = (1/2)] =(2/(Sqrt[π]))Exp[i( θ - kr Cos[θ] - (π/4))] (∫_0^((((2kr))^(1/2))Cos[θ/2]))Exp[i (t^2)]dt.

Cross section

From eq.(8,9) and (17) the differential cross section is easily obtained:

dσ = (((Sin^2)[π α])/(2π k (Cos^2)[θ/2]))dθ.

Finite cylinder

In the case of finite cylinder radius the solution is easily obtained:

ψ[r,θ] = ψ = (∑_(m=-∞)^∞)(e^(i m θ)) (((-i))^(|m-α|))((Subscript[J, |m-α|])[kr]+ (Subscript[b, m])((Subscript[H, |m-α|])^((1)))[kr]);           (Subscript[b, m])= - (((Subscript[J, |m-α|])[kR])/(((Subscript[H, |m-α|])^((1)))[kR])).

In effect (21) has the same convergent cylindrical wave as (17), satisfying in this way asymptotic condition. (H^((1)))are Hanckel functions of the first kind, representing asymptotically a divergent wave:

((Subscript[H, ν])^((1)))[z] = (Subscript[J, ν])[z] + i (Subscript[Y, ν])[z];              ((Subscript[H, ν])^((1)))[z]∼ (Sqrt[(2/(π z))])Exp[i( z - (π/2)ν - (π/4))];

The coefficients (Subscript[b, m]) are fixed by the boundary condition ψ = 0 on the boundary.

Numerical Solution

For our numerical study we choose units such that k=1. In this way R >> 1 means small wavelength, R << 1 long wavelengths. The limiting case R→0 belongs to the long wavelength regime.

For numerical purposes we define a maximum allowed value for the angular momentum, mMax, and compute together all functions needed, using Mathematica function Table:

phase[θ_,α_,mMax_]:= Table[Exp[I m θ -I (π/2)Abs[m-α]] ,{m,-mMax,mMax}];
vecJ[r_,α_,mMax_]:= Table[BesselJ[Abs[m-α],r],{m,-mMax,mMax}];
vecY[r_,α_,mMax_]:= Table[BesselY[Abs[m-α],r],{m,-mMax,mMax}];
vecH1[r_,α_,mMax_]:= Table[BesselJ[Abs[m-α],r]+I  BesselY[Abs[m-α],r],{m,-mMax,mMax}];

The coefficients (Subscript[b, m]) of eq.(22) are

bm[α_,R_,mMax_]:= -vecJ[R,α,mMax]/vecH1[R,α,mMax];

The wave function for the general case and for the simplest situation R=0 are respectively

ψ[r_,θ_,α_,R_,mMax_]:= If[r< R,0.,
Plus@@(phase[θ,α,mMax](vecJ[r,α,mMax]+bm[α,R,mMax]*vecH1[r,α,mMax]))]

ψAB[r_,θ_,α_,mMax_]:= If[r==0 ,0.,
Plus@@(phase[θ,α,mMax]*vecJ[r,α,mMax])]

As an example we plot the intensity of the wave at fixed distance d along x axis (simulating a screen at distance d). The distance from the origin is r = -d/Cos[θ] in the region θ> π/2. We give on the same plot a graph of |ψ((|)^2) and a ContourPlot to give a flavour of the intensity.

plots[ψ_,d_,mMax_]:= Table[Show[Plot[1+((Abs[Evaluate[ψ[-(d/(Cos[θ])),θ,0.1k,mMax]]])^2),{θ,(π/2)+1,((3π)/2)-1},PlotRange-> All,Ticks-> {{4/5π,π,6/5 π},{}},PlotPoints-> 25,
MaxRecursion-> 0,
PlotLabel-> "α="<>ToString[NumberForm[0.1k,{3,2}]]<>
"  M="<>ToString[NumberForm[mMax,{3,0}]]],
 ContourPlot[((Abs[Evaluate[ψ[-(d/(Cos[θ])),θ,0.1k,mMax]]])^2),{θ,(π/2)+1,((3π)/2)-1},{y,0,1},ContourLines->False, 
(* PerformanceGoal->"Speed", *)
PlotPoints-> {20,3},MaxRecursion-> 0,
ColorFunction->GrayLevel,
Frame-> False]],
{k,2,10}]

GraphicsGrid[Partition[plots[ψAB,30,100],3],ImageSize-> Medium]//Timing

The reader can explore different possibilities, changing the value of mMax and of d. We leave as an exercise the comparison with the expression (20) for the particular case α = 1/2.

A faster procedure

Drawing several figures for high values of mMax can be a quite long procedure, the time is spent mostly in the computation of Bessel functions. An alternative is to compute Bessel functions using the following recursion relation

(Subscript[Z, ν+1])[x]-((2ν)/x)(Subscript[Z, ν]) + (Subscript[Z, ν-1])= 0.

It is known, see ref.[3], that for numerical stability the recurrence relation must be used in increasing order for bessel functions of the second kind (Subscript[Y, ν])and in decreasing order for Bessel functions (Subscript[J, ν]) (for positive ⋁).

For a given value of mMax our sequences contains (2 mMax+1) terms. We can take Floor[α]+1 as a separator between the 2 subsequences, distinguishing the case m>α and m<α.
1) Floor[α]+1 ≤ m ≤ mMax;  we have ind1 = mMax - Floor[α] terms. The index of the Bessel functions runs as 
   1+Floor[α] - α = 1-frac[α] ≤ m-α ≤ mMax - α
2) -mMax ≤ m ≤ Floor[α], we have ind2 = Floor[α] + mMax + 1 terms. The index of the Bessel functions runs as
     frac[α] ≤ α - m ≤  mMax + α

frac[α] ≡ α - Floor[α] and we assume α > 0.

It is now easy to write the needed sequences for the phases and for the Bessel functions:

phaseL[θ_,α_,mMax_]:=Flatten[ {Table[Exp[I m θ -I (π/2)(m-α)] ,{m,Floor[α]+1,mMax}],
Table[Exp[I m θ -I (π/2)(α-m)] ,{m,Floor[α],-mMax,-1}]}]

jL[r_,α_,mMax_]:= 
Module[{frac = α-Floor[α],ind1=mMax-Floor[α],ind2=mMax+Floor[α]+1},
JList1=Table[1.,{ind1}];JList2=Table[1.,{ind2}];
JList1[[ind1]]= BesselJ[mMax-α,r];
JList1[[ind1-1]]= BesselJ[mMax-α-1,r];
Do[ JList1[[i]] = -JList1[[i+2]]+(2 (i+1-frac))JList1[[i+1]]/r;,{i,ind1-2,1,-1}];
JList2[[ind2]]= BesselJ[mMax+α,r];
JList2[[ind2-1]]= BesselJ[mMax+α-1,r];
Do[ JList2[[k]] = -JList2[[k+2]]+(2 (k+1-1+frac))JList2[[k+1]]/r;,{k,ind2-2,1,-1}];Flatten[{JList1,JList2}]];

yL[r_,α_,mMax_]:= 
Module[{frac = α-Floor[α],ind1=mMax-Floor[α],ind2=mMax+Floor[α]+1},
YList1=Table[1.,{ind1}];YList2=Table[1.,{ind2}];
YList1[[1]] = BesselY[1-frac,r];
YList1[[2]] = BesselY[2-frac,r];
Do[ YList1[[i]] = -YList1[[i-2]]+(2 (i-1-frac))YList1[[i-1]]/r;,{i,3,ind1}];
YList2[[1]] = BesselY[frac,r];
YList2[[2]] = BesselY[frac+1,r];
Do[ YList2[[i]] = -YList2[[i-2]]+(2 (i-2+frac))YList2[[i-1]]/r;,{i,3,ind2}];
Flatten[{YList1,YList2}]];

Hanckel's functions and coefficients (Subscript[b, m]) are given by

hL[r_,α_,mMax_]:= jL[r,α,mMax]+ I yL[r,α,mMax];
bmL[α_,R_,mMax_]:= -jL[R,α,mMax]/hL[R,α,mMax];

The AB-function and the ψ function for the case R > 0 are written as above:

ψABL[r_,θ_,α_,mMax_]:= If[r== 0,0,
Plus@@(phaseL[θ,α,mMax]*jL[r,α,mMax])]

ψL[r_,θ_,α_,R_,mMax_]:= If[r< R,0,
Plus@@(phaseL[θ,α,mMax](jL[r,α,mMax]+bmL[α,R,mMax]*hL[r,α,mMax]))]

The reader is warned that this recursive computation is instable for too high values of r. We leave as an exercise the modification of ψABL and ψL to take into account this effect (e.g. one can use asymptotic expansion in this region).

The previous graph is obtained much faster.

GraphicsGrid[Partition[plots[ψABL,30,100],3],ImageSize-> Medium]//Timing

Contour plots

Using the previous defined functions it is quite easy to plot the intensity of the wave in all the  x-y plane

xmax = 20; ymax = 20;

makeContour[α_,R_,mMax_]:=
DensityPlot[Abs[Evaluate[ψL[(Sqrt[x^2+y^2]),ArcTan[x,y],α,R,mMax]]],{x,-xmax,xmax},{y,-ymax,ymax},
 PlotPoints-> 40 {1,1},MaxRecursion-> 0,  
ColorFunction->GrayLevel,
AspectRatio-> Automatic,
PlotLabel-> "α="<>ToString[NumberForm[α,{2,1}]]<>" R="<>ToString[NumberForm[R,{3,2}]]<>
" m="<>ToString[NumberForm[mMax,3]],
Epilog->{{GrayLevel[0],Disk[{0,0},R]}}
];

scatteringPlot[R_,mMax_]:=Table[makeContour[0.1*k,R,mMax],{k,1,9}];

We present as an example plots for R=2

GraphicsGrid[Partition[scatteringPlot[2,30],3],ImageSize-> Medium]//Timing

Plane of constant phase

In ref.[4] has been pointed out that while ψ must be single valued, its phase can have jumps. A singularity in the phase can occur only for |ψ|=0, because ψ itself is continuous. If we analyze constant phase planes, i.e. front-waves, we can discover such critical points, one of which is the origin, for the AB-case.

To give a graphical view of the phenomenon we can plot, as an example, the line of constant Arg[Exp[-i α π/2] ψ]. The function Arg is discontinuous, if we try to select a particular numerical value by Arg = value we usually found multiple solutions, one correct, the other is an artifact of the discontinuity. Take as an example the asymptotic form of ψ:

ψas[x_,y_,α_]:= Module[{r,θ},
r = Sqrt[(x^2)+(y^2)];θ = ArcTan[x,y];
αI = IntegerPart[α];
Exp[-I x + I α θ]+ (((E^(I r))Sin[π α])/((Sqrt[2π I r])Cos[θ/2]))(((-1))^(αI))Exp[I(αI + 1/2)θ]
]

For a fixed value of y, eg. 3, we have

Plot[Arg[(E^(-I β π/2))ψas[x,3,β]]/.β-> 1,{x,-10,10}]

The vertical spurious vertical part of the plot intercept a value between -3 and 3 giving rise to a spurious solution. The discontinuity of Arg is at π so we adopt the following trick: multiply the function with Exp[i π] and we look for solutions of Arg = 0.95 π. We find two very close solutions (corresponding to turning points of previous graph) which are almost indistinguishable on the graphs.

The computation with the full form of ψ is quite time consuming, and for graphics illustration is sufficient to work with its asymptotic form. A simple possibility is to substitute the asymptotic form of ψ everywhere except near θ = π, where it breaks down.

ψ1[x_,y_,α_,mMax_]:=If[x^2+y^2 ≥ 0.1 && π - Abs[ArcTan[x,y]]≥ 4/Sqrt[x^2+y^2], Arg[ (E^(-I α π/2+I π))ψas[x,y,α]],Arg[ (E^(-I α π/2+I π))ψABL[(Sqrt[(x^2)+(y^2)]),ArcTan[x,y],α,mMax]]];

β = 1; mMax = 30;

We use Mathematica  function ContourPlot to select the desired value of the phase.

ContourPlot[ψ1[x,y,β,mMax]== 0.95π ,{x,-xmax,xmax},{y,-xmax,xmax},PlotPoints-> 150{1,1},MaxRecursion-> 1,ContourStyle-> Thickness[.01],Epilog-> {Thickness[.01],Line[{{-xmax,0},{xmax,0}}]}]//Timing

We draw an animation to show the methamorphosis of the waves as α change. First we select the desired values for α:

alpha = Table[1+0.05 *k,{k,0,29}];

then use ListAnimate. The computation is quite long, the reader can use a shorter list for α if a faster result is required.

ListAnimate[
Table[ContourPlot[ψ1[x,y,alpha[[k]],mMax]== 0.95π ,{x,-xmax,xmax},{y,-xmax,xmax},PlotPoints-> 150{1,1},MaxRecursion-> 1,ContourStyle-> Thickness[.005],Epilog-> {Thickness[.01],Line[{{-xmax,0},{xmax,0}}]}],{k,1,Length[alpha]}],
AnimationRunning-> False]//Timing

References