0
$\begingroup$

The goal is to display the relation between K and E in the following equation:

$$\cos(K) = 5 \text{sinc}( 5.12\sqrt{E} ) - \cos(5.12\sqrt{E} )$$

The intended result is:

enter image description here

However, using the following code:

k = linspace(0, 3*pi, 99);  % 0 < k < 3pi , crystal momentum.
E = linspace(0.1,  3.5, 99);% Energy.
[X, Y] = meshgrid(k, E);

P = 5;                      % 2mU_0b/hbar^2.
Z = cos(X)  - ( P * sinc( 5.12 * sqrt(Y) ) )  - cos( 5.12 * sqrt(Y) );

contour(X, Y, Z);

xlabel('k')
ylabel('E(k)')
legend('cos(K) - P * sin(E) - cos(E)')
title('Dispersion relation for periodic potential in extended zone.')

what I currently get unfortunately bears no resemblance at all:

enter image description here

I just don't see how to separate the graph into three (discontinuous) regions, $K = [0, \pi], [\pi, 2\pi], [2\pi, 3\pi]$, and what value of $E$ should be used for each of these regions.

How to plot the equation the same way as in the first graph?


The first plot is from here, where they do some modular arithmetic trick by adding $i*floor(k ./ \pi)$ on the l.h.s. and $i*floor(5.12*\sqrt{E} ./ \pi)$ on the r.h.s. There is something similar here.

$\endgroup$
  • 1
    $\begingroup$ I tried to find the solution using bisection and am now convinced that the plot and the closed-form dispersion relation do not correspond to each other. Highlight the zero level-set in your contour plot and you'll see why, e.g., v = [0,0]; contour(X,Y,Z,v). $\endgroup$ – Biswajit Banerjee Sep 14 '18 at 22:59
3
$\begingroup$

First, you originally wrote that your equation is $\cos(K)=5 \text{ sinc}(5.12\sqrt{E})- \cos(5.12\sqrt{E})$, but you clearly meant $\cos(K)=5 \text{ sinc}(5.12\sqrt{E}) + \cos(5.12\sqrt{E})$

Second, the $\text{sinc}()$ function has differing conventions. The other question is based on Mathematica, which uses the convention $\text{sinc}(x) \equiv \frac{\sin x}{x}$, while Matlab/Octave use the convention $\text{sinc}(x) \equiv \frac{\sin (\pi x)}{\pi x}$. (To be precise of course we need to separately define $\text{sinc}(0)\equiv 1$.) This is why your first contour plot shows "no resemblance at all".

Compensating for these differing conventions and altering your code slightly so that it only shows the contours for Z=0 (as suggested in a comment by Biswajit Banerjee) yields

k = linspace(0, 3*pi, 99);  % 0 < k < 3pi , crystal momentum.
E = linspace(0.1,  3.5, 99);% Energy.
[X, Y] = meshgrid(k, E);

P = 5;                      % 2mU_0b/hbar^2.
Z = cos(X)  - ( P * pi * sinc( 5.12 * sqrt(Y)/pi ) )  - cos( 5.12 * sqrt(Y) );
contour(X, Y, Z,[0,0],'linewidth',2,'linecolor','k');
for i=1:3
  line(i*[pi,pi],[0,3.5],'color','k','linestyle','--')
end

xlabel('k')
ylabel('E(k)')
legend('cos(K) - P * sin(E) - cos(E)')
title('Dispersion relation for periodic potential in extended zone.')

which produces the plot

contour plot

The best way I could come up with to isolate sections of these curves is to simply run contour for the separate regions and plot them together. The following code accomplishes this and generalizes it to N sections.

P = 5;                      % 2mU_0b/hbar^2.
a=5.12;
N=3;

hold on;


for i=1:N
    k = linspace((i-1)*pi, i*pi, 1000);  % 0 < k < 3pi , crystal momentum.
    E = linspace(((i-1)*pi/a)^2,  (i*pi/a)^2, 1000);% Energy.
    [X, Y] = meshgrid(k, E);

    Z = cos(X)  - ( P * pi * sinc( a * sqrt(Y)/pi ) )  - cos( a * sqrt(Y) );
    contour(X, Y, Z,[0;0],'linewidth',2,'linecolor','k');

    line(i*[pi,pi],[0,(N*pi/a)^2+0.2],'color','k','linestyle','--')
end

xlabel('k')
xlim([0,N*pi+0.2])
ylim([0,(N*pi/a)^2+0.2])
ylabel('E(k)')

contour plot sectioned

$\endgroup$
0
$\begingroup$

Here is something near to the intended goal:

enter image description here

But the code is ugly:

k = linspace(0,  3*pi, 100);   

hbar = 1;
m = 1;
E =  ( (hbar^2) .* (k.^2) ) ./  (2*m);

E(33:66) = E(33:66) + 3;
E(66:100) = E(66:100) + 6;

E(33) = E(32);
E(66) = E(65);

E(34) = NaN;
E(67) = NaN;

hold on
plot(k,  E, [pi pi], [0 60], '--k', [2*pi 2*pi], [0 60], '--k');
plot(k(34), E(33), 'ob', 'MarkerFaceColor','b',...
        k(67), E(66), 'ob',  'MarkerFaceColor', 'b',...
        k(34), E(35), 'ob', 'MarkerFaceColor', 'b',...
        k(67), E(68), 'ob' , 'MarkerFaceColor', 'b');
hold off

xlabel('k')
ylabel('E(k)')
legend('E = hbar^2 k^2 / 2m', 'x = pi', 'x = 2pi')
title('Dispersion relation for periodic potential in extended zone.')

txt = '<- Band Gap';
text(k(35), E(33) +2, txt)

$\endgroup$

Your Answer

By clicking "Post Your Answer", you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.