Problem plotting Kronig Penney Model dispersion curves

[maverick280857]

To the moderator: I'm not sure if this should go here or in the Computational Physics forum. Please shift it there if you think that's the appropriate place for it.

Hi everyone

Merry Christmas!

I'm writing a computer program in C, to explicitly compute the band structures for a 1D crystal modeled as an infinite array of rectangular potential barriers of width 2a and distance between nearest walls of two successive barriers 2b, i.e. the Kronig Penney Model.

My program correctly plots the energy versus (normalized) wavenumber graph in the extended wavenumber scheme, but when I try to superimpose the free particle parabola on this graph (i.e. the parabolic graph of energy versus wavenumber for a free particle), it is not found to intersect the band curves at points of the form n\pi where n is an integer. (I call gnuplot from within my C program to plot this curve.)

The plot generated by the program is as shown in the attachment.

I think this is because of some numerical approximations and propagation errors, but I haven't been able to fix it. As you can see from the plot, the two curves do intersect at \pi (or very nearly so) but that's just it.

I would appreciate any ideas or suggestions. I can post my algorithm/code here if you wish.

Thanks in advance,
Cheers
Vivek.

 

[mercurial]

Hi Vivek,

I'd love to help but need a bit more information about how you're calculating energy from the wavenumber. Are you using an explicit formula? I'm not a physicist by training, so please excuse me if there is a closed form for the energy as a function of parameters that I just don't know about.

It may also be helpful if you post the code.

Jerry

 

[maverick280857]

Ok this is how I do it. The condition for existence of bands is

\cos\chi = f(\epsilon)

where \chi is proportional to the wavenumber \kappa, and \epsilon is the ratio of the energy of the particle E to the height of the potential barrier V_{0}:

\epsilon = \frac{E}{V_{0}}

Now, I iterate over \epsilon. For every value of \epsilon, I compute f(\epsilon). If f(\epsilon) \in [-1, 1], then a real value of \chi exists and I find \chi = \cos^{-1}[f(\epsilon)]. (V_{0} is fixed)

By definition of the arccos function, this value of \chi lies in the principal value interval [0, \pi], so it has to be adjusted according to the band it belongs to, for purposes of plotting the curve in the extended wavenumber scheme (otherwise everything gets mapped to the principal value interval as \chi is not monotonically increasing).

The parabola is \epsilon = C\chi^2 (the proportionality constant C is related to the height of the well, the mass of the particle and the ratio of the a/b, where a and b are as defined in my first post.)

 

[mercurial]

Hi Maverick,

Thanks for the clarification. If I understand it right, the pairs (\chi, \epsilon) satisfying the equation

\cos \chi = f(\epsilon)

should lie on the parabola \epsilon = C\chi^2. So far so good? The jpg you attached indeed shows a discrepancy much higher than one would expect due to numerical error for computing arc-cosine of f(\epsilon). The most obvious question is then: how do you compute f(\epsilon)?

If this depends on, say, the eigenvalue from some boundary-value problem that you in turn find via the shooting method, some matrix calculation, etc., the problem may lie in that your computed value of f(\epsilon) is too far off the mark. To correct this you may need to work on your numerical method, find a better discretization, etc.

Please do supply more information regarding this point.

Jerry

 

[maverick280857]

Hi Maverick,

Thanks for the clarification. If I understand it right, the pairs (\chi, \epsilon) satisfying the equation

\cos \chi = f(\epsilon)

should lie on the parabola \epsilon = C\chi^2. So far so good?

No, that's not right. These are two independent curves:

C_{1}: (\chi, \epsilon) = (\cos^{-1}[f(\epsilon)], \epsilon)

C_{2}: (\chi, \epsilon) = (\chi, C\chi^2)\qquad{ \mbox{ (C = constant) } }[/itex]<br /> <br /> It is only at points of the form \chi = n\pi where n is an integer, do the two curves intersect. This is as far as the theory tells us.<br /> <br /> <blockquote data-attributes="" data-quote="" data-source="" class="bbCodeBlock bbCodeBlock--expandable bbCodeBlock--quote js-expandWatch"> <div class="bbCodeBlock-content"> <div class="bbCodeBlock-expandContent js-expandContent "> The most obvious question is then: how do you compute f(\epsilon)? </div> </div> </blockquote><br /> f(\epsilon) has piecewise continuous definitions for 0 \leq \epsilon \leq 1 and \epsilon \geq 1. It is generally of the form<br /> <br /> f(\epsilon) = \cosh(k_{1}(1-\epsilon)^{1/2})cos(k_{2}\epsilon^{1/2}) + \frac{1-2\epsilon}{2\epsilon^{1/2}(1-\epsilon)^{1/2}}\sinh(k_{1}(1-\epsilon)^{1/2})sin(k_{2}\epsilon^{1/2})<br /> <br /> where the k_{i}&#039;s are constants. When \epsilon &amp;gt; 1 the expression for f(\epsilon) is obtained by replacing the square root of 1-\epsilon by the square root of \epsilon-1 and cosh replaced by cos (I will post the exact expressions later as I don&#039;t have them before me right now--but they&#039;re just the boundary conditions for the wavefunction).

 

[maverick280857]

I just had a thought. When these potential barriers become delta functions in the limit a --> 0 and V0 --> infinity with V0*2a = constant, the matching condition reduces to a simpler form and it is possible to show analytically that the two curves C1 and C2 intersect at points of the form \chi = n\pi.

But does this have to be true for a rectangular potential barrier array as well? In other words, do the band curves and the free particle parabola have to intersect each other at the end of every band?

This question is more appropriate for the Quantum Physics forum I guess.

 

[mercurial]

Maybe this thread is helping you talk yourself through the points until you locate the source of the problem. I really should read up on the Kronig Penny model myself. Do the parameters a, V_0 appear in the constants k_i in your expression for f(\epsilon), or do their values "fall out" in the limit subject to replacing the potential barriers with delta functions?

 

[maverick280857]

The full function is:

<br /> f(\epsilon) = \left\{ \begin{array}{ccc}<br /> \cosh[rA] + \frac{A}{2}\sinh[rA] &amp; \mbox { for } \epsilon = 0\\<br /> <br /> \cosh[rA(1-\epsilon)^{1/2}]\cos[A\epsilon^{1/2}] + \frac{(1-2\epsilon)}{2\epsilon^{1/2}(1-\epsilon)^{1/2}}\sinh[rA(1-\epsilon)^{1/2}]\sin[A\epsilon^{1/2}] &amp; \mbox { for } 0 &lt; \epsilon &lt; 1\\<br /> <br /> \cos[rA(\epsilon-1)^{1/2}]\cos[A\epsilon^{1/2}] + \frac{(1-2\epsilon)}{2\epsilon^{1/2}(\epsilon-1)^{1/2}}\sin[rA(\epsilon-1)^{1/2}]\sin[A\epsilon^{1/2}] &amp; \mbox { for } \epsilon &gt; 1\\<br /> <br /> <br /> \end{array}<br /> \right.<br />

\chi &amp; = &amp; \kappa \times 2b(1+r)
\epsilon &amp; = &amp; \frac{E}{V_{0}}
r &amp; = &amp; \frac{a}{b}
A &amp; = &amp; 2b\sqrt{\frac{2mV_{0}}{\hbar^2}}

The matching condition is

f(\epsilon) = \cos[\kappa\times 2b(1+r)] = \cos \chi