130B Quantum Physics

Part 2. Piecewise Potentials

Basic Piecewise Potentials

General Setup

Schrödinger Equation

(OverHat[H])=(((OverHat[p])^2)/(2m))+V((OverHat[x])).

(OverHat[x]): position operator,

(OverHat[p])=-i ℏ(∂_x): momentum operator.

(OverHat[H]) ψ(x)=E ψ(x),

This is an eigen equation for the Hamiltonian operator (OverHat[H]).

Each solution ψ(x) is an energy eigenstate with eigen energy E.

Stationary State

ψ(x,t)=ψ(x)(e^(-i E t/ℏ)).

The time dependence is only a global phase rotation.

The probability density ρ(x) and the probability current density j(x) both remain time-independent:

ρ(x)=((|ψ(x,t)|)^2)=((|ψ(x)|)^2),
j(x)=(1/(2m))((SuperStar[ψ(x)])(OverHat[p]) ψ(x)-ψ(x) (OverHat[p]) (SuperStar[ψ(x)])).

For any physical observable (OverHat[O]) that does not explicitly depend on time,

its expectation value ⟨O⟩ will be invariant on stationary state,

its measurement outcome probability distribution will be invariant on stationary state.

Stationary states are energy eigenstates.

Free Particle (Flat Potential)

Wave Function

(OverHat[H])=(((OverHat[p])^2)/(2m))=-((ℏ^2)/(2m))(∂_x^2).

Energy Spectrum

Energy spectrum is continuous. (∵ k takes continuous value, ∴ E is also continuous)

Level degeneracy = 2 (except for the k=0 state).

Travelling and Standing Wave

Travelling wave (momentum eigenstate, ψ(x)=(e^(i k x)))

k>0: right-moving, i.e. k=+(ℏ^(-1))(Sqrt[2m E]).

k<0: left-moving, i.e. k=-(ℏ^(-1))(Sqrt[2m E]).

They are orthogonal states, related by the reflection symmetry x->-x.

If two states are related by symmetry, their energy must be the same. => They are degenerated.

For 1D free particle with reflection symmetry, every energy level is two-fold degenerated (except for k=0).

Sanding wave (equal weight superposition of travelling wave)

even parity:(1/2)((e^(i k x))+(e^(-i k x)))=cos k x,
odd parity:(1/(2i))((e^(i k x))-(e^(-i k x)))=sin k x.

Linear superposition of degenerated eigenstates are still eigenstates.

Generic superposition of travelling waves

Ψ(x,t)=ψ(x)(e^(-i ω t)),
ψ(x)=A (e^(i k x))+B (e^(-i k x)),

Normalization condition: ((|A|)^2)+((|B|)^2)=1.

DynamicModule[{A=1,B=0,t=0},
Column[Row/@{{Style["\!\(A:\) ","InlineFormula"],DynamicWrapper[Manipulator[Dynamic[A],{0,1},ImageSize->180],B=Sqrt[1-A^2]]},
{Style["\!\(B:\) ","InlineFormula"],DynamicWrapper[Manipulator[Dynamic[B],{0,1},ImageSize->180],A=Sqrt[1-B^2]]},
{Style["\!\(t:\) ","InlineFormula"],Animator[Dynamic[t],{0,1},AnimationRate->0.2,AnimationRunning->False,ImageSize->{200,16}]},
{Dynamic[Plot[Evaluate@ReIm[(A Exp[2π I x]+B Exp[-2π I x])Exp[-2π I t]],{x,0,3},FrameLabel->{"\!\( x\)","\!\(ψ(x)\)"},
PlotStyle->{Blue,Yellow},PlotLegends->(Style[#,"Text"]&/@{"Re","Im"}),PlotRange->{-1.5,1.5},ImageSize->200]]}}]]

Wave node and antinode.

Probability distribution

ρ(x)=((|ψ(x)|)^2)=((|A (e^(i k x))+B (e^(-i k x))|)^2).

DynamicModule[{A=1,B=0},
Column[Row/@{{Style["\!\(A:\) ","InlineFormula"],DynamicWrapper[Manipulator[Dynamic[A],{0,1},ImageSize->180],B=Sqrt[1-A^2]]},
{Style["\!\(B:\) ","InlineFormula"],DynamicWrapper[Manipulator[Dynamic[B],{0,1},ImageSize->180],A=Sqrt[1-B^2]]},
{Dynamic[Plot[Abs[(A Exp[2π I x]+B Exp[-2π I x])]^2,{x,0,3},FrameLabel->{"\!\( x\)","\!\(ρ(x)\)"},Filling->Axis,FillingStyle->LightGray,PlotStyle->Black,PlotRange->{0,2},PlotRangePadding->0.1,ImageSize->200]]}}]]

Hard-Wall Potential

Wave Function

V(x)=({{{, ({{0, x>0,}, {∞, x<0.}})}})

Particle can not enter the negative half-space (x<0) => no probability to appear in the negative half-space => wavefunction must vanish for x≤0 => the hard-wall condition

ψ(x≤0)=0.

In the positive half-space (x>0), Hamiltonian still looks like Eq. () => general solution is given by Eq. ()

ψ(x≥0)=A (e^(i k x))+B (e^(-i k x)).

Energy Spectrum

Energy spectrum is still continuous. E=(ℏ^2)(k^2)/(2m)

Level degeneracy = 1 => only sin k x states survives. cos k x states are not consistent with the hard-wall condition ψ(x≤0)=0.

WaveExplorer[LeftWall[x],{x,-1,10},3,{-0.5,5}]

Infinite Square Well

Wave Function

V(x)=({{{, ({{0, 0<x<a,}, {∞, elsewhere.}})}})

ψ(x=0)=ψ(x=a)=0.

Energy Spectrum

Energy spectrum becomes discrete

n labels the energy level, so it is also called the level index.

Level degeneracy = 1 (non-degenerated).

WaveExplorer[LeftWall[x]+RightWall[x-1],{x,-0.1,1.1},5,{-5,60},GridLines->{None,π^2/2Range[1,3]^2}]

As we squeeze the potential (a smaller), energy will increase. => an example of the uncertainty principle. The particle exerts a pressure to the walls.

Probability Distribution

ρ(x)=({{{, ({{(2/a)(sin^2) ((n π)/a) x, 0<x<a,}, {0, elsewhere.}})}})

Probability must add up to 1 => normalization condition:

Check the normalization condition Eq. ().

Simplify[Integrate[(2/a)Sin[n π x/a]^2,{x,0,a}],n∈Integers]

The prefactor in the wavefunction is called the normalization factor.

Bound state wavefunctions are normalizable.

Time Evolution of Generic States

The initial state can be decomposed on the energy eigen basis.

ψ(x,0)=(∑_n)(Subscript[c, n])(Subscript[ψ, n])(x).

Coefficient is given by the inner product (calculated as integral)

(Subscript[c, n])=∫dx (Subscript[ψ, n]^( * ))(x)ψ(x,0).

Each energy eigen state evolves independently (phase angle rotates according to the energy)

ψ(x,t)=(∑_n)(Subscript[c, n])(Subscript[ψ, n])(x)(e^(-i (Subscript[E, n])t)).

DynamicModule[{w=0.5,xc=0.1,nmax=10,cs,gc},
cs=Normalize@Table[NIntegrate[Sin[n π x]x(1-x)Exp[-((x-xc)/w)^2],{x,0,1},AccuracyGoal->5],{n,1,nmax}];
gc=ListPlot[Abs[cs]^2,PlotStyle->Black,Filling->Axis,FillingStyle->Black,PlotRange->{0,1},PlotRangePadding->Scaled[0.05],FrameLabel->{"\!\( n\)","\!\(\[LeftBracketingBar]c\_n\[RightBracketingBar]\^2\)"},ImagePadding->{{44, 2}, {40, 5}},ImageSize->180];
Manipulate[Column@{gc,Plot[Abs[cs.Table[Sqrt[2]Sin[n π x]Exp[-I 2 π n^2 t],{n,1,nmax}]]^2,{x,0,1},FrameLabel->{"\!\( x\)","\!\(ρ(x)\)"},PlotStyle->Black,Filling->Axis,PlotRange->{0,7},PlotRangePadding->{0.04,0.2},ImagePadding->{{44, 2}, {40, 5}},ImageSize->180]},{t,0,1,ImageSize->Small,AnimationRate->0.1},Paneled->False]]

Step Potential

Wave Function

V(x)=(Subscript[V, 0]) Θ(x)=({{{, ({{0, x<0,}, {(Subscript[V, 0]), x>0.}})}})

E=((p^2)/(2m))+V(x)=(((ℏ^2)(k^2))/(2m))+({{{, ({{0, x<0,}, {(Subscript[V, 0]), x>0.}})}})

(Subscript[k, 1]^2)=((2m E)/(ℏ^2)), (Subscript[k, 2]^2)=((2m (E-(Subscript[V, 0])))/(ℏ^2)).

ψ(x)=({{{, ({{(e^(i (Subscript[k, 1]) x))+r (e^(-i (Subscript[k, 1])x)), x<0,}, {t (e^(i (Subscript[k, 2])x)), x>0.}})}})

r=(((Subscript[k, 2])-(Subscript[k, 1]))/((Subscript[k, 1])+(Subscript[k, 2]))), t=((2(Subscript[k, 1]))/((Subscript[k, 1])+(Subscript[k, 2]))).

Scattering State (E>(Subscript[V, 0]))

Travelling wave

Level degeneracy = 2 (left- and right-moving)

WaveExplorer[UnitStep[x],{x,-5,5},2,{-0.2,3}]

Reflection coefficient

Transmission coefficient

Conservation of probability implies that R+T=1.

Total Reflection State (0<E<(Subscript[V, 0]))

Standing wave

Level degeneracy = 1

WaveExplorer[UnitStep[x],{x,-5,5},0.5,{-0.2,3}]

ψ(x)=({{{, ({{(e^(i (Subscript[k, 1]) x))+r (e^(-i (Subscript[k, 1])x)), x<0,}, {t (e^(-(Subscript[κ, 2])x)), x>0.}})}})

The particle can penetrate into the potential barrier, but is wave amplitude decays exponentially. Penetration depth (a length scale)

ξ=(1/(Subscript[κ, 2]))=(ℏ/(Sqrt[2m ((Subscript[V, 0])-E)])).

Solution for amplitudes r and t

r=((i (Subscript[κ, 2])-(Subscript[k, 1]))/((Subscript[k, 1])+i (Subscript[κ, 2]))), t=((2(Subscript[k, 1]))/((Subscript[k, 1])+i (Subscript[κ, 2]))).

Phase shift of the reflection wave

δ=arg r=2arccot ((Subscript[κ, 2])/(Subscript[k, 1])).

Reflection and transmission coefficients

Total reflection (no transmission).

Square Potential

General Discussion

V(x)=({{\[Piecewise], ({{(Subscript[V, 0]), \[LeftBracketingBar]x\[RightBracketingBar]<a/2,}, {0, \[LeftBracketingBar]x\[RightBracketingBar]>a/2.}})}})

E=(((ℏ^2)(k^2))/(2m))+({{\[Piecewise], ({{(Subscript[V, 0]), \[LeftBracketingBar]x\[RightBracketingBar]<a/2,}, {0, \[LeftBracketingBar]x\[RightBracketingBar]>a/2.}})}})

(Subscript[k, 1]^2)=((2m E)/(ℏ^2)), (Subscript[k, 2]^2)=((2m (E-(Subscript[V, 0])))/(ℏ^2)),

(Subscript[k, 1]^2)-(Subscript[k, 2]^2)=((2m (Subscript[V, 0]))/(ℏ^2)).

The potential is symmetric under reflection: V(x)=V(-x), so every eigenstate has a definite parity quantum number (this is a (Subscript[ℤ, 2]) symmetry).

Even parity ψ(-x)=ψ(x)

(Subscript[ψ, e])(x)=({{{, ({{(Subscript[A, e])(e^(i (Subscript[k, 1])x))+(Subscript[B, e])(e^(-i (Subscript[k, 1])x)), x>a/2}, {cos (Subscript[k, 2])x, |x|<a/2}, {(Subscript[A, e])(e^(-i (Subscript[k, 1])x))+(Subscript[B, e])(e^(i (Subscript[k, 1])x)), x<-a/2}})}}).

Odd parity ψ(-x)=-ψ(x)

(Subscript[ψ, o])(x)=({{{, ({{(Subscript[A, o])(e^(i (Subscript[k, 1])x))+(Subscript[B, o])(e^(-i (Subscript[k, 1])x)), x>a/2}, {sin (Subscript[k, 2])x, |x|<a/2}, {-(Subscript[A, o])(e^(-i (Subscript[k, 1])x))-(Subscript[B, o])(e^(i (Subscript[k, 1])x)), x<-a/2}})}}).

Bound State in Potential Well ((Subscript[V, 0])<E<0)

The Logarithmic Trick

Even parity ψ(-x)=ψ(x)

(Subscript[ψ, e])(x)=({{{, ({{(Subscript[A, e])(e^(-(Subscript[κ, 1])x)), x>a/2}, {cos (Subscript[k, 2])x, |x|<a/2}, {(Subscript[A, e])(e^((Subscript[κ, 1])x)), x<-a/2}})}}).

Odd parity ψ(-x)=-ψ(x)

(Subscript[ψ, o])(x)=({{{, ({{(Subscript[A, o])(e^(-(Subscript[κ, 1])x)), x>a/2}, {sin (Subscript[k, 2])x, |x|<a/2}, {-(Subscript[A, o])(e^((Subscript[κ, 1])x)), x<-a/2}})}}).

(∂_x)ln (Subscript[ψ, e])(a/2+(SubPlus[0]))=(∂_x)ln (Subscript[ψ, e])(a/2-(SubPlus[0])),
(∂_x)ln (Subscript[ψ, o])(a/2+(SubPlus[0]))=(∂_x)ln (Subscript[ψ, o])(a/2-(SubPlus[0])).

(Subscript[A, e])=(e^((Subscript[κ, 1])a/2))cos (((Subscript[k, 2])a)/2), (Subscript[A, o])=(e^((Subscript[κ, 1])a/2))sin (((Subscript[k, 2])a)/2).

Pictorial Solutions

DynamicModule[{kmax=4π,plt,lbl},
plt=Plot[{k2 Tan[k2/2],-k2 Cot[k2/2]},{k2,0,kmax},
PlotRange->{{0,kmax},{0,kmax}},GridLines->{Range[0,kmax,π],None},FrameTicks->{{Range[0,kmax,π],None},{Range[0,kmax,π],None}},
FrameLabel->{"\!\(k\_2 a\)","\!\(κ\_1 a\)"},AspectRatio->1,ImageSize->160];
lbl=Grid[{Graphics[{#1,Line[{{0,0},{1,0}}]},AspectRatio->0.2,ImageSize->15],Style[#2,"Text"]}&@@@{{Blue,"even"},{Red,"odd"}}];
Manipulate[Row[{Show[{plt,Graphics[{Circle[{0,0},ka,{0,π/2}]}]}],lbl}],{{ka,2π,Style["\!\(k a\)","InlineFormula"]},0,kmax,ImageSize->Small},Paneled->False],SaveDefinitions->True]

The k a->0 limit (meaning (Subscript[V, 0])(a^2)->0), physically it can be

a shallow well: (Subscript[V, 0])->0 with a fixed,

a narrow well: a->0 with (Subscript[V, 0]) fixed,

a sharp well: a->0, (Subscript[V, 0])->-∞, with (Subscript[V, 0])a=γ fixed.

In either cases, there is one (and only one) bound state solution remaining in the even parity sector. All the higher energy bound states are pushed out of the well.

WaveExplorer[-10UnitStep[1/2-Abs[x]],{x,-1,1},-5,{-10.5,1}]

Bound state energy levels are discrete and non-degenerated

(Subscript[E, n])=-((ℏ^2)/(2m))(Subscript[κ, 1,n]^2).

The number of bound states are given by

(Subscript[n, max])=⌊(((((2m (Subscript[V, 0])))^(1/2))a)/(π ℏ))⌋+1.

In the sharp well limit (a->0 with (Subscript[V, 0])a=γ fixed), the bound state is unique and its energy is given by

(Subscript[E, 1])=-(m/(2(ℏ^2)))((((Subscript[V, 0])a))^2)=-((m (γ^2))/(2(ℏ^2))).

k a∝((|(Subscript[V, 0])|)^(1/2))a=((|γ|)^(1/2))(a^(1/2))->0.

((Subscript[κ, 1])a)=((Subscript[k, 2])a) tan (((Subscript[k, 2])a)/2)->(1/2)((((Subscript[k, 2])a))^2),
(((k a))^2)=((((Subscript[κ, 1])a))^2)+((((Subscript[k, 2])a))^2)=((((Subscript[κ, 1])a))^2)+2((Subscript[κ, 1])a).

(Subscript[E, 1])=-((ℏ^2)/(2m))(Subscript[κ, 1]^2)=-((ℏ^2)/(2m))(((k^4)(a^2))/4)=-((ℏ^2)/(2m))(1/4)(((-((2m (Subscript[V, 0]))/(ℏ^2))))^2)(a^2)=-(m/(2(ℏ^2)))((((Subscript[V, 0])a))^2).

Scattering State (E>0 regardless of (Subscript[V, 0]))

The energy levels should be continuous for scattering states, so any E>0 is admissible (there is no need to solve for (Subscript[k, 1]) and (Subscript[k, 2])).

Every level is doubly degenerated (degeneracy = 2), corresponding to (Subscript[ψ, e]) and (Subscript[ψ, o]).

ψ(x)=({{{, ({{t (e^(i (Subscript[k, 1])x)), x>a/2}, {…, |x|<a/2}, {r (e^(-i (Subscript[k, 1])x))+(e^(i (Subscript[k, 1])x)), x<-a/2}})}}),

r=(1/2)(((Subscript[A, e])/(Subscript[B, e]))+((Subscript[A, o])/(Subscript[B, o]))),
t=(1/2)(((Subscript[A, e])/(Subscript[B, e]))-((Subscript[A, o])/(Subscript[B, o]))).

The transmission coefficient

Derive Eq. ().

Module[{ψ1,ψ2e,ψ2o,r,t,invT},
ψ1[x_]:=A Exp[I k1 x]+B Exp[-I k1 x];
ψ2e[x_]:=Cos[k2 x];
ψ2o[x_]:=Sin[k2 x];
{{r},{t}}=Simplify[({{1,1},{1,-1}}/2).((A/B)/.Table[Solve[{ψ1[a/2]==ψ2[a/2],ψ1'[a/2]==ψ2'[a/2]},{A,B}],{ψ2,{ψ2e,ψ2o}}])];
invT=Simplify@ComplexExpand[1/Abs[t]^2];
1/(Simplify[Simplify[Numerator@invT,k1^2==2 m Ε&&k2^2==2m(Ε-V0)]/(Denominator@invT/.{k1^2->2 m Ε,k2^2->2m(Ε-V0)})])]

where (Subscript[k, 2]) is set by (E-(Subscript[V, 0])) according to Eq. ().

In the sharp well limit (a->0 with (Subscript[V, 0])a=γ fixed), Eq. () becomes

T=(((1+((m (γ^2))/(2(ℏ^2)E))))^(-1)).

Prove Eq. ().

(sin^2) (Subscript[k, 2])a→(Subscript[k, 2]^2)(a^2)=((2m (E-(Subscript[V, 0])))/(ℏ^2))(a^2),

T=(((1+((Subscript[V, 0]^2)/(4E(E-(Subscript[V, 0]))))((2m (E-(Subscript[V, 0])))/(ℏ^2))(a^2)))^(-1))
=(((1+((m (Subscript[V, 0]^2)(a^2))/(2(ℏ^2)E))))^(-1))
=(((1+((m (γ^2))/(2(ℏ^2)E))))^(-1)).

Scattered by a Potential Well ((Subscript[V, 0])<0)

WaveExplorer[-10UnitStep[1/2-Abs[x]],{x,-4,4},5,{-10.5,10.5}]

Resonant transmission: perfect transmission due the the cancelation (by interference) of reflection waves form the front and back interfaces.

Manipulate[
Module[{V0=-10^lgV0},Plot[{(1+V0^2/(4f Abs[V0](f Abs[V0]-V0))Sin[Sqrt[2(f Abs[V0]-V0)]]^2)^(-1),(4f(f-Sign[V0]))/(2f-Sign[V0])^2},{f,0,2},FrameLabel->{"\!\(E/\[LeftBracketingBar]V\_0\[RightBracketingBar]\)","\!\(T \)"},PlotStyle->{Black,Directive[AbsoluteThickness[0.8],Gray,Dashed]},PlotLabel->"\!\(V\_0="<>ToString[V0]<>"\)",PlotRange->{0,1},ImageSize->180]],{{lgV0,0,Style["\!\(log \[LeftBracketingBar]V\_0\[RightBracketingBar]\)","InlineFormula"]},-3,4,ImageSize->Small},SaveDefinitions->True,Paneled->False]

Resonance happens at (Subscript[k, 2])a=n π (n∈ℤ), where sin (Subscript[k, 2])a=0 and hence T=1.

Energy associated to the resonance forms discrete level

(Subscript[E, n])=(((π^2)(ℏ^2))/(2m (a^2)))(n^2) +(Subscript[V, 0]),

for n≥(Subscript[n, max]), where (Subscript[n, max]) is the number of bound states given in Eq. (). Compare with Eq. (), these levels locate at the same energy as the bound state level of a infinite depth well of the same width.

Lower bound of the transmission coefficient is set by

(Subscript[T, low])=((4E(E-(Subscript[V, 0])))/(((2E-(Subscript[V, 0])))^2)).

Scattered by a Potential Barrier ((Subscript[V, 0])>0)

WaveExplorer[10UnitStep[1/2-Abs[x]],{x,-4,4},5,{-10.5,10.5}]

Quantum tunnelling: Even if the energy is bellow the barrier (0<E<(Subscript[V, 0])), wave can still tunnel through.

But the tunneling rate is exponentially suppressed.

T=((16 E((Subscript[V, 0])-E))/(Subscript[V, 0]^2))exp(-((2a)/ℏ)(Sqrt[2m ((Subscript[V, 0])-E)])).

Heavier particles (larger m) are harder to tunnel.

The wider and higher the barrier is, the harder to tunnel.

Manipulate[
Module[{V0=10^lgV0},Plot[{(1+V0^2/(4f Abs[V0](f Abs[V0]-V0))Sin[Sqrt[2(f Abs[V0]-V0)]]^2)^(-1),(4f(f-Sign[V0]))/(2f-Sign[V0])^2},{f,0,2},FrameLabel->{"\!\(E/\[LeftBracketingBar]V\_0\[RightBracketingBar]\)","\!\(T \)"},PlotStyle->{Black,Directive[AbsoluteThickness[0.8],Gray,Dashed]},PlotLabel->"\!\(V\_0="<>ToString[V0]<>"\)",PlotRange->{0,1},ImageSize->180]],{{lgV0,0,Style["\!\(log \[LeftBracketingBar]V\_0\[RightBracketingBar]\)","InlineFormula"]},-3,4,ImageSize->Small},SaveDefinitions->True,Paneled->False]

Above the barrier (E>(Subscript[V, 0])): resonant tunneling happens at

(Subscript[E, n])=(((π^2)(ℏ^2))/(2m (a^2)))(n^2) +(Subscript[V, 0]),

for n=1,2,3,….

Dirac Potential

General Discussion

V(x)=γ δ(x).

(∂_x^2) ψ(x)=((2m)/(ℏ^2)) (γ δ(x)-E) ψ(x).

(∂_x) ψ((SubPlus[0]))-(∂_x) ψ((SubMinus[0]))=((2m)/(ℏ^2))(∫_(SubMinus[0])^(SubPlus[0]))dx (γ δ(x)-E) ψ(x)=((2m γ)/(ℏ^2))ψ(0).

Matching Condition: Eq. () requires the first-order derivative to jump by

Continuity Condition: Wavefunction must be continuous

Bound State (γ<0, E<0)

Attractive impurity traps bound states.

WaveExplorer[-2DiracDelta[x],{x,-2,2},-2.1,{-4,2},"Scale"->0.2]

There is a unique bound state trapped by the Dirac potential

Note that because γ<0, the scattering length (Subscript[a, s]) is negative.

The bound state is of even parity. There is no odd parity state that can be trapped by Dirac potential.

The corresponding eigen energy is

which is consistent with Eq. ().

Scattering State (E>0)

Both attractive and repulsive impurity scatters the particle around.

ψ(x)=({{{, ({{(e^(i k x))+r (e^(-i k x)), x<0,}, {t (e^(i k x)), x>0,}})}})

r=-(1/(1-i k (Subscript[a, s]))),
t=(1/(1+i /(k (Subscript[a, s])))).

WaveExplorer[-5DiracDelta[x],{x,-4,4},5,{-5,10.5}]

WaveExplorer[5DiracDelta[x],{x,-4,4},5,{-2,10.5}]

The reflection and transmission coefficients

They sum up to unity R+T=1. This is consistent with Eq. () for the sharp square potential. The reflection is enhanced (or the transmission is reduced)if

the Dirac potential is strong (|γ| large),

the energy is low (E small),

the particle is massive (m large). Intuitive argument: heavy particle => large momentum => short wave length => hard to diffract through the impurity, scatters strongly.

Scattering phase shift: the reflected and transmitted wave are phase shifted from the incident wave.

Weak potential limit (γ->0, (Subscript[a, s])->∞): (Subscript[Δφ, t])=0, plane wave transmitted through the potential without phase shift (as if nothing happened).

Strong potential (γ->∞, (Subscript[a, s])->0): (Subscript[Δφ, r])=π, all incident wave are reflected with half-wave loss.

Double Dirac Potential

General Discussion

V(x)=γ (δ(x-a/2)+δ(x+a/2)).

Define momentum k and scattering length (Subscript[a, s])

(k^2)=((2m E)/(ℏ^2)), (Subscript[a, s])=((ℏ^2)/(m γ)).

Bound state (E<0) wavefunction (k=i κ in this case)

Even parity

ψ(x)=({{{, ({{A (e^(κ x)), x<-a/2,}, {cosh κ x, -a/2<x<a/2,}, {A (e^(-κ x)), x>a/2.}})}})

Odd parity

ψ(x)=({{{, ({{-A (e^(κ x)), x<-a/2,}, {sinh κ x, -a/2<x<a/2,}, {A (e^(-κ x)), x>a/2.}})}})

Scattering state (E>0) wavefunction

ψ(x)=({{{, ({{(e^(i k x))+r (e^(-i k x)), x<-a/2,}, {A (e^(i k x))+B (e^(-i k x)), -a/2<x<a/2,}, {t (e^(i k x)), x>a/2.}})}})

Bound State (E<0)

Depending on the scattering length (Subscript[a, s])=(((ℏ^2)(m γ))^(-1)):

(Subscript[a, s])>0: no bound state,

-a<(Subscript[a, s])<0: one bound state of even parity,

(Subscript[a, s])<-a<0: two bound states of both parities.

WaveExplorer[-2(DiracDelta[x-0.5]+DiracDelta[x+0.5]),{x,-2,2},-2.5,{-3,1}]

By matching wavefunction, κ can be determined from

Manipulate[Plot[Exp[κ a](1-κ),{κ,0,2},PlotStyle->Black,GridLines->{{1,2,If[a>1,1-1/a,Nothing],{If[a>0,1+ProductLog[a Exp[-a]]/a,2],Directive[AbsoluteThickness[1],Opacity[1],Blue]},If[a>1,{1+ProductLog[-a Exp[-a]]/a,Directive[AbsoluteThickness[1],Opacity[1],Red]},Nothing]},{{-1,Directive[AbsoluteThickness[1],Opacity[1],Blue]},{1,Directive[AbsoluteThickness[1],Opacity[1],Red]}}},FrameLabel->{"\!\(κ \[LeftBracketingBar]a\_s\[RightBracketingBar]\)","\!\(E\^\(κ a\)(1+κ a\_s)\)"},PlotLabel->"\!\(a/\[LeftBracketingBar]a\_s\[RightBracketingBar]="<>ToString[a]<>"\)",PlotRange->{-3,3},ImageSize->200],{{a,1,Style["\!\( a:\)","InlineFormula"]},0,5,ImageSize->Small},Paneled->False]

Even parity: , odd parity: .

When a<\[LeftBracketingBar](Subscript[a, s])\[RightBracketingBar], impurities are too close to host the odd parity solution.

When a=0, two impurities merge to one with doubled strength, so κ=2|(Subscript[a, s])|.

The bounding energies are given by

Plot[{-If[a>0,1+ProductLog[a Exp[-a]]/a,2]^2,-If[a>1,1+ProductLog[-a Exp[-a]]/a,Indeterminate]^2},{a,0,5},GridLines->{None,{-1}},FrameLabel->{"\!\(a/\[LeftBracketingBar]a\_s\[RightBracketingBar]\)","\!\(E/E\_b\)"},PlotRange->{-4,0},ImageSize->200]

The energy E is in unit of the bound state energy -(Subscript[E, b]) of a single Dirac well, with

as given in Eq. ().

As the impurities binds together (a->0), the binding energy (Subscript[E, b])=3|(Subscript[E, 0])| will be released.

Force between impurities are given by

F=-((∂E)/(∂a))

Force can be mediated by exchanging particles.

Bonding orbital: impurity attracts by exchanging particles in even parity state.

Anti-bounding orbital: impurity repels by exchanging particles in odd parity state.

Scattering State (E>0)

T=(((1+((4(1+(k^2)(Subscript[a, s]^2)))/((k^4)(Subscript[a, s]^4)))(sin^2)(k a+φ)))^(-1)),
φ=arctan k (Subscript[a, s]).

Solve the scattering state wavefunction and show Eq. ().

t=((((k (Subscript[a, s])))^2)/((e^(2i k a))-(((1-i k (Subscript[a, s])))^2)))

Block[{ψ1,ψ2,ψ3,sol},
ψ1[x_]:=Exp[I k x]+r Exp[-I k x];
ψ2[x_]:=A Exp[I k x]+B Exp[-I k x];
ψ3[x_]:=t Exp[I k x];
sol=Solve[{ψ1[-a/2]==ψ2[-a/2],ψ2'[-a/2]-ψ1'[-a/2]==2/as ψ2[-a/2],ψ2[a/2]==ψ3[a/2],ψ3'[a/2]-ψ2'[a/2]==2/as ψ2[a/2]},{r,t,A,B}];
t/.sol]

(1/(Sqrt[1+(k^2)(Subscript[a, s]^2)]))=cos φ, ((k (Subscript[a, s]))/(Sqrt[1+(k^2)(Subscript[a, s]^2)]))=sin φ, k (Subscript[a, s])=tan φ,

((1-i k (Subscript[a, s]))/(Sqrt[1+(k^2)(Subscript[a, s]^2)]))=cos φ-i sin φ=(e^(-i φ)).

t=((((k (Subscript[a, s])))^2)/((e^(2i k a))-(1+(k^2)(Subscript[a, s]^2))(e^(-2i φ))))
=(e^(-2i k a))((((k (Subscript[a, s])))^2)/(1-(1+(k^2)(Subscript[a, s]^2))(e^(-2i (k a+φ)))))

T=((|t|)^2)=((((k (Subscript[a, s])))^4)/((1-(1+(k^2)(Subscript[a, s]^2))(e^(2i (k a+φ))))(1-(1+(k^2)(Subscript[a, s]^2))(e^(-2i (k a+φ))))))
=(((k (Subscript[a, s])))^4)(((1-(1+(k^2)(Subscript[a, s]^2))(e^(2i (k a+φ)))-(1+(k^2)(Subscript[a, s]^2))(e^(-2i (k a+φ)))+(((1+(k^2)(Subscript[a, s]^2)))^2)))^(-1))
=(((k (Subscript[a, s])))^4)(((1-(1+(k^2)(Subscript[a, s]^2))(e^(2i (k a+φ)))-(1+(k^2)(Subscript[a, s]^2))(e^(-2i (k a+φ)))+1+2(k^2)(Subscript[a, s]^2)+(k^4)(Subscript[a, s]^4)))^(-1))
=(((k (Subscript[a, s])))^4)((((k^4)(Subscript[a, s]^4)-(1+(k^2)(Subscript[a, s]^2))((e^(2i (k a+φ)))-2+(e^(-2i (k a+φ))))))^(-1))
=(((k (Subscript[a, s])))^4)((((k^4)(Subscript[a, s]^4)-(1+(k^2)(Subscript[a, s]^2))((((e^(i (k a+φ)))-(e^(-i (k a+φ)))))^2)))^(-1))
=(((k (Subscript[a, s])))^4)((((k^4)(Subscript[a, s]^4)-(1+(k^2)(Subscript[a, s]^2))(((2i sin(k a+φ)))^2)))^(-1))
=(((k (Subscript[a, s])))^4)((((k^4)(Subscript[a, s]^4)+4(1+(k^2)(Subscript[a, s]^2))(sin^2)(k a+φ)))^(-1))
=(((1+((4(1+(k^2)(Subscript[a, s]^2)))/((k^4)(Subscript[a, s]^4)))(sin^2)(k a+φ)))^(-1)).

Manipulate[If[as!=0,Plot[Evaluate@Block[{kas,φ},
kas=Sign[as]Sqrt[En];φ=ArcTan[kas];1/(1+4(1+kas^2)/kas^4 {1,Sin[kas/as+φ]^2})],{En,0,16},PlotLabel->StringTemplate["\!\(a\_s/a=`1`\)"][as],FrameLabel->{"\!\(E/E\_b\)","\!\( T\)"},PlotStyle->{Directive[Dashed,AbsoluteThickness[1],Gray],Black},PlotRange->{0,1},PlotRangePadding->Scaled[.02],ImageSize->200],
Plot[Block[{kas,φ},
kas=Sqrt[En];φ=ArcTan[kas];1/(1+4(1+kas^2)/kas^4)],{En,0,16},PlotLabel->StringTemplate["\!\(a\_s/a=`1`\)"][as],FrameLabel->{"\!\(E/E\_b\)","\!\( T\)"},PlotStyle->Black,Filling->1,FillingStyle->Black,PlotRange->{0,1},PlotRangePadding->Scaled[.02],ImageSize->200]],{{as,0.2,Style["\!\(a\_s/a\)","InlineFormula"]},-2,2},SaveDefinitions->True,Paneled->False]

Resonant tunneling (T=1) when (k a+φ) vanishes, which happens at

with n=(0),1,2,…. (n can take 0 only if γ<0 and a<\[LeftBracketingBar](Subscript[a, s])\[RightBracketingBar]), and the corresponding energy is

(Subscript[E, n])=(((ℏ^2)(Subscript[k, n]^2))/(2m)).

Weak potential limit γ->0 (or (Subscript[a, s])→±∞)

Weak well limit γ->(SubMinus[0])

(Subscript[E, n])=(((π^2)(ℏ^2))/(2m (a^2)))(((n+(1/2)))^2), (n=0,1,2,…).

Weak barrier limit γ->(SubPlus[0])

(Subscript[E, n])=(((π^2)(ℏ^2))/(2m (a^2)))(((n-(1/2)))^2), (n=1,2,…).

Strong potential limit γ->∞ (or (Subscript[a, s])->0)

(Subscript[E, n])=(((π^2)(ℏ^2))/(2m (((a+(Subscript[a, s])))^2)))(n^2), (n=1,2,…).

Quasi-bound states formed at these levels, when

n≪|(a/(π (Subscript[a, s])))|.

(Such that the resonance is sharp).

Note that quasi-bound state are still scattering state. The resonant levels are still doubly degenerated (left- and right-moving).

When (Subscript[a, s])=0, one of the quasi-bound states becomes the bound state inside a pair of infinite walls, and the other becomes a resonant state outside the walls.

Periodic Dirac Potential

General Discussion

V(x)=(∑_(l=-∞)^∞)γ δ(x-l a),

γ - the strength of each Dirac potential, which defines the scattering length (Subscript[a, s])

(Subscript[a, s])=((ℏ^2)/(m γ)),

a - lattice spacing (periodicity in space),

l∈ℤ - lattice index.

(K^2)=((2m E)/(ℏ^2)),

ψ(x)=({{{, ({{⋮, ⋮}, {(Subscript[A, l])(e^(i K x))+(Subscript[B, l])(e^(-i K x)), l a<x<(l+1)a,}, {⋮, ⋮}})}})

Bloch’s Theorem

(Subscript[u, k])(x) is periodic: (Subscript[u, k])(x+a)=(Subscript[u, k])(x) (of the same periodicity)

k∈[-π/a,π/a) is the quasi-momentum.

(Subscript[(OverHat[T]), a]^†)(OverHat[p])(Subscript[(OverHat[T]), a])=(OverHat[p]),
(Subscript[(OverHat[T]), a]^†)(OverHat[x])(Subscript[(OverHat[T]), a])=(OverHat[x])-a.

Show Eq. () based on the definition Eq. ().

(OverHat[x]) ((OverHat[p])^n)=((OverHat[p])^n)(OverHat[x])+i ℏ n ((OverHat[p])^(n-1)).

When n=1, (OverHat[x]) (OverHat[p])=(OverHat[p])(OverHat[x])+i ℏ is satisfied by applying the commutation relation directly.

Assuming the relation holds for n, we aim to prove that it holds for n+1,

(OverHat[x]) ((OverHat[p])^(n+1))=(OverHat[x])((OverHat[p])^n)(OverHat[p])
=(((OverHat[p])^n)(OverHat[x])+i ℏ n ((OverHat[p])^(n-1)))(OverHat[p])
=((OverHat[p])^n)(OverHat[x])(OverHat[p])+i ℏ n ((OverHat[p])^n)
=((OverHat[p])^n)((OverHat[p])(OverHat[x])+i ℏ)+i ℏ n ((OverHat[p])^n)
=((OverHat[p])^(n+1))(OverHat[x])+i ℏ (n+1) ((OverHat[p])^n).

Thus by mathematical induction, the relation holds for all natural number n.

f((OverHat[p]))=(∑_(n=0)^∞)(Subscript[f, n])((OverHat[p])^n),

(OverHat[x])f((OverHat[p]))=(∑_(n=0)^∞)(Subscript[f, n])(OverHat[x])((OverHat[p])^n)
=(∑_(n=0)^∞)(Subscript[f, n])(((OverHat[p])^n)(OverHat[x])+i ℏ n ((OverHat[p])^(n-1)))
=(∑_(n=0)^∞)(Subscript[f, n])((OverHat[p])^n)(OverHat[x])+i ℏ (∑_(n=0)^∞)(Subscript[f, n])n ((OverHat[p])^(n-1))
=f((OverHat[p]))(OverHat[x])+i ℏ f'((OverHat[p]))

(Subscript[(OverHat[T]), a]^†)(OverHat[x])(Subscript[(OverHat[T]), a])=(e^(-i (OverHat[p])a / ℏ))(OverHat[x])(e^(i (OverHat[p])a / ℏ))
=(e^(-i (OverHat[p])a / ℏ))((e^(i (OverHat[p])a / ℏ))(OverHat[x])+i ℏ (((i a)/ℏ)(e^(i (OverHat[p])a / ℏ))))
=(e^(-i (OverHat[p])a / ℏ))(e^(i (OverHat[p])a / ℏ))((OverHat[x])-a)
=(OverHat[x])-a.

(Subscript[(OverHat[T]), a]^†)(OverHat[H])(Subscript[(OverHat[T]), a])=(OverHat[H]),

(OverHat[T])(Ket[ψ])=(e^(i k a))(Ket[ψ]).

ψ(x+a)=(e^(i k a))ψ(x).

(Subscript[u, k])(x+a)=(e^(-i k (x+a)))ψ(x+a)=(e^(-i k x))ψ(x)=(Subscript[u, k])(x)

Kronig-Penney Relation

x mod a ensures the function (Subscript[u, k])(x+a)=(Subscript[u, k])(x) to be periodic.

The momentum K±k ensures that in the region between Dirac potentials, the wave momentum should match the itinerant momentum K set by the energy E in Eq. ().

ψ(l a+(SubPlus[0]))=ψ(l a+(SubMinus[0])),
(∂_x) ψ(l a+(SubPlus[0]))-(∂_x) ψ(l a+(SubMinus[0]))=2/(Subscript[a, s])ψ(l a),

(e^(i k l a))(({{1-(e^(i (K-k)a)), 1-(e^(-i (K+k)a))}, {-2(Subscript[a, s]^(-1))+i K (1-(e^(i (K-k)a))), -2(Subscript[a, s]^(-1))-i K (1-(e^(-i (K+k)a)))}}))(({{A}, {B}}))=(({{0}, {0}})).

Derive Eq. () and Eq. ().

Block[{ψ,M},
ψ[l_,x_]:=Exp[I k x](A Exp[I (Κ-k)(x-l a)]+B Exp[-I (Κ+k)(x-l a)]);
M=Simplify@Normal@Last@CoefficientArrays[{ψ[l,l a]-ψ[l-1,l a],Derivative[0,1][ψ][l,l a]-Derivative[0,1][ψ][l-1,l a]-2/as ψ[l,l a]},{A,B}];
Column@(Row/@{{"\!\(M=\)",MatrixForm[M]},
{"\!\(det(E\^\(-I k l a\)M)=\)",FullSimplify[Det[Exp[-I k l a]M]]}})]

K (Subscript[a, s])(cos k a-cos K a)-sin K a=0
=>K (Subscript[a, s])(cos k a-cos K a)=sin K a
=>cos k a-cos K a=(1/(K (Subscript[a, s])))sin K a
=>cos k a=cos K a+(1/(K (Subscript[a, s])))sin K a

Bound states: For E<0, the itinerant momentum becomes imaginary K=i Q, i.e.

(Q^2)=-((2m E)/(ℏ^2)),

the Kronig-Penney relation for Q reads

Energy Bands and Gaps

Manipulate[Block[{f,Κ,root,max,min,edges},
f[Κ_]:=If[as==0,Indeterminate,Piecewise[{{Cos[Κ]+1/(Κ as)Sin[Κ],Κ>0},{Cosh[Κ]+1/(Κ as)Sinh[Κ],Κ<0}},1+1/as]];
root[y_,x__]:=Κ/.FindRoot[f[Κ]==y,{Κ,x}];
max[x__]:=Κ/.Last@FindMaximum[f[Κ],{Κ,x}];
min[x__]:=Κ/.Last@FindMinimum[f[Κ],{Κ,x}];
edges=Which[as<0,{If[as<-0.1,root[1,1/as,-∞,0],1/as(1+Csch[1/as]/Tanh[1/as])],If[as<-0.1,root[-1,1/as,If[as<-3/π^2,π(3+as π^2)/(2+as π^2),0.9/as]],1/as(1-Csch[1/as]/Tanh[1/as])],π,root[1,π,max[π,2π]],2π,root[-1,2π,min[2π,3π]],3π,root[1,3π,max[3π,4π]]},as>0,{root[1,0,π],π,root[-1,min[π,2π],2π],2π,root[1,max[2π,3π],3π],3π,root[-1,min[3π,4π],4π],4π},True,{}];Plot[f[Κ],{Κ,-5,4π},FrameLabel->{""},FrameTicks->{{Automatic,Automatic},{{{-4,"\!\(4I\)"},{-2,"\!\(2I\)"},0,π,2π,3π,4π},{}}},GridLines->{edges,{-1,1}},Epilog->{{Opacity[0.3],Red,Rectangle[{-5,1},{4π,3}],Rectangle[{-5,-1},{4π,-3}]},{Opacity[0.3],Green,Rectangle[##]&@@@Partition[Thread[{edges,Sign[as](-1)^Range@Length@edges}],2]},{Arrowheads[{-0.04,0.04}],AbsoluteThickness[1],{Arrowheads[{0.04}],Arrow[{#,-4.1}&/@{-4,π}],Dashing[{0.01,0.018}],Line[{#,-4.1}&/@{-5.5,-4}],Text["\!\(I Q=K\_1\)",{-0.2,-4.7}]},Table[{Arrow[{#,-4.1}&/@{(n-1)π,n π}],Text[StringTemplate["\!\(K\_`1`\)"][n],{n π-π/2,-4.7}]},{n,2,4}]},Text["\!\(cos k a\)",{-7,0},{0,-1},{0,1}],Text["\!\(K a\)",{4π,-3},{-1.3,0}]},ImagePadding->{{32,25}, {35, 0}},PlotStyle->Black,PlotRange->{{-5,4π},{-3,3}},PlotRangeClipping->False,PlotLabel->StringTemplate["\!\(a\_s/a=`1`\)"][as],ImageSize->220]],{{as,-0.3,Style["\!\(a\_s/a\)","InlineFormula"]},-2,2},SaveDefinitions->True,Paneled->False]

Horizontal axis:

To the right: real K, corresponding to E>0 scattering states.

To the left: imaginary K=i Q, corresponding to E<0 bound states.

Vertical axis:

White zone: feasible region [-1,1]

Red zone: infeasible region

Energy bands: energy values corresponding to the allowed solutions that fall in the green bands, where the right-hand-side function fall in the feasible region [-1,1], such that real solutions of k exist for cos k a.

Energy gaps: energy values that are forbidden, as they result in the right-hand-side function to enter the infeasible region where no real solutions of k exist for cos k a (an imaginary solution of k would indicate a decaying wave that can not propagate into the lattice bulk).

Band Structure

The index n=1,2,3,… labels the energy bands, and is called the band index.

(Subscript[K, n])(k+((2π)/a))=(Subscript[K, n])(k) is a periodic function of k with 2π/a periodicity, because k only enters the Kronig-Penney relation through cos k a, which has that periodicity. Therefore, it is sufficient to restrict k within the first Brillouin zone:

k∈[π/a,-π/a).

Each band index n and quasi-momentum k jointly labels a state, with

Energy:

It is often measured in unit of the recoil energy (Subscript[E, R]):=(((ℏ^2)(π/a))^2)/(2m).

Wavefunction: takes the Bloch wavefunction form

where (Subscript[A, n k]), (Subscript[B, n k]) are null vector solutions of Eq. () (with the solution K=(Subscript[K, n])(k) plugged in).

DynamicModule[{nmax=4},Manipulate[DynamicModule[{f,root,scan,sol},
f[K2_]:=If[as==0,Indeterminate,Cos[Sqrt[K2]]+1/as Sinc[Sqrt[K2]]];
root[y_,x__]:=Block[{K2},Re[K2]/.FindRoot[f[K2]==y,{K2,x},PrecisionGoal->6]];
scan[x0_]:=Block[{dy=0.01,ys},
ys=Range[#,-#,-# dy]&@Re@f[x0];
Interpolation[Thread@{ys,FoldList[root[#2,#1]&,x0,Rest@ys]},InterpolationOrder->1]];
sol=Association@Which[as<0,
Join[{1->If[as<-0.05,scan[root[1,-1/as^2]],Interpolation[{#,-1/as^2}&/@Range[-1,1,0.01],InterpolationOrder->1]]},Table[n->scan[((n-1)π)^2],{n,2,nmax}]],
as>0,
Table[n->scan[(n π)^2],{n,nmax}],
True,{}];
Row@{Plot[Evaluate@Table[sol[n][Cos[k]]/π^2,{n,nmax}],{k,-π,π},PlotRange->{-5,16},PlotLabel->StringTemplate["\!\(a\_s/a=`1`\)"][as],Epilog->{{Red,Opacity[0.5],AbsoluteThickness[1],InfiniteLine[{k0,#}&/@{0,1}]},{Red,AbsolutePointSize[5],Point[{k0,sol[n0][Cos[k0]]/π^2}]}},FrameLabel->{"\!\(k a\)","\!\(E\_n(k)/E\_R\)"},Axes->False,GridLines->{{-π,0,π},Range[0,4]^2},PlotStyle->Normal[SparseArray[n0->Red,{nmax},Black]],AspectRatio->1.8,ImageSize->{Automatic,220}],
DynamicModule[{Κ,A,B},
Κ=Sqrt[sol[n0][Cos[k0]]];
{A,B}=Last@Eigenvectors[{{1-Exp[I(Κ-k0)],1-Exp[-I(Κ+k0)]},{-2/as+I Κ(1-Exp[I(Κ-k0)]),-2/as-I Κ(1-Exp[-I(Κ+k0)])}}];
Manipulate[Plot[Re[Exp[I k0 x-I 2 π t](A Exp[I(Κ-k0)Mod[x,1]]+B Exp[-I(Κ+k0)Mod[x,1]])],{x,-5,5},GridLines->{Range[-5,5],{}},PlotRange->{-2,2},PlotStyle->Red,FrameLabel->{"\!\( x\)","\!\(ψ\_\(n k\)(x)\)"},ImageSize->180],{{t,1/8,Style["   \!\(t/T\)","InlineFormula"]},0,1,ImageSize->Small},Paneled->False]]}
],{{as,-0.3,Style["\!\(a\_s/a\)","InlineFormula"]},-2,2},{{n0,1,Style["\!\( n\)","InlineFormula"]},Range[nmax]},{{k0,0,Style["\!\( k\)","InlineFormula"]},-π,π,ImageSize->Small},SaveDefinitions->True,Paneled->False]]

The quasi-momentum k characterizes the momentum of the envelope wave of the wavefunction, omitting the intricate and sharp variations in the wavefunction due to the potential scattering. It reflects the large-scale coherent propagation of a quantum particle in the lattice, capturing its effective momentum.

The energy gap opens at k=0 or k=π/a, where Bloch waves become standing waves, which have nodes (points of zero amplitude) and anti-nodes (points of max amplitude).

If the standing wave places its anti-nodes on the potential sites, the particle experiences the potential interaction with maximal likelihood. Depending on whether the potential is attractive or repulsive, the energy will be lowered or raised.

Conversely, if the standing wave places its nodes on the potential sites, the particle effectively avoids the potential — as if the periodic potential is invisible — and its energy remains unaffected.

These two standing wave configurations share the same quasi-momentum but experience different energy shifts due to their different spatial alignment with the potential. This energy difference defines the energy gap, which quantifies the strength of potential scattering.

General Piecewise Potentials

Transfer Matrix Method

Generic Potential Profile

Within each region, the potential is constant

(Subscript[V, Θ])(x)=(Piecewise[{{(Subscript[V, 1]), x∈(-∞,(Subscript[x, 1]))}, {(Subscript[V, i]), x∈((Subscript[x, i-1]),(Subscript[x, i]))}, {(Subscript[V, n]), x∈((Subscript[x, n-1]),+∞)}}]).

At the region interface, we may have δ-function potentials

(Subscript[V, δ])(x)=(∑_(i=1)^(n-1))(Subscript[γ, i])δ(x-(Subscript[x, i])).

V(x)=(Subscript[V, Θ])(x)+(Subscript[V, δ])(x).

Interface Conditions

(-((ℏ^2)/(2m))(∂_x^2)+V(x)) ψ(x)=E ψ(x).

Prove Eq. ().

V(x∼a)=(Subscript[V, 0])+W Θ(x-a)+γ δ(x-a).

(∂_x) ψ(a+(SubPlus[0]))-(∂_x) ψ(a-(SubPlus[0]))
=(∫_(a-(SubPlus[0]))^(a+(SubPlus[0])))2((Subscript[V, 0])-E+W Θ(x-a)+γ δ(x-a))ψ(x)dx
=2γ ψ(a).

(∂_x) ψ(x∼a)=(Subscript[ψ, 0]^′)+2γ ψ(a)Θ(x-a),

ψ(a+(SubPlus[0]))- ψ(a-(SubPlus[0]))
=(∫_(a-(SubPlus[0]))^(a+(SubPlus[0])))((Subscript[ψ, 0]^′)+2γ ψ(a)Θ(x-a))dx
=0.

Matching Wave Functions

(({{1, 1}, {i (Subscript[k, i+1]), -i (Subscript[k, i+1])}}))(({{(e^(i (Subscript[k, i+1])(Subscript[x, i]))), 0}, {0, (e^(-i (Subscript[k, i+1])(Subscript[x, i])))}}))(({{(Subscript[A, i+1])}, {(Subscript[B, i+1])}}))
=(({{1, 1}, {i (Subscript[k, i])+2(Subscript[γ, i]), -i (Subscript[k, i])+2(Subscript[γ, i])}}))(({{(e^(i (Subscript[k, i])(Subscript[x, i]))), 0}, {0, (e^(-i (Subscript[k, i])(Subscript[x, i])))}}))(({{(Subscript[A, i])}, {(Subscript[B, i])}})).

(({{(Subscript[A, i+1])}, {(Subscript[B, i+1])}}))=(Subscript[S, i+1←i])(({{(Subscript[A, i])}, {(Subscript[B, i])}})),

(Subscript[S, i+1←i])=(1/(2(Subscript[k, i+1])))(({{((Subscript[k, i+1])+(Subscript[k, i])-2i (Subscript[γ, i]))(e^(i ((Subscript[k, i])-(Subscript[k, i+1]))(Subscript[x, i]))), ((Subscript[k, i+1])-(Subscript[k, i])-2i (Subscript[γ, i]))(e^(-i ((Subscript[k, i])+(Subscript[k, i+1]))(Subscript[x, i])))}, {((Subscript[k, i+1])-(Subscript[k, i])+2i (Subscript[γ, i]))(e^(i ((Subscript[k, i])+(Subscript[k, i+1]))(Subscript[x, i]))), ((Subscript[k, i+1])+(Subscript[k, i])+2i (Subscript[γ, i]))(e^(-i ((Subscript[k, i])-(Subscript[k, i+1]))(Subscript[x, i])))}})).

Derive Eq. ().

(Subscript[S, i+1←i])=(((({{(e^(i (Subscript[k, i+1])(Subscript[x, i]))), 0}, {0, (e^(-i (Subscript[k, i+1])(Subscript[x, i])))}})))^(-1))(((({{1, 1}, {i (Subscript[k, i+1]), -i (Subscript[k, i+1])}})))^(-1))(({{1, 1}, {i (Subscript[k, i])+2(Subscript[γ, i]), -i (Subscript[k, i])+2(Subscript[γ, i])}}))(({{(e^(i (Subscript[k, i])(Subscript[x, i]))), 0}, {0, (e^(-i (Subscript[k, i])(Subscript[x, i])))}})).

Inverse[DiagonalMatrix[{Exp[I k2 x],Exp[-I k2 x]}]].Inverse[{{1,1},{I k2,-I k2}}].{{1,1},{I k1+2γ,-I k1+2 γ}}.DiagonalMatrix[{Exp[I k1 x],Exp[-I k1 x]}]//Simplify

(({{(Subscript[A, i-1])}, {(Subscript[B, i-1])}}))=(Subscript[S, i-1←i])(({{(Subscript[A, i])}, {(Subscript[B, i])}})),

(Subscript[S, i-1←i])=(1/(2(Subscript[k, i-1])))(({{((Subscript[k, i-1])+(Subscript[k, i])+2i (Subscript[γ, i-1]))(e^(i ((Subscript[k, i])-(Subscript[k, i-1]))(Subscript[x, i]))), ((Subscript[k, i-1])-(Subscript[k, i])+2i (Subscript[γ, i-1]))(e^(-i ((Subscript[k, i])+(Subscript[k, i-1]))(Subscript[x, i])))}, {((Subscript[k, i-1])-(Subscript[k, i])-2i (Subscript[γ, i-1]))(e^(i ((Subscript[k, i])+(Subscript[k, i-1]))(Subscript[x, i]))), ((Subscript[k, i-1])+(Subscript[k, i])-2i (Subscript[γ, i-1]))(e^(-i ((Subscript[k, i])-(Subscript[k, i-1]))(Subscript[x, i])))}}))

(Subscript[S, i←c])=({{{, ({{(Subscript[S, i←i+1])(Subscript[S, i+1←i+2])… (Subscript[S, c-1←c]), if c>i,}, {(Subscript[S, i←i-1])(Subscript[S, i-1←i-2])… (Subscript[S, c+1←c]), if c<i.}})}})

Boundary Regions

If (Subscript[V, 1])>E, momentum (Subscript[k, 1])=i (Subscript[κ, 1])∈𝕀 is imaginary => ψ(x<(Subscript[x, 1]))=(Subscript[A, 1])(e^(-(Subscript[κ, 1])x))+(Subscript[B, 1])(e^((Subscript[κ, 1])x))=> (Subscript[A, 1])=0.

(({{1, 0}}))(({{(Subscript[A, 1])}, {(Subscript[B, 1])}}))=(Subscript[v, l])(({{(Subscript[A, l])}, {(Subscript[B, l])}}))=0,

with (Subscript[v, l])=(({{1, 0}})).

If (Subscript[V, 1])=∞ => (Subscript[A, 1])=(Subscript[B, 1])=0; left-most region shifts to l=2

(({{(e^(i (Subscript[k, 2])(Subscript[x, 1]))), (e^(-i (Subscript[k, 2])(Subscript[x, 1])))}}))(({{(Subscript[A, 2])}, {(Subscript[B, 2])}}))=(Subscript[v, l])(({{(Subscript[A, l])}, {(Subscript[B, l])}}))=0.

with (Subscript[v, l])=(({{(e^(i (Subscript[k, 2])(Subscript[x, 1]))), (e^(-i (Subscript[k, 2])(Subscript[x, 1])))}})).

If (Subscript[V, n])>E, momentum (Subscript[k, n])=i (Subscript[κ, n])∈𝕀 is imaginary => ψ(x>(Subscript[x, n-1]))=(Subscript[A, n])(e^(-(Subscript[κ, n])x))+(Subscript[B, n])(e^((Subscript[κ, n])x))=> (Subscript[B, n])=0.

(({{0, 1}}))(({{(Subscript[A, n])}, {(Subscript[B, n])}}))=(Subscript[v, r])(({{(Subscript[A, r])}, {(Subscript[B, r])}}))=0,

with (Subscript[v, r])=(({{0, 1}})).

If (Subscript[V, n])=∞ => (Subscript[A, n])=(Subscript[B, n])=0; right-most region shifts to r=n-1

(({{(e^(i (Subscript[k, n-1])(Subscript[x, n-1]))), (e^(-i (Subscript[k, n-1])(Subscript[x, n-1])))}}))(({{(Subscript[A, n-1])}, {(Subscript[B, n-1])}}))=(Subscript[v, r])(({{(Subscript[A, r])}, {(Subscript[B, r])}}))=0.

with (Subscript[v, r])=(({{(e^(i (Subscript[k, n-1])(Subscript[x, n-1]))), (e^(-i (Subscript[k, n-1])(Subscript[x, n-1])))}})).

The boundary conditions can be in general expressed as

(Subscript[v, l])(({{(Subscript[A, l])}, {(Subscript[B, l])}}))=(Subscript[v, l])(Subscript[S, l←c])(({{(Subscript[A, c])}, {(Subscript[B, c])}}))=0,
(Subscript[v, r])(({{(Subscript[A, r])}, {(Subscript[B, r])}}))=(Subscript[v, r])(Subscript[S, r←c])(({{(Subscript[A, c])}, {(Subscript[B, c])}}))=0,

which combines into a single condition P·((Subscript[A, c]),(Subscript[B, c]))ᵀ=0, where P is a 2×2 matrix given by

Solutions for ((Subscript[A, c]),(Subscript[B, c])) span the null space of P.

Once the solutions of ((Subscript[A, c]),(Subscript[B, c])) is found in the reference region, ((Subscript[A, i]),(Subscript[B, i])) in any other regions can be obtained by applying the transfer matrix, following Eq. (). Then the wave function can be reconstructed from Eq. ().

If P is full rank, we will need to tune the energy E to make

such that E will correspond to an eigen energy level, and the null vectors of P(E) will form the solutions of ((Subscript[A, c]),(Subscript[B, c])), which can be used to construct the eigen wave function.

Applications

Harmonic Oscillator

V(x)=(1/2)m (ω^2)(x^2).

Approximate V(x) by a piecewise potential:

V(x)=({{{, ({{+∞, x∈(-∞,(Subscript[x, 1]))}, {(1/2)((((((Subscript[x, i-1])+(Subscript[x, i]))/2)))^2), x∈((Subscript[x, i-1]),(Subscript[x, i]))}, {+∞, x∈((Subscript[x, n-1]),+∞)}})}})

Block[{Vmax=7,dV=0.5,xs,Vdat},
xs=Sqrt[2Range[0,Vmax,dV]];
Vdat=Quiet@PotentialData[Piecewise[Flatten[{{((#1+#2)/2)^2/2,#1<=x<#2},{((#1+#2)/2)^2/2,-#2<x<=-#1}}&@@@Thread@{Most[xs],Rest[xs]},1],∞],x];
Graphics[{PotentialPlot[Vdat,{-4,4},{-0.5,Vmax}],{Red,Opacity[0.5],Line[{#,#^2/2}&/@Range[-4,4,0.1]]}},PlotRange->{{-4,4},{-0.5,8}},PlotRangeClipping->True,PlotRangePadding->0,Frame->True,FrameLabel->{"\!\( x\)","\!\(V(x)\)"},AspectRatio->0.7,ImageSize->180]]

Scan |det P(E)| (in logarithmic scale) to find its zeros, see Eq. (), which corresponds to the eigen energies in the potential well.

Block[{Vmax=7,dV=0.5,xs,Vdat},
xs=Sqrt[2Range[0,Vmax,dV]];
Vdat=Quiet@PotentialData[Piecewise[Flatten[{{((#1+#2)/2)^2/2,#1<=x<#2},{((#1+#2)/2)^2/2,-#2<x<=-#1}}&@@@Thread@{Most[xs],Rest[xs]},1],∞],x];
ListLinePlot[Table[{Ε,Log@Abs@Quiet@Det[TransferMatrix[Vdat,Ε]["P"]]},{Ε,0,5,0.02}],
Axes->False,PlotStyle->Black,GridLines->{Range[0,4]+1/2,None},FrameLabel->{"\!\( E/ℏ ω\)","\!\(log\[LeftBracketingBar]det P(E)\[RightBracketingBar]\)"},ImageSize->180]]

Energy levels appear approximately at

(Subscript[E, n])≈(n+(1/2))ℏ ω,

matching the expectation (with ℏ ω=1).

Wave functions can be reconstructed at eigen energy levels, demonstrating how any 1D wave function can be approximated by piecing exponential functions (e^(±i k x)) or (e^(±κ x)) together.

WaveExplorer[Block[{Vmax=7,dV=0.5,xs},
xs=Sqrt[2Range[0,Vmax,dV]];
Piecewise[Flatten[{{((#1+#2)/2)^2/2,#1<=x<#2},{((#1+#2)/2)^2/2,-#2<x<=-#1}}&@@@Thread@{Most[xs],Rest[xs]},1],∞]],{x,-4,4},0.5,{-0.5,8}]

Crystal and Band Structure

V(x)=({{{, ({{(Subscript[V, 0]), a i-b/2<x<a i+b/2 for i∈ℤ}, {0, otherwise}})}})

a - lattice spacing: distance between adjacent potentials,

b  - potential width: width of each single potential,

(Subscript[V, 0]) - potential height: energy inside the potential region.

Block[{V0=10,a=1,b=0.4,L=5},
WaveExplorer[Piecewise[Table[{V0,i a-b/2<x<i a+b/2},{i,-L,L}],0],{x,-a (L+3),a(L+3)},3.165,{-0.1V0,1.5V0}]]

Incident wave energy lies in an energy band => wave can propagate through the crystal => the crystal appears transparent.

Incident wave energy falls in an energy gap => wave propagation is forbidden in the crystal => the crystal becomes opaque.

Block[{V0=10,a=1,b=0.4,L=5,Vdat,T},
Vdat=Quiet@PotentialData[Piecewise[Table[{V0,i a-b/2<x<i a+b/2},{i,-L,L}],0],x];
T[Ε_?NumericQ]:=Module[{TM0},
TM0=TransferMatrix[Vdat,Ε];
If[PossibleZeroQ@Abs@Det[TM0@"P"],
1/Abs[First@Fold[#2.#1&,{1,0},TM0@"Srs"]]^2,0]];
Plot[T[Ε],{Ε,0,15},PlotRange->{0,1},PlotRangePadding->0.05,PlotStyle->Black,MaxRecursion->7,FrameLabel->{"\!\( E\)","\!\(T(E)\)"},ImageSize->200]]

Disorder and Localization

Block[{V0=10,a=1,b=0.4,L=5},
SeedRandom[6];
WaveExplorer[Piecewise[Table[{RandomReal[{-V0,V0}]+V0,i a-b/2<x<i a+b/2},{i,-L,L}],0],{x,-a (L+3),a(L+3)},4.6,{-0.1V0,2V0}]]

Transmission coefficient is largely suppressed in the presence of disorder.

Block[{V0=10,a=1,b=0.4,L=5,Vdat,T},
SeedRandom[6];
Vdat=Quiet@PotentialData[Piecewise[Table[{RandomReal[{-V0,V0}]+V0,i a-b/2<x<i a+b/2},{i,-L,L}],0],x];
T[Ε_?NumericQ]:=Module[{TM0},
TM0=TransferMatrix[Vdat,Ε];
If[PossibleZeroQ@Abs@Det[TM0@"P"],
1/Abs[First@Fold[#2.#1&,{1,0},TM0@"Srs"]]^2,0]];
Plot[T[Ε],{Ε,0,15},PlotRange->{0,1},PlotRangePadding->0.05,PlotStyle->Black,MaxRecursion->7,FrameLabel->{"\!\( E\)","\!\(T(E)\)"},ImageSize->200]]

Localization: A wave become localized if it decays exponentially away from some local region, rather than spreading throughout the space.

|ψ(x)|∼(e^(-|x-(Subscript[x, 0])|/ξ)),

where ξ is the localization length.

Disorder and Anderson localization: in 1D, even weak disorder can dramatically affect wave propagation,

Every potential scatters part of the wave,

Randomness causes phases to misalign, suppressing constructive interference,

Over long distance, scattered waves interfere destructively with the forward propagating wave, effectively trapping the wave.

Code Section

LeftWall[x_]:=Piecewise[{{∞,x<0}}];
RightWall[x_]:=Piecewise[{{∞,x>0}}];
SyntaxInformation[PotentialData]={
"ArgumentsPattern"->{_,_},
"LocalVariables"->{"Solve",{2}}};
PotentialData[expr_,x_]:=Module[{par,xs,Vs,γs,leg},
par[(UnitStep|HeavisideTheta|DiracDelta)[v_,___]]:=
With[{sol=Solve[v==0,x,Reals]},
If[sol=!={},Sow[x/.sol]];];
par[HeavisidePi[v_,___]]:=
With[{sol=Solve[v==1/2||v==-1/2,x,Reals]},
If[sol=!={},Sow[x/.sol]];];
par[HoldPattern[Piecewise[lst_,___]]]:=Module[{conds,cpar},
conds=Reduce[#,x]&/@Last/@lst;
cpar[(Greater|Less|GreaterEqual|LessEqual)
[a___,x,b___]]:=Sow[{a,b}];
cpar[HoldPattern[Inequality[a_,_,x,_,b_]]]:=Sow[{a,b}];
cpar[p_]:=cpar/@p;
cpar/@conds];
par[p_]:=par/@p;
xs=Sort@DeleteDuplicates@Flatten@Last@Reap[par@expr];
Vs=expr/.x->#&/@Join[{First@xs-1},
Mean/@Partition[xs,2,1],{Last@xs+1}];
γs=Coefficient[Chop[expr/.x->#],DiracDelta[0]]&/@xs;
leg[Vdat_]:=Which[
MatchQ[Vdat["Vs"],{∞,∞,__}],
Rest/@Vdat,
MatchQ[Vdat["Vs"],{__,∞,∞}],
Most/@Vdat,
MatchQ[Vdat["Vs"],{∞..}],
<|"xs"->{},"Vs"->{0},"γs"->{}|>,
True,Vdat];
FixedPoint[leg,<|"xs"->xs,"Vs"->Vs,"γs"->γs|>,Length@xs]];

arr[x_,{y1_,y2_},z_]:=Which[
z>0,Translate[{Arrow[Offset/@{{0,0},{0,20}}],
Text[Style["+"<>ToString[Abs[z]],10],Offset[{0,20}],{0,-1}]},{x,Max[y1,y2]}],
z<0,Translate[{Arrow[Offset/@{{0,0},{0,-20}}],
Text[Style["-"<>ToString[Abs[z]],10],Offset[{0,-20}],{0,1}]},{x,Min[y1,y2]}],
True,Nothing];
PotentialPlot[Vdat_,{xmin_,xmax_},{Εmin_,Εmax_}]:=
Module[{Vs,xs,pts},
Vs=Clip[Vdat@"Vs",{Εmin,Εmax+Abs[Εmax-Εmin]}];
xs=Join[{Min[xmin,Min@Vdat@"xs"]},Vdat@"xs",{Max[xmax,Max@Vdat@"xs"]}];
pts=First@Last@Reap[Nest[(Sow[First/@#];If[Length[#[[1]]]>Length[#[[2]]],{Rest[#[[1]]],#[[2]]},{#[[1]],Rest[#[[2]]]}])&,{xs,Vs},Length[xs]+Length[Vs]-1]];
{{LightGray,Polygon[Join[{{pts[[1,1]],Εmin}},pts,{{pts[[-1,1]],Εmin}}]]},{Black,Line[pts]},
{Arrowheads[Medium],arr@@@Thread[{Vdat@"xs",
Partition[Vdat@"Vs",2,1],Vdat@"γs"}]}}];

TM=Compile[{{x,_Real},{γ,_Real},{k1,_Complex},{k2,_Complex}},
1/(2k2)DiagonalMatrix[Exp[I k2 x{-1,1}]].
{{k2+k1-2I γ,k2-k1-2I γ},{k2-k1+2I γ,k2+k1+2I γ}}.
DiagonalMatrix[Exp[I k1 x{1,-1}]]];
invTM=Compile[{{x,_Real},{γ,_Real},{k1,_Complex},{k2,_Complex}},
1/(2k1)DiagonalMatrix[Exp[I k1 x{-1,1}]].
{{k1+k2+2I γ,k1-k2+2I γ},{k1-k2-2I γ,k1+k2-2I γ}}.
DiagonalMatrix[Exp[I k2 x{1,-1}]]];
TransferMatrix[Vdat_,Ε_?NumericQ,x0_:Automatic]:=Module[{n,V1,Vn,l,r,c,vl,vr,ks,xdat,Sls,Srs,P},
n=Length@Vdat["Vs"];
{V1,Vn}=Vdat["Vs"][[{1,n}]];
{l,r}={1,0};
{vl,vr}={{UnitStep[Re[V1-Ε]],0},{0,UnitStep[Re[Vn-Ε]]}};
ks=Sqrt[2(Ε-Vdat@"Vs")];
If[V1===∞,l=2;vl=Exp[I ks[[2]]Vdat["xs"][[1]]{1,-1}]];
If[Vn===∞,r=-1;vr=I Exp[I ks[[-2]]Vdat["xs"][[-1]]{1,-1}]];
If[x0===Automatic,c=l-1,c=Length[Cases[Vdat@"xs",x_/;x<x0]]];
xdat=Thread[{Vdat@"xs",Vdat@"γs",Most@ks,Rest@ks}];
Sls=invTM@@@Reverse@Take[xdat,{l,c}];
Srs=TM@@@Take[xdat,{c+1,r-1}];
P={vl.Fold[#2.#1&,IdentityMatrix[2],Sls],vr.Fold[#2.#1&,IdentityMatrix[2],Srs]};
<|"n"->n,"l"->l,"r"->r,"c"->c,"ks"->ks,"P"->P,"Sls"->Sls,"Srs"->Srs|>];

Off[Det::luc];
EigenPair[Vdat_,Ε0_,x0_:Automatic]:=Module[{δ=1*^-10,ϵ=1*^-3,extended,TM0,Ε1,TM1,nsp,mkWF},
TM0=TransferMatrix[Vdat,Ε0+δ I,x0];
extended=Abs@Det[TM0@"P"]<δ;
If[extended,
{Ε1,TM1}={Ε0,TM0},
Module[{Εmin,Εmax,logdet,Ε},
Εmin=Min[Min@Vdat@"Vs",Min[# Abs[#]&/@Vdat@"γs"]];
Εmax=Min@Vdat["Vs"][[{1,-1}]];
logdet[Ε_?NumericQ]:=
Log@Abs@Det[TransferMatrix[Vdat,Ε+δ I,x0]@"P"];
Ε1=Ε/.Last@Quiet@FindMinimum[logdet[Ε],{Ε,Clip[Ε0,{Εmin,Εmax}],Εmin,Εmax},Method->"PrincipalAxis"];
TM1=TransferMatrix[Vdat,Ε1+δ I,x0]]];
Module[{maxP,Pvals,Pvecs},
maxP=Max@Flatten@Abs[TM1@"P"];
{Pvals,Pvecs}=Eigensystem[TM1@"P"];
nsp=Pick[Pvecs,(Abs[#]<ϵ Max[1,maxP])&/@Pvals,True]];
mkWF[TM_,nv_]:=Module[{Als,Ars,As,Ws,Rs,WF},
Als=FoldList[#2.#1&,nv,TM@"Sls"];
Ars=FoldList[#2.#1&,nv,TM@"Srs"];
As=Join[Reverse@Rest[Als],Ars];
Block[{x},
Ws=Chop[MapThread[Dot,{As,Table[Exp[I k x{1,-1}],
{k,N[TM@"ks"][[TM@"l";;TM@"r"-1]]}]}],100δ];
Rs=Join[{x≤First@Vdat@"xs"},
#1≤x≤#2&@@@Partition[Vdat@"xs",2,1],
{x≥Last@Vdat@"xs"}][[TM@"l";;TM@"r"-1]];
WF=Piecewise[Thread[{Ws,Rs}],0];
Function[Evaluate[WF/.x->Slot[1]]]]];
<|"Ε"->Ε1,"ψs"->(mkWF[TM1,#]&/@nsp),"extended"->extended|>];

WFNormalize[ψ_,{xmin_,xmax_}]:=Module[{Z1,Z2},
Block[{x},
Z1=Max@Table[Abs[ψ[x]],{x,Range[xmin,xmax,(xmax-xmin)/32]}];
Z2=FindMaxValue[Abs[ψ[x]],{x,xmin,xmax}]];
ψ[#]/Max[Z1,Z2]&];

Options[WaveExplorer]=Join[Options[Graphics],
{"Scale"->0.1}];
SyntaxInformation[WaveExplorer]={
"ArgumentsPattern"->{_,{_,_,_},_,{_,_},OptionsPattern[]},
"LocalVariables"->{"Manipulate",{2,∞},"Lexical"}};
WaveExplorer[expr_,{x_,xmin_,xmax_},Ε0_,{Εmin_,Εmax_},
opt:OptionsPattern[WaveExplorer]]:=WaveExplorer[expr,{{x,Automatic},xmin,xmax},Ε0,{Εmin,Εmax},
opt];
WaveExplorer[expr_,{{x_,x0_:Automatic},xmin_,xmax_},Ε0_,{Εmin_,Εmax_},
opt:OptionsPattern[WaveExplorer]]:=
DynamicModule[{Vdat,Vplt,Ε=Ε0,t=0,EP},
Vdat=Quiet@PotentialData[expr,x];
Vplt=PotentialPlot[Vdat,{xmin,xmax},{Εmin,Εmax}];
EP=EigenPair[Vdat,Ε,x0];Ε=EP@"Ε";
Column[Row/@{
{Style["\!\(E:\) ","InlineFormula"],Manipulator[Dynamic[Ε,{Automatic,
Module[{},EP=EigenPair[Vdat,#1,x0];Ε=EP@"Ε";]&}],
{Max[Εmin,Min[Vdat@"Vs"]+1.*^-10+
Total[-#^2UnitStep[-#]/2&/@Vdat@"γs"]],
Εmax},ImageSize->180]},
{Style["\!\(t:\) ","InlineFormula"],Animator[Dynamic[t],{0,1},
AnimationRate->0.5,AnimationRunning->False,
ImageSize->{200,16}]},
{Dynamic[
Graphics[
{Vplt,
{AbsoluteThickness[1],If[EP@"extended",Green,Red],
Line[{#,EP@"Ε"}&/@{xmin,xmax}]},
First@Plot[Evaluate@Table[
Re[OptionValue["Scale"](Εmax-Εmin)Exp[-I 2π t]
WFNormalize[ψ,{xmin,xmax}][x]]+EP@"Ε",
{ψ,EP@"ψs"}],{x,xmin,xmax},PlotRange->All,PlotStyle->{Yellow,Blue}]},
PlotRange->{{xmin,xmax},{Εmin,Εmax}},
FilterRules[Flatten@{opt},Options[Graphics]],
PlotRangeClipping->True,
PlotLabel->Style[StringTemplate["\!\(E=`1`\)"][NumberForm[Ε,{3,2}]],Black,14],
Frame->True,FrameLabel->{"\!\( x\)","\!\(V(x)\)"},
AspectRatio->0.7,ImageSize->200]
]}}],
SaveDefinitions->True];