Quantum Scattering in 1D

Mads Brandbyge, mbr@mic.dtu.dk
Mikroelectronik Centret (MIC)
Danish Technical University,  Build. 345 east
DK-2800 Lyngby, DENMARK

Some unimportant stuff - just run it once

Off[Inverse::"luc"]; 
    Off[PseudoInverse::"cfail"];
<< Graphics`FilledPlot`;
<< Graphics`Animation`;

Instruction:

You do the following:
[(0) Aswer "Yes" to initialize initial cells when asked in the beginning.]
(1) First build your potential landscape by writing a list of potential jumps as shown below.
(2) Then plot the potential.
(4) Evaluate "Calculate" and  see the Transmission function.
(5) Plot the scattering wavefunction at a certain energy, see the animation. 
(6) Change energy, change potential, play!

\(Define\ the\ potential\ in\ terms\ of\ "jumps" : \ {x\_\(\(jump\)\(\ \)\), \ V\_jump}\)

Define the part of the x-axis we want to consider: From L1 to L2

L1 = -5; L2 = 15;

potdiskont = {{0. , 1. }, {2, 0.1}, {4. , 1. }, {6, 0.1}, {8, 
          1. }, {10, 0. }}

Plot the potential

Clear[x]; 
    rules = {x > #[[1]], #[[2]]} & /@ Reverse[potdiskont]; 
    rules =  Append[rules, {1 == 1, 0. }]; 
    pot[x_] = Apply[Which[##] &, Flatten[rules]]; 
    potplot = 
      FilledPlot[pot[x], {x, L1, L2}, PlotPoints -> 50, 
        AxesLabel -> {"x", "V(x)"}, 
        Fills -> {{{1, Axis}, GrayLevel[ .8]}}]

Show[potplot]

Calculate

LocE[en_] = 
      Table[ If[
          i == 0, (en), (en - potdiskont[[i, 2]])], {i, 0, 
          Length[potdiskont]}]; K[en_] := Sqrt[2*LocE[en]] // N; 
    d[k1_, k2_, x_] =  .5*
        DiagonalMatrix[{Exp[-I*k1*x], 
              Exp[I*k1*x]}] . 
{{1 + (k2/k1), 
              1 - (k2/k1)}, {1 - (k2/k1), 
              1 + (k2/k1)}} . 
DiagonalMatrix[{Exp[
                I*k2*x], 
              Exp[-I*k2*
                  x]}];

Clear[TM]; Clear[TMList]; Clear[x];

\(\(getit[en_] := 
        Module[{TM, TMList, xLst, solution, B, F, CoeffList, FPList, r, t, R, 
            T, PsiIn, Psi}, 
TM = IdentityMatrix[2]; 

          TMList = {TM}; 

          TM = d[\ \ \(K[en]\)[\([1]\)], \(K[en]\)[\([2]\)], 
                potdiskont[\([1, 1]\)]\ \ \ ] // N; 
AppendTo[TMList, TM]; 

          For[i = 1, 
            i < Length[potdiskont], \(i++\), 
TM = 
              TM . d[\(K[en]\)[\([i + 1]\)], \(K[en]\)[\([i + 2]\)], \ 
                    potdiskont[\([i + 1, 1]\)]\ \ ] // 
                N; 
AppendTo[TMList, 
              TM];
]; 
Clear[xLst]; 
          xLst = \(Transpose[potdiskont]\)[\([1]\)]; 
Clear[B]; Clear[F]; 
          solution = Solve[TM . {F, 0} == {1, B}, {F, B}] // Flatten; \ 
          B = B /. solution; 

          invTMList = \(PseudoInverse[#] &\) /@ TMList; 

          CoeffList = 
            Table[Inverse[TMList[\([i]\)]] . {1, B}, {i, Length[TMList]}]; 

          FPList = 
            Table[{\ Exp[\ \ \ I*\((x*\(K[en]\)[\([i]\)])\)\ \ \ ]*
                      If[i == 1, 0, 
                        1], 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ Exp[\ \(-I\)*\((x*\(K[en]\)[\([i]\)]\ )\)\ \ ]} . 
                  CoeffList[\([i]\)], {i, 1, Length[K[en]]}] // 
              N; 

rules = \({x < \ #[\([1]\)], #[\([2]\)]} &\) /@ 
              Transpose[\ {Append[xLst, 10000], 
                  FPList}]; 

PsiIn[x_] = 
            If[x < xLst[\([1]\)], Exp[I*\((x\ *\(K[en]\)[\([1]\)])\)], 
              0. ]; 
Psi[x_] = 
            Apply[Which[##] &, 
              Flatten[rules, 1]]; 

r = 
            B /. solution; \ t = F /. solution; 
R = r*Conjugate[r] // Re; 

          T = t*Conjugate[t]*Re[\(K[en]\)[\([\(-1\)]\)]/\(K[en]\)[\([1]\)]] // 
              Re; 
Return[{R, T, PsiIn[x], 
              Psi[x]}];

];\)\)
R[en_] := getit[en][[1]]; T[en_] := getit[en][[2]]; 
    PsiIn[x] := getit[en][[3]]; 
    Psi[x] := getit[en][[4]];
\(\(Tplot = 
        Plot[T[en], {en, 0, 1.5}, PlotRange -> All, 
          PlotStyle -> {{RGBColor[0, 0, 0], 
                Thickness[0.01]}, {RGBColor[0, 0, 1], Thickness[0.01]}}, 
          AxesLabel -> {"E", "T(E)"}, 
          DefaultFont -> {"Times-Roman", 
              14}];\)


\)

Plot Transmission and Reflection for the potential

Show[Tplot, PlotRange -> {{0.3, 1.2}, All}]

Plot the wavefunction (time dependent) for a certain energy:

energy = 0.503;

Animation: After you have generated a time series select all images (outer tag) and press CLTR-Y for animation

Plot generation: The incoming wave is blue.

Clear[PP]; Clear[PPin];
PP[x_, t_] = 0.15*Re[getit[energy][[4]]*Exp[-I*2  Pi*t]]; 
    PPin[x_, t_] = 
      0.15*Re[getit[energy][[3]]*
            Exp[-I*2  Pi*t]];
tmpplot = 
      Table[FilledPlot[{pot[x], PP[x, t] + energy, PPin[x, t] + energy, 
            energy}, {x, L1, L2}, PlotRange -> {-0.5, 1.5}, 
          PlotPoints -> 50, AxesLabel -> {"x", None}, 
          Ticks -> {Automatic, None}, 
          Fills -> {{{1, Axis}, GrayLevel[ .8]}}, Curves -> Front, 
          PlotStyle -> {Thickness[0.012], RGBColor[1, 0, 0], 
              RGBColor[0, 0, 1], Thickness[0.005]}], {t, 0, 
          1.9,  .1}]

ShowAnimation[tmpplot]