Mach-Zehnder Interferometer

Zeitentwickelte Darstellung eines Wellenpakets in einem Potentialtopf mit Deltapotentialbarriere

von Thomas Lang

Initialisierung

ClearAll["Global`*"];
<< VQM`ArgColorPlot`;

Vorwort

Projektarbeit zur Vorlesung "Visuelle Quantenmechanik" SS 2001 bei Herrn Dr., Ao.Univ.-Prof., Univ.-Doz. f. Mathematik Bernd Thaller. 

Das Mach-Zehnder Interferometer basiert darauf das ein kolimierter Strahl durch einen Halbspiegel geteilt, an Vollspiegeln reflektiert und erneut an einem Halbspiegel wieder zusammengührt wird. Dabei tritt der zusammengeführte Strahl ausschließlich in der ursprünglichen Einfallsrichtung aus. Dieses Experiment kann auch mit einzelnen Teilchen bzw. Wellenpaketen durchgeführt werden, wobei scheinbar ebenso destruktive Interferenz in eine der beiden Ausgangsrichtungen auftritt was auch als "Single Particle Intereference" bezeichnet wird [ClEkHe99].



In diesem Notebook soll die Zeitentwicklung eines Wellenpakets dargestellt werden, dass an einem Potential gestreut und wieder zusammengeführt wird. Das Problem wird auf eine Dimension beschränkt gelöst.

Delta-Potential

Um einen Halbspiegel nachzustellen wird eine Potentialbarriere verwendet wobei dabei die Wahl auf ein Delta-Potential fällt, da es zahlreichen notwendigen Eigenschaften genügt. Ziel ist es ein Wellenpaket durch Streuung am Potential in zwei möglichst idente aber in Bewegungsruchtung spiegelverkehrte Teilpakete aufzuspalten. Eine beliebige Potentialbarriere würde die Impulskomponenten ungleich verteilen, da schnellere Anteile die Potentialbarriere wahrscheinlicher durchtunneln als langsamere Anteile. Da eine beliebige Potentialbarriere eine endliche Breite aufweist, bildet sich das Wellenpaket in dieser langsamer fort als außerhalb. Dadurch "hinkt" der transmittierte Anteil dem Reflektierten nach und die gewünschte Auslöschung wäre von unterschiedlichen Laufzeiten und daher Wegstrecken abhängig. 

Das Delta-Potential hingegen weißt antisymmetrische Reflexions- und Transmissionskoeffizienten auf - also eine antisymmetrische Impulverteilung wie gewünscht. Da es idealerweise nur punktbreit ist tretten keinerlei Laufzeitverzögerungen noch Mehrfachreflexionen auf.

opacity = 8;
FilledPlot[opacity^2/(k^2 + opacity^2), {k, 0, 100}, 
        Fills -> {GrayLevel[ .8]}, 
        Frame -> 
          True, 
PlotLabel -> "Reflexionsintensit\
ät", Axes -> {True, False}, 
        AxesLabel -> {"k", ""},  
        PlotRange -> {-0.03, 1}];
FilledPlot[k^2/(k^2 + opacity^2), {k, 0, 100}, 
        Fills -> {GrayLevel[ .8]}, 
        Frame -> 
          True, 
PlotLabel -> "Transmissionsintensit\
ät", Axes -> {True, False}, 
        AxesLabel -> {"k", ""},  
        PlotRange -> {-0.03, 1}];

Stationäre Lösung

Das Problem auf eine Dimension beschränkt entspricht einem unendlich hohen Potentialtopf mit einem Delta-Potential das leicht beschrieben werden kann durch

		V(
          x) = (ℏ/m) Ω δ(x)	und der vereinfachten Schrödingergleichung		u'' + [(k^2) - 
                2 Ω δ(x)]  u = 0.
		
Die Eigenfunktionen links und rechts vom Delta-Potential an der Stelle x=0 angesetzt, müssen unseren Randbedingen u(±a)=0 und natürlich auch den Stetigkeitsbedingungen der ersten Ableitung entsprechen 

		(Subscript[u, (x < 0)]):		(Subscript[A, 1]) Sin k(x+a) = (e^( i k x)) + 
        R (e^(-i k x))		-a<x≤0
		(Subscript[u, (x > 0)]):		(Subscript[A, 2]) Sin k(x-a) = (1 + 
            T) (e^( i k x))			0<x≤a
		
		(Subscript[(u'), (x < 0)]):		k (Subscript[A, 1]) Cos k(x+a) = \(TraditionalForm\`i\ k\ e\^\(\(\ \)\(i\ k\ x\)\) - i\ k\ R\ e\^\(\(-i\)\ k\ x\)\)	-a<x≤0
		(Subscript[(u'), (x > 0)]):		k (Subscript[A, 2]) Cos k(x-a) = i k (1 + 
            T) (e^( i k x))		0<x≤a

mit T als Transmissions- und R als Reflexionskoeffizient. Sowohl die Randbedingungen als auch die angesetzte Differentialgleichung sind invariant gegenüber dem Paritätsoperator dessen Eigenfunktionen also Lösungen der Gleichung sein müssen. Daher stellt der Sin(k(x±a)) eine ungerade Lösung dar, da er an den den Stellen x=±a und bei x=0 den Wert 0 hat und eine Verschwindene Ableitung. Die geraden Lösungen durch einsetzen und Elemenieren ergeben sich zu

		(Subscript[u, (x < 0)]):		A Sin k(x+a)		-a<x≤0
		(Subscript[u, (x > 0)]):		-A Sin k(x-a)		0<x≤a

wobei k-> (SuperPlus[(Subscript[k, n])]) mit Möglichen Werten (n-(1/2)) π < (SuperPlus[(Subscript[k, n])]) a < nπ  mit n ϵ Ν und den Koeffizienten der Eigenfunktionen 

		|A\( | \^2\) = ((2 (SuperPlus[(Subscript[k, n])]))/(2 (SuperPlus[(Subscript[k, n])]) a - Sin 2 (SuperPlus[(Subscript[k, n])]) a)) .

Integration der Schrödingergleichung über die Umgebung der Barriere führt zu

		u' (+0) - u' (-0) = 
        2  Ω u(0)		mit 	u(+0) = u(-0).

D.h. die Ableitung der geraden Eigenfunktionen hat eine Sprungstelle um 2Ω was innerhalb des Codes als opacity bezeichnet wird und der Durchlässigkeit der Potentialbarriere entspricht. Mit dem Ansatz R=T und Einsetzen der Ableitungen erhält man für Koeffizienten im Betragsquadrat

		|R\( | \^2\) = (((Ω^2))/((Ω^2) +  (k^2)))		|T\( | \^2\) = (((k^2))/((Ω^2) +  (k^2)))		|R\( | \^2\) = 1-|1 + T\( | \^2\)

Für unser Problem soll das Wellenpaket geteilt werden und suchen daher eine Lösung für ein Ω von (1/2):

Solve[k^2/(k^2 + opac) == 0.5, opac]

Mit den bekannten Beziehungen 

		E = ((p^2)/(2  m))  ->  (p^2) = 
          (E/(2  m))	vereinfacht mit m=1	und	(k^2)  =  2  E
		
wissen wir, dass der mittlere Impuls unsere Wellenpakets gleich Ω sein muß, damit eine Aufspaltung in zwei gleiche Teile gewährleistet ist.

Zeitentwickelte Lösung

Durch Projektion einer Gausschen Spektralfunktion auf die Eigenfunktionen wird ein Wellenpaket geformt. Dabei wird nur eine endliche Anzahl an Koeffizienten berechnet um den späteren Rechenaufwand zu beschränken. Zunächst die notwendigen Definitionen des Potentialtopfs, des mittleren Impulses u. Startposition des Wellenpakets und der Spektralfunktion. Anm.: alle Werte (besonders die Breite des Potentialtops) bestimmen die Genauigkeit und den Rechenaufwand der Darstellung!

a = 5;      
koeffNr = 30;     
p0 = 8;     
x0 = -2.5;     
opacity = p0;     

α = 1;     

g[x_] := (α/Pi)^(1/4)*
          Exp[-α*((x - x0)^2/2) + 
              I*p0*x];
\(\(\(ArgColorPlot[Evaluate[g[x]], {x, \(-a\), a}, 
        PlotRange -> {\(-0.1\), 1}, Axes -> {True, False}, 
        Frame -> 
          True, 
FrameTicks -> {{\(-5\), \(-2.5\), 0, 
              2.5, 5}, {0,  .25, 0.5, 0.75, 1}, None, 
            None}, 
PlotLabel -> 
          StringForm["Wellenpacket im Ortsraum"], 
Epilog -> {GrayLevel[ .8], 
            Line[{{0, 0}, {0, 
                  1}}]}
];\)\(
\)
    \)\)
gk[k_] := (1/Sqrt[2*Pi])*(1/α*Pi)^(1/4)*
          Exp[-(k - p0)^2/(2*α) + 
              I*x0*k];
ArgColorPlot[Evaluate[gk[k]], {k, 0, 2.5*a}, 
        PlotRange -> {-0.1, 1}, Axes -> {True, False}, 
        Frame -> 
          True, 
FrameTicks -> {{0, 2, 4, 6, 8, 10, 
              12, 14}, {0,  .25, 0.5, 0.75, 1}, None, 
            None}, 
PlotLabel -> 
          StringForm["Impulverteilung"]
];

Die Farbcodierung gibt die Phase des Wellenpakets wider und verdeutlicht hier seinen Impuls in positive x-Richtung. Weiterführende Informationen dazu bei [Th00].

Nun werden in einem Zug zunächst aus den Koeffizienten und Eigenwerten der geraden Lösungen die entsprechenden Eigenfunktionen gebildet.

knEven = {};
enEven = {};
cnEven = {};
\(\(\(AnEven = {};\)\(
\)
    \)\)
\(\(\(Table[
kn = 
          FindRoot[
            Cot[k*a]*k == \(-opacity\), {k, n*\((3/5)\) + 0.01*n}, 
            AccuracyGoal -> 8, MaxIterations -> 90, 
            WorkingPrecision -> \((2*$MachinePrecision)\)]; 
knEven = 
          Append[knEven, kn[\([1, 2]\)]]; 
enEven = 
          Append[enEven, kn[\([1, 2]\)]^2/2]; 
AnEven = 
          Append[AnEven, 
            Sqrt[Abs[\(2*kn[\([1, 2]\)]\)\/\(2*kn[\([1, 2]\)]*a - Sin[2*kn[\([1, 2]\)]*a]\)]]]; 
Print["n=", 
          n, \*""    \!\(\(k\_n\^+\)\)="", 
          N[kn[\([1, 2]\)], 4], \*""    \!\(\(A\_n\^+\)\)="", 
          N[AnEven[\([n]\)], 4], \*""    \!\(\(E\_n\^+\)\)="", 
          N[\((kn[\([1, 2]\)]^2)\)/2, 4]], 
{n, 1, 
          koeffNr}
];\)\(
\)
    \)\)
unEven[n_, x_] := 
        If[x < 0, 
          AnEven[[n]] Sin[
              knEven[[n]] (x + a)], -AnEven[[n]] Sin[
              knEven[[n]] (x - a)]];

Zur Veranschaulichung die ersten drei Eigenfunktionen unEven

border = 0.9;
xleft = -a - border;
xright = a + border;
lower = -0.5;
upper = -lower;

    
\(\(Table[
        Plot[unEven[n, x], {x, xleft, 
            xright - 0.1}, 
Axes -> {True, False}, 
          Frame -> True, 
          PlotRange -> {lower, 
              upper}, 
FrameTicks -> {{\(-5\), \(-2.5\), 0, 2.5, 5}, {}, None, None}, 
PlotLabel -> 
            StringJoin["n =" <> 
                ToString[n] <> \*"" at (Subscript[E, n])="" <> 
                ToString[
                  N[enEven[\([n]\)], 
                    2]]], 
Prolog -> {GrayLevel[ .8], 
              Line[{{0, lower}, {0, upper}}], 
              GrayLevel[
                0]}, 
Epilog -> {GrayLevel[ .7], 
Rectangle[{xleft, lower}, {\(-a\), 
                  99*upper/100}], 
Rectangle[{a, 
                  lower}, {xright, 
                  99*upper/100}], 
GrayLevel[
                0], 
Line[{{\(-a\), lower}, {\(-a\), 
                    upper}}], 
Line[{{a, lower}, {a, 
                    upper}}]}
], {n, 1, 3}];\)\)

Anmerkung: Da die Startstellen und Intervallbereiche des Newtonschen Verfahren von FindRoot nicht allgemein zu definieren sind, sollten sie bei Veränderung des ursprünglich gegebenen Werte überprüft werden! Dazu werden hier die Nullstellen und die zu lösende Funktion des Cotangens dargestellt.

gra = 
        Plot[Cot[k*a]*k + opacity, {k, 0, 4*a}];
grb = 
        ListPlot[Table[{knEven[[n]], 0}, {n, 1, koeffNr}], 
          PlotStyle -> {Hue[0.6], 
              PointSize[0.015]}];
Show[gra, grb];

Selbiges Verfahren für die eher unspektakulären Lösungen des einfachen Sinus - den ungeraden Lösungen unOdd

knOdd = {};
enOdd = {};
cnOdd = {};
\(\(\(AnOdd = Sqrt[Abs[1\/a]];\)\(
\)
    \)\)
\(\(\(Table[
kn = \(n\ π\)\/a; 
knOdd = Append[knOdd, kn]; 
enOdd = 
          Append[enOdd, kn^2/2]; 
Print["n=", 
          n, \*""    \!\(\(k\_n\^-\)\)="", 
          kn, \*""    \!\(\(A\_n\^-\)\)="", 
          AnOdd, \*""    \!\(\(E\_n\^-\)\)="", 
          kn^2/2], 
{n, 1, 
          koeffNr}
];\)\(
\)
    \)\)
unOdd[n_, x_] := 
        If[x < 0, -AnOdd*Sin[knOdd[[n]]*(x + a)], -AnOdd*
            Sin[knOdd[[n]]*(x - a)]];

    
\(\(Table[
        Plot[unOdd[n, x], {x, xleft, 
            xright - 0.1}, 
Axes -> {True, False}, 
          Frame -> True, 
          PlotRange -> {lower, 
              upper}, 
FrameTicks -> {{\(-5\), \(-2.5\), 0, 2.5, 5}, {}, None, None}, 
PlotLabel -> 
            StringJoin["n =" <> 
                ToString[n] <> \*"" at (Subscript[E, n])="" <> 
                ToString[
                  N[enOdd[\([n]\)], 
                    2]]], 
Prolog -> {GrayLevel[ .8], 
              Line[{{0, lower}, {0, upper}}], 
              GrayLevel[
                0]}, 
Epilog -> {GrayLevel[ .7], 
Rectangle[{xleft, lower}, {\(-a\), 
                  99*upper/100}], 
Rectangle[{a, 
                  lower}, {xright, 
                  99*upper/100}], 
GrayLevel[
                0], 
Line[{{\(-a\), lower}, {\(-a\), 
                    upper}}], 
Line[{{a, lower}, {a, 
                    upper}}]}
], {n, 1, 3}];\)\)

Die Spektralfunktion g auf die Eigenfunktionen unEven abgebildet liefert die nötigen Koeffizienten cnEven und 
auf die Eigenfunktionen unOdd abgebildet liefert die nötigen Koeffizienten cnOdd des Wellenpakets.

For[n = 1, 
        n <= koeffNr, n++, 
cn = 
          NIntegrate[unEven[n, x]*g[x], {x, -a, a}, 
            AccuracyGoal -> 8, 
            MaxRecursion -> 14]; 
cnEven = 
          Append[cnEven, cn];
];

    
For[n = 1, 
        n ≤ koeffNr, n++, 
cn = 
          NIntegrate[unOdd[n, x]*g[x], {x, -a, a}, 
            MaxRecursion -> 12, 
            AccuracyGoal -> 8]; 
cnOdd = 
          Append[cnOdd, cn];
];

Das Wellenpaket mit dem Zeitentwicklungsoperator (e^(-i E t))wird nun beschrieben durch  psi[x,t]

psi[x_, t_] := 
        Sum[cnEven[[n]]*unEven[n, x]*
              Exp[-I*enEven[[n]]*t]
 + 
            cnOdd[[n]]*unOdd[n, x]*Exp[-I*enOdd[[n]]*t], {n, 1, 
            koeffNr}];

Die Zeitenwicklung wird in folgender Zelle evaluiert. Bei den ursprünglichen Angaben betrug die Rechenzeit auf einem PIII-450Mhz rund 1.5 Stunden - die Schrittweite dt sollte daher angepasst werden. Ebenso kann man statt der Funktion ArgColorPlot das Standart Plot verwenden, muß aber den Betrag von psi[x,t] darstellen lassen.

t0 = 0. ;   
t1 = 
        2.25;    
dt = 0.025;    

  
border = 0.9;
xleft = -a - border;
xright = a + border;
lower = -0.1;
\(\(\(upper = 1;\)\(
\)
    \)\)

Animation

Do[ 
  QArgColorPlot[Evaluate[psi[x,t]],{x,xleft,xright-0.1},
    Axes->{True,False},Frame->True,PlotRange->{lower,upper},QSquared->True,
    FrameTicks->{{-5,-2.5,0,2.5,5},{0,.25,0.5,0.75,1},None,None},
    PlotLabel->StringForm["t =``",PaddedForm[N[t],{3,2}]],
    Prolog->{GrayLevel[.8],Line[{{0,0},{0,upper}}],GrayLevel[0]},
    Epilog->{GrayLevel[.7],
        Rectangle[{xleft,lower},{-a,99*upper/100}],
        Rectangle[{a,lower},{xright,99*upper/100}],
        GrayLevel[0],
        Line[{{-a,lower},{-a,upper}}],
        Line[{{a,lower},{a,upper}}]}
    ]
  ,{t,t0,t1,dt}];

(* output cells deleted *)

Diskussion

Der erwartete Effekt der Auschlöschung ist deutlich zur Zeit t= 2.2 zu erkennen, wenn auch es nicht ganz perfekt klappt - ein Effekt der ohne das experimentelle Vorwissen nicht zu erwarten gewesen wäre (zumindest meinerseits ; )). 

Variert man die Durchlässigkeit des Delta-Potential ein wenig läßt sich die Auslöschung noch erahnen. Im Extremfall wirkt das DeltaPotential wie ein unedlichhohe Potentialwand oder vollkommen durchlässig.

Fehlerquellen zur unvollständigen Auslöschung:

 -  Eine Steigerung der Koeffizienten verstärkt die Auslöschung unwesentlich - hinzukommt die schwere Kontrolle der Nullstellen wie bereits bemerkt. 

 -  Da sich das Wellenpaket mit der Zeit zunehmend verbreitert und schon von Anfang an die Grenzen des Potentialtops berührt trägt es natürlich zu den Restbeträgen bei. Eine andere Wahl des oben verwendeten Faktros α verschmälert das Paket zunächst, bedingt aber eine größere Impulsunschärfe und damit eine schnellere Verbreiterung. Zusätzliche Skalierungen in beiden Darstellungräumen brachten ebenfalss keine deutlichen Verbesserungen. 
 
-  Die Rechengenauigkeit beläuft sich durchwegs auf 8 Stellen nach dem Komma und scheidet zunächst als Ursache aus, sofern die Darstellung nicht detailgetreuer festgelegt wird.

Die, auf die hier verzichtete Impulsraumdarstellung, lieferte keine Besonderheiten, sondern den schon bekannten Wechsel der anfänglichen Impulse vom Zentrum um +(Subscript[p, 0]) nach -(Subscript[p, 0]) und umgekehrt.

Eine fortführende Arbeit am Problem wäre die Abhängigkeit der Auslöschung von der Laufzeit der Teilpakete zu visualieren:
 -  durch unterschiedliche Wegstrecken der beiden Seiten umd das Potential
 -  durch eine zusätzlichen Potentialstufe ähnlich dem i.d. LV vorgestellten Dekohärenzexperiment mit Neutronen

Literaturverzeichnis

[Th00]   Thaller Bernd, Visual Quantum Mechanics
	(Springer Verlag, August 2000)
	ISBN: 0387989293
	Bernd Thaller
	
[Fl98]   Siegfried Flügge, Practical Quantum Mechanics
	(Springer Verlag, Dezember 1998)
	ISBN: 3540650350

[ClEkHe99]    R.Cleve, A.Ekert, L. Henderson, On Quantum Algorithms
	(Preprints-archive quant-ph/9903061)