plotting - how to plot the dispersion curve in first brillouin zone?

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

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}]