This notebook contains additional material to the Mathematica GuideBooks  (Springer-Verlag, New York, 2004) http://www.MathematicaGuideBooks.org . It uses the stylesheet of the Mathematica GuideBooks, which is embedded.

Symbolics Volume, Chapter 1

The Green’s function (Subscript[G, k])(y,x;(V^((n)))) for the time-independent Schrödinger equation with a sum of delta function potentials (V^((n)))(z)=(∑_(k=1)^n)(Subscript[V, k])δ(z-(Subscript[z, k])) is defined through the equation -(∂^2)(Subscript[G, k])(y,x;(V^((n))))/∂(x^2)+V(x)(Subscript[G, k])(y,x;(V^((n))))=(k^2)(Subscript[G, k])(y,x;(V^((n)))) [✶] and fulfills the recursion [✶], [✶], [✶]

(Subscript[G, k]^((0)))(y,x;(V^((0))))=-(i/k)(e^(i k |y-x|))
(Subscript[G, k]^((j)))(y,x;(V^((j))))=(Subscript[G, k]^((j-1)))(y,x;(V^((0))))+(((Subscript[V, j])(Subscript[G, k]^((j-1)))(y,(Subscript[z, j]);(V^((j-1))))(Subscript[G, k]^((j-1)))((Subscript[z, j]),x;(V^((j-1)))))/(1-(Subscript[V, j])(Subscript[G, k]^((j-1)))((Subscript[z, j]),(Subscript[z, j]);(V^((j-1))))))
(Subscript[G, k])(y,x;(V^((0))))=(Subscript[G, k]^((n)))(y,x;(V^((0)))).

The poles of (Subscript[G, k])(y,x;(V^((n)))) in the complex k-plane indicate bound and resonant states. Make an animation that shows how the poles move in the complex k-plane for (V^((8)))(z) with the parameters (Subscript[V, j])=(e^(i ω)), 0≤ω≤2π, (Subscript[z, j])=(j-1)λ, λ∈ℝ as a function of ω.

G[0][y_, x_] := -I/k Exp[I k Abs[y - x]]
    
(* summed Dyson equation *)
G[n_][y_, x_] := G[n - 1][y, x] + V[n] G[n - 1][y, z[n]] G[n - 1][z[n], x]/
                 (1 - V[n] G[n - 1][z[n], z[n]])

We specialize to the case of equidistant potential spikes (with distance λ) of identical strength 𝒱 and cache togethered forms of the Green’s function. (Calculating these togethered forms instead of the direct ones is of course slower, but avoids an exponential increase of the size of the Green’s function.)

(* specialize to equal distance and identical strength *)
z[n_] := λ (n - 1); V[n_] := 𝒱
(* recursively cache calculated Green's functions *)
Do[G[n][y_, x_] = Together[G[n][y, x]], {n, 1, 8}]

This is the simplified form of the denominator of (Subscript[G, k]^((8)))(y,x;(V^((8)))). Its poles in the complex k-plane are the bound and resonant states of the potential. It factors into two parts.

den[k_, 𝒱_, λ_] = (If[Head[#] === Times, Collect[#, k, Factor]& /@ #, #]& @
                   Factor[Denominator[Together[G[8][0, 0]]]]) /. Abs[λ] -> λ

At k=0 , the denominator has a root of order 7 for all generic values of 𝒱 and λ.

Series[den[k, 𝒱, λ], {k, 0, 8}] // Factor

The following graphic shows the negative logarithm of the absolute value of the denominator over the i κ,λ-plane for attractive 𝒱=-1. Poles in the denominator show up as humps. For large λ (meaning the single delta potentials) are widely separated we have 8 bound states (κ<0) with exponential (in λ) separation. For smaller λ, some of the bound states merge into the “zero energy state” until only one bound state survives.

Plot3D[Evaluate[-Log @ Abs @ Re[den[I κ, -1, λ]]],
       {κ, -3, 3},  {λ, 0, 3}, PlotPoints -> 360, 
       Mesh -> False, BoxRatios -> {1, 1/2, 1/4}, 
       PlotRange -> All, ViewPoint -> {0, -3, 2.6},
       AxesLabel -> {"κ", "λ", None}]

The following graphic shows the logarithm of the absolute value of the denominator over the complex k-plane for 𝒱=i and λ = 1. Most of the resonance states lie near the real axis.

Plot3D[Evaluate[{Log @ Abs @ Re[den[kx + I ky, I, 1]],
                   Hue[(Cos[Arg[den[kx + I ky, I, 1]]] + 1)/2]}],
       {kx, -6, 6}, {ky, -2, 2},  PlotPoints -> 360 {1, 1/3}, 
       Mesh -> False, BoxRatios -> {1, 1/2, 1/4}, PlotRange -> All,
       ViewPoint -> {1, 3, 2}]

For a speedy generation of the animation frames, we generate a compiled version of the denominator. We use λ=1 in the following.

denCF = Compile[{{k, _Complex}, {𝒱, _Complex}, λ}, Evaluate[den[k, 𝒱, λ]]];

The function polePlot generates one frame of the animation for a given 𝒱=exp(i ω). The graphic shows the curves of vanishing real part in red and the curves of vanishing imaginary part in blue. The intersections are the bound and resonant states. We mark the poles by black dots. To find the intersection programmatically, we calculate the imaginary part of the denominators along the curves of vanishing real part. Because the intersections of the real and imaginary part are perpendicular (it is an analytic function in k), this will work fine. Using the endpoints of the line segments where the imaginary part changes sign as the starting points for FindRoot allows to calculate the poles. (We could have also calculated a differential equation k'(ω)=… by differentiating the denominator with respect to ω and solve this differential equation for once-calculated starting pole position for a fixed (Subscript[ω, 0]).)

polePlot[𝒱_, λ_, L_, pp_] := 
Module[{data, cpPlots, cps, segments, startPoints},
(* values of the denominator over the complex k-plane *)
data = Table[denCF[kx + I ky, 𝒱, λ], 
             {ky, -L/3, L/3, 2/3 L/(pp/3)}, {kx, -L, L, 2L/pp}];
(* contour plots of vanishing real and imaginary part *)
cpPlots = Show[cps = 
ListContourPlot[#1[data], Contours -> {0}, ContourShading -> False,
            MeshRange -> L {{-1, 1}, {-1, 1}/3}, AspectRatio -> Automatic,
            DisplayFunction -> Identity, ContourStyle -> {Hue[#2]}]& @@@
            {{Re, 0}, {Im, 0.76}}];
(* line segments of vanishing real part over the complex k-plane *)
segments = Map[{1, I}.#&, Flatten[Partition[First[#], 2, 1]& /@ 
                          Cases[Graphics[cps[[1]]], _Line, Infinity], 1], {-2}];
(* start points for the numerical root finding *)
startPoints = Plus @@@ segments[[Flatten[Position[Times @@@ 
                   Map[Sign[Im[denCF[#, 𝒱, λ]]]&, segments, {-1}], -1]]]]/2;
(* the poles *)
roots[𝒱] = Cases[FindRoot[den[k, 𝒱, λ], {k, #},
                     PrecisionGoal -> 5, AccuracyGoal -> 5]& /@ startPoints,
                   _Complex, {-1}];
(* show contour plots and poles *)
Show[{cpPlots, Graphics[{PointSize[0.01], GrayLevel[0], 
                    Point[{Re[#], Im[#]}]& /@ roots[𝒱]}]},
       DisplayFunction -> $DisplayFunction]]

Here are some of the resulting plots for various 𝒱.

polePlot[#, 1, 6, 360]& /@ {1, -I, -1}

And here is the animation. It shows how the 2×8 poles nearest to the origin take part in an exchange of their position with respect to the axis Re(k)=0. While doing this, the curves of vanishing real and imaginary parts perform a relatively complicated splitting and reuniting process. The other poles are largely uninfluenced by the changing potential.

frames = 121; 
Do[polePlot[Exp[I ω], 1, 6, 360], {ω, 0, 2Pi (1 - 1/frames), 2Pi/frames}]

The movement of the poles as the potential changed along the unit circle is shown in the following.

(* calculate zeros to be displayed in the next graphic *)
frames = 121; 
Do[Block[{$DisplayFunction = Identity},
          polePlot[Exp[I ω], 1, 6, 360]], 
         {ω, 0, 2Pi (1 - 1/frames), 2Pi/frames}]

Show[Graphics[{Hue[(Arg[#[[1, 1, 1]]] + Pi)/(2Pi)], 
               Point[{Re[#], Im[#]}]& /@ #[[2]]}& /@ DownValues[roots]],
     PlotRange -> All, AspectRatio -> Automatic]