The plots below demonstrate the band edge wavefunctoins for the dirac barrier lattice as a function of the strength of the barriers.
	The slider at the top left controls the strength of the barriers.

The first plot presents the band edge energies as horizontal lines, filling in the bands corresponding to allowed states with a rainbow of colors.
	Note: MMatica’s algorithm for filing barfs when the bands are narrow -- improving that slows things down, so let’s not bother.

The next three plots present |ψ((|)^2) for the bottom of the band in the first three bands (blue lines) over two wells in the lattice.  
	Observe that the weak-barrier state has a maximum at the deltas (x=-1,0,1,...), but develops kinks at the deltas for finite strength, which kinks become nodes at strong coupling.
	The barriers are indicated by the black lines; the opacity of the lines reflects the strength, (Subscript[g, o]), of the barriers.

The last three plots show a single barrier+well’s worth of the probability density |ψ((|)^2)
	Note: Included are |ψ((|)^2) for the free wavefunctions corresponding to the top (yellow) and bottom(red)  of the bands.  
	Note: the bottom of the band (blue) interpolates between them as you sweep go.  Why?

For stupid technical reasons, you must move the slider very, very slowly or the code will panic and plot crazy things.  
	If the plots freak out, run the command again and try more slowly.  
	DO NOT REMOVE THE ko=0.01; LINE!  It is necessary initialization.

Finally, if you’d like to see ten wells worth of wavefunctions, comment out the PlotBandEdges and uncomment the PlotBandEdgeWFs command, then run again.
	DO NOT REMOVE THE ko=0.01; LINE!  It is again necessary initialization.

Code: Run the entire notebook to initialize

Analytics: Solve the eigenvalue equation and determine the form of the wavefunction and band equation.

Input the equations we derived in lecture determining the boundary conditions at the delta function at the origin:

EQs={A+B ==  A (E^(I (k-q) L))+ B (E^(-I (k+q) L)),
I k(A-B )- I k(A (E^(I (k-q) L))- B (E^(-I (k+q) L)))==(go/L)(A+B)}/.{L->1,A->1};

Note that we set L=1.  This is equivalent to measuring all lengths in units of L and momenta in units of (1/L).  We also set A=1 to simplify algebra and normalization.
We’re also setting ℏ and m to 1, corresponding to measuring all energies and momenta in terms of dimensionless variables divided by suitable powers of ℏ and m.

Use the first equation to solve for B.  This lets us write the wavefunction in closed form.

BRule = Solve[EQs[[1]],{B}][[1]]//Simplify

Plug this into the second equation, then solve for Cos[q] to get the band energy condition:

EnergyEQ=Cos[q]==(Cos[q]/.Solve[  EQs[[2]]/.BRule //FullSimplify  ,Cos[q]][[1]])//Simplify

Now let’s build the wavefunction at energy E=(ℏ^2)(k^2) and crystal momentum ℏ q.

ψbot[k_,q_][x_]:=(E^(I q Floor[x]))(1/2)((E^(I k Mod[x,1]))+(((E^(2 I k))-(E^(I (k+q))))/(-1+(E^(I (k+q))))) (E^(-I k Mod[x,1])));

Here, Floor[x] counts the integer part of x, and thus determines how many lattice sites we are away from the origin, while Mod[x,1] is the position within the given well.
The overall phase (e^(i q Floor[x])) thus takes care of the Bloch periodicity condition.

We’ll also need the even and odd free-particle wavefunctions corresponding to the resonant wavelengths (Subscript[k, n]) = (n/π).  Since the top of the band corresponds to q = 0, π, 0, π, 0,... in alternating sequence, we need to keep track of the phase.

ψtop[n_][x_]:=(1/2)((E^(I n π x))+(((-1))^n) (E^(-I n π x)));

Similarly, the state that would have been at the bottom of the (n^th) band is the even/odd compliment of the top of the (((n-1))^th) band:

ψtopm[n_][x_]:=(1/2)((E^(I n π x))+(((-1))^(n+1)) (E^(-I n π x)));

Numerics: Plotting the energies and eigenfunctions

First, let’s set some parameters that we’ll need later:

AR = 8; (* aspect ratio of many-well wavefunction plots *)
dpi=52; (* good for a projector at 1024x768 *)
dpi=72; (* good for a big display *)
gmin=0.01;gmax=5.5;

PlotBandEdges : A macro to plot the band edge energies and eigenfunctions as a function of barrier strength

TwoWells:=Line[{{-1,1.2},{-1,0},{0,0},{0,1.2},{0,0},{1,0},{1,1.2}}];
TwoWellOptions:={PlotRange->{{-1,1},{-1.2,1.2}},Axes->False,Frame->True,ImageSize->{dpi*3*GoldenRatio,dpi*3}};

PlotBandEdges:=Manipulate[
(*First, find the energies (Subscript[E, go])~(k^2) for a given barrier height. 
 Several things to note:
a) the manipulated variable is go.  
this means the slider controls go linearly.
this is annoying, since we want to spend more time loooking at the small-go behavior
but still see the results for very large go.
we thus replace (go/2) with ((go^2)+(go^3)) instead of  to  the manipulated variable go.
b) note that Cos[qL]=(((-1))^n) at the bottom of the (n^th) band.
c) finally, note that find-root requires an initial guess.  this makes sense: there are many roots; the 
initial guess should help find-root find the desired root.  we thus need good guesses.  our
first guess comes from an input value of ko, plus nπ for the nth band minimum.
thereafter, the previous value of ko, again plus nπ, is used.  
working through why these are good guesses is a good check on your understanding.
*)
ko=k/.FindRoot[Cos[k]+(((go^2)+(go^3))/k)Sin[k]==1,{k,ko}]; (* qL = 0 *)
k1=k/.FindRoot[Cos[k]+(((go^2)+(go^3))/k)Sin[k]==-1,{k,π+ko}];(* qL = π *)
k2=k/.FindRoot[Cos[k]+(((go^2)+(go^3))/k)Sin[k]==1,{k,2*π+ko}];(* qL = 2π ~ 0 *)
k3=k/.FindRoot[Cos[k]+(((go^2)+(go^3))/k)Sin[k]==-1,{k,3*π+ko}];(* qL = 3π ~ π *)
(* Now that we have the roots, let's plot the energies and wavefunctions in a big list: *)
{(* 
First, let's plot the energies of the band bottoms, the band tops, and shade in the bands between them. 
*)
Show[
Plot[{(ko^2),(k1^2),(k2^2),(k3^2),0,(π^2),(((2*π))^2),(((3*π))^2),(((4*π))^2)},{x,0,1}],
RegionPlot[(ko^2)<y<(π^2)||(k1^2)≤y≤(((2*π))^2)||(k2^2)≤y≤(((3*π))^2)||(k3^2)≤y≤(((4*π))^2),
{x,0,1},{y,-0.1,(((4*π))^2)+0.1},Mesh->None,BoundaryStyle->None,ColorFunction->"Rainbow",PlotStyle->Opacity[0.2],PerformanceGoal->"Quality"],
Axes->False,Frame->True,ImageSize->{dpi*3*GoldenRatio,dpi*3}
],
(* 
Now let's plot the wavefunctions of the states at the bottom of the wells.  Things to note: 
1) the wavefunction as written
above is not pure real -- it has an overall phase which is constant in space (independent of x), but depends on k.  rather than adjust the wavefunction to be pure real, we can just remove the overall phase explicitly. hence the factors of (E^(-I Arg[ψbot[ko,0][0]])) in each Plot[] below.
2) We're also plotting the potential, with opacity fixed by the (square root of the) strength of the potential.
*)
Show[
Plot[Re[(E^(-I Arg[ψbot[ko,0][0]]))ψbot[ko,0][x]],{x,-1,1}],
Graphics[{Black,Opacity[1],Line[{{-1,0},{1,0}}],PlotStyle->Thick}],
Graphics[{Opacity[(Sqrt[(go/gmax)])-(Sqrt[(gmin/gmax)])],TwoWells,PlotStyle->Thick}],TwoWellOptions
],
Show[
Plot[Re[-(E^(-I Arg[ψbot[k1,π][0]]))ψbot[k1,π][x]],{x,-1,1}],
Graphics[{Black,Opacity[1],Line[{{-1,0},{1,0}}],PlotStyle->Thick}],
Graphics[{Opacity[(Sqrt[(go/gmax)])-(Sqrt[(gmin/gmax)])],TwoWells,PlotStyle->Thick}],TwoWellOptions
],
Show[
Plot[Re[(E^(-I Arg[ψbot[k2,0][0]]))ψbot[k2,0][x]],{x,-1,1}],
Graphics[{Black,Opacity[1],Line[{{-1,0},{1,0}}],PlotStyle->Thick}],
Graphics[{Opacity[(Sqrt[(go/gmax)])-(Sqrt[(gmin/gmax)])],TwoWells,PlotStyle->Thick}],TwoWellOptions
],
Show[
Plot[Re[-(E^(-I Arg[ψbot[k3,π][0]]))ψbot[k3,π][x]],{x,-1,1}],
Graphics[{Black,Opacity[1],Line[{{-1,0},{1,0}}],PlotStyle->Thick}],
Graphics[{Opacity[(Sqrt[(go/gmax)])-(Sqrt[(gmin/gmax)])],TwoWells,PlotStyle->Thick}],TwoWellOptions
],
(* 
Now let's plot the probability density of the state at the bottom of the first few bands,
including the prob density of the top of the band and of the top of the prev band, too 
*)
Show[
Plot[{Abs[ψbot[ko,0][x]],1,Abs[ψtop[1][x]]},{x,-(1/2),(1/2)},ImageSize->{dpi*3*GoldenRatio,dpi*3},PlotRange->{{-(1/2),(1/2)},{-0.2,1.2}},Axes->False,Frame->True],
Graphics[{Black,Opacity[1],Line[{{-1/2,0},{1/2,0}}],PlotStyle->Thick}],
Graphics[{Opacity[(Sqrt[(go/gmax)])-(Sqrt[(gmin/gmax)])],Line[{{0,0},{0,1.2}}],PlotStyle->Thick}]
],
Show[
Plot[{Abs[ψbot[k1,π][x]],Abs[ψtopm[1][x]],Abs[ψtopm[2][x]]},{x,-(1/2),(1/2)},ImageSize->{dpi*3*GoldenRatio,dpi*3},PlotRange->{{-(1/2),(1/2)},{-0.2,1.2}},Axes->False,Frame->True],
Graphics[{Black,Opacity[1],Line[{{-1/2,0},{1/2,0}}],PlotStyle->Thick}],
Graphics[{Opacity[(Sqrt[(go/gmax)])-(Sqrt[(gmin/gmax)])],Line[{{0,0},{0,1.2}}],PlotStyle->Thick}]
],
Show[
Plot[{Abs[ψbot[k2,0][x]],Abs[ψtop[2][x]],Abs[ψtop[3][x]]},{x,-(1/2),(1/2)},ImageSize->{dpi*3*GoldenRatio,dpi*3},PlotRange->{{-(1/2),(1/2)},{-0.2,1.2}},Axes->False,Frame->True],
Graphics[{Black,Opacity[1],Line[{{-1/2,0},{1/2,0}}],PlotStyle->Thick}],
Graphics[{Opacity[(Sqrt[(go/gmax)])-(Sqrt[(gmin/gmax)])],Line[{{0,0},{0,1.2}}],PlotStyle->Thick}]
]
Show[
Plot[{Abs[ψbot[k3,π][x]],Abs[ψtopm[3][x]],Abs[ψtopm[4][x]]},{x,-(1/2),(1/2)},ImageSize->{dpi*3*GoldenRatio,dpi*3},PlotRange->{{-(1/2),(1/2)},{-0.2,1.2}},Axes->False,Frame->True],
Graphics[{Black,Opacity[1],Line[{{-1/2,0},{1/2,0}}],PlotStyle->Thick}],
Graphics[{Opacity[(Sqrt[(go/gmax)])-(Sqrt[(gmin/gmax)])],Line[{{0,0},{0,1.2}}],PlotStyle->Thick}]
]
},{go,gmin,gmax}];

PlotBandEdgeWFs : A macro to plot the band edge eigenfunctions over 10 wells as a function of barrier strength

Note that both PlotBandEdges and PlotBandEdgeWFs use the same global vaiable, so only run one at once!!!!!!

TenWells=Line[{{0,1.2},{0,0},{1,0},{1,1.2},{1,0},{2,0},{2,1.2},
{2,0},{3,0},{3,1.2},{3,0},{4,0},{4,1.2},{4,0},{5,0},{5,1.2},{5,0},{6,0},{6,1.2},{6,0},
{7,0},{7,1.2},{7,0},{8,0},{8,1.2},{8,0},{9,0},{9,1.2},{9,0},{10,0},{10,1.2}}];

TenWellOptions := {PlotRange->{{0,10},{-1.2,1.2}},Axes->False,Frame->True,ImageSize->{dpi*2*AR,dpi*2},AspectRatio->(1/AR)};
PlotBandEdgeWFs:=Manipulate[
(*First, find the *)
ko=k/.FindRoot[Cos[k]+(((go^2)+(go^3))/k)Sin[k]==1,{k,ko}]; (* qL = 0 *)
k1=k/.FindRoot[Cos[k]+(((go^2)+(go^3))/k)Sin[k]==-1,{k,π+ko}];(* qL = π *)
k2=k/.FindRoot[Cos[k]+(((go^2)+(go^3))/k)Sin[k]==1,{k,2*π+ko}];(* qL = 2π ~ 0 *)
k3=k/.FindRoot[Cos[k]+(((go^2)+(go^3))/k)Sin[k]==-1,{k,3*π+ko}];(* qL = 3π ~ π *)
{Show[
Plot[Re[(E^(-I Arg[ψbot[ko,0][0]]))ψbot[ko,0][x]],{x,0,10}],
Graphics[{Opacity[(Sqrt[(go/gmax)])-(Sqrt[(gmin/gmax)])],TenWells,PlotStyle->Thick}],
TenWellOptions
],
Show[
Plot[Re[-(E^(-I Arg[ψbot[k1,π][0]]))ψbot[k1,π][x]],{x,0,10}],
Graphics[{Opacity[(Sqrt[(go/gmax)])-(Sqrt[(gmin/gmax)])],TenWells,PlotStyle->Thick}],
TenWellOptions
],
Show[
Plot[Re[(E^(-I Arg[ψbot[k2,0][0]]))ψbot[k2,0][x]],{x,0,10}],
Graphics[{Opacity[(Sqrt[(go/gmax)])-(Sqrt[(gmin/gmax)])],TenWells,PlotStyle->Thick}],
TenWellOptions
],
Show[
Plot[Re[-(E^(-I Arg[ψbot[k3,π][0]]))ψbot[k3,π][x]],{x,0,10}],
Graphics[{Opacity[(Sqrt[(go/gmax)])-(Sqrt[(gmin/gmax)])],TenWells,PlotStyle->Thick}],
TenWellOptions
]},{go,0.01,5.5}];

ko=0.01;
PlotBandEdges
(*PlotBandEdgeWFs*)