how to plot the dispersion curve in first brillouin zone?

Question

Those who know some solid state physics should know what the first brillouin zone is. How do I plot the dispersion relation in the 1st brillouin zone so that the curves can 'fold back'?

For instance the code I have now is:

OmegaP=8; 
epM[Omega_]:=1-OmegaP^2/Omega^^2; 
epD[Omega_]:= 1; 

e2alpha[Omega_]:=(epM[Omega]-epD[Omega])/(epM[Omega]+epD[Omega])    
ContourPlot[{Exp[0.5*ak]==e2alpha[Omega],Exp[0.5*ak]==-e2alpha[Omega]},{ak,0,0.5},{Omega,0.05,8}];

I have rescaled the axis so that 1st Brillouin zone is 0-0.5, 2nd Brillouin zone is 0.5-1.0, and so on. The above above only generates the lower part of the whole dispersion relation. If I extend the range of ak to, say {ak,0,2.5}, then I get a lot more missing curves, but how to I 'fold' the part from 0.5 to 2.5 into the 1st Brillouin zone?

PS: Please click here (last page, Fig.16) for an image showing what the '1st Brillouin zone' is.

Answer 1

The "folding in" of the reduced zone scheme is achieved by translating the center of the energy function, placing one center on each brillouin zone edge:

Plot[
 {x^2, (x - Pi)^2, (x + Pi)^2, (x - 2 Pi)^2, (x + 2 Pi)^2},
 {x, -Pi/2, Pi/2},
 AspectRatio -> GoldenRatio,
 PlotStyle -> Black
 ]

I plotted here a function proportional to $x^2$ as would be the case for free electrons. I used AspectRatio only to compensate for not having the right constant of proportionality, $\frac{\hbar^2}{2m}$. The centers were taken to be $\pm \pi n$ but should be $\pm \frac{\pi n}{a}$ for the actual case, where $a$ is the lattice constant.

To plot an arbitrary number of energy bands without having to type in the functions manually you can generate them using Table:

Plot[
 Flatten@Table[{(x + Pi n)^2, (x - Pi n)^2}, {n, 0, 10}],
 {x, -Pi/2, Pi/2},
 AspectRatio -> GoldenRatio,
 PlotStyle -> Black
 ]

For your dispersion relation, maybe something like this could work:

eqns = Flatten@Table[{
     Exp[0.5*(ak - n Pi)] == e2alpha[Omega],
     Exp[0.5*(ak - n Pi)] == -e2alpha[Omega],
     Exp[0.5*(ak + n Pi)] == e2alpha[Omega],
     Exp[0.5*(ak + n Pi)] == -e2alpha[Omega]
     }, {n, 0, 5}];

ContourPlot[Evaluate[eqns], {ak, 0, 0.5}, {Omega, 0.05, 8}]