Potencial periódico

Modelo de Kronig-Penney

Consideremos el movimiento de una partícula en un potencial periódico de periodo l=a+b, formado por pozos de potencial de anchura a y profundidad V0 y barreras de potencial de anchura b.

La función de onda tiene la forma

Ψ(x)=u(x)exp(ikx)

Obtenemos la imagen de las funciones de onda considerando que u(x) se asemeja a la función de onda de los átomos aislados (de un pozo de potencial)

y reemplazando exp(ikx) por las funciones de onda de una partícula libre en una caja de potencial. Esto es lo que hemos observado al visualizar las funciones de onda de los primeros niveles de energía de un sistema de pozos de potencial.

Como el potencial es periódico con periodo l=a+b, suma de la anchura a del pozo y de la separación b entre pozos, se deberá cumplir que

u(x+l)=u(x)

Ambas expresiones constituyen el teorema de Bloch.

En la figura, se muestra tres regiones en las que vamos a obtener la solución de la ecuación de Schrödinger.

El punto x en la región 3 se corresponde con el punto x-l en la región 1, de modo que u(x)=u(x-l).

A 1 exp( i k 1 (xl) )+ B 1 exp( i k 1 (xl) )=u(x)exp( ik(xl) )

Despejando u(x) e introduciendo dicha expresión en la función de onda en la tercera región.

Ψ 3 (x)=exp(ikl)( A 1 exp( i k 1 (xl) )+ B 1 exp( i k 1 (xl) ) )

Escribiremos ahora las condiciones de continuidad de la función de onda y de su derivada primera en los puntos x=0, y x=a.

Tenemos un sistema homogéneo de cuatro ecuaciones con cuatro incógnitas, el determinante de los coeficientes debe ser cero. Para el caso en el que E<V0, k1 es una cantidad imaginaria, llamemos k1=ik3.

cos(kl)=cos( k 2 a)cosh( k 3 b) + k 3 2 k 2 2 2 k 2 k 3 sin( k 2 a)sinh( k 3 b)

Ecuación que nos da la relación entre la energía E y el número de onda k y que representaremos en la ventana gráfica. Ya que el módulo del coseno no puede ser mayor que la unidad, obtenemos así la condición impuesta a k3 y a k2 y por tanto a la energía E.

1cos( k 2 a)cosh( k 3 b) + k 3 2 k 2 2 2 k 2 k 3 sin( k 2 a)sinh( k 3 b)1

Esta condición define las bandas de energía permitidas.

Puede observarse que la energía presenta una discontinuidad, para números de onda k múltiplos de π/l, donde l=a+b es el periodo de la red lineal. Estos valores representan fronteras entre zonas de Brillouin contiguas. Los intervalos de existencia de cada una de las zonas medidas en el eje vertical representan las bandas de energía, que se muestran como rectángulos de color azul en la parte derecha de la ventana.

Cálculo de las bandas de energía

function res =kp(a, b, H, cte, E)  
   q=sqrt(E);
   k=sqrt(H-E);
   res=cos(q*a)*cosh(k*b)+(k^2-q^2)/(2*q*k)*sin(q*a)*sinh(k*b)+cte;
end
a=2;  %anchura del pozo
b=1; %separación entre pozos
H=5.0; %profundidad del pozo

f=@(E) kp(a,b,H, -1.0, E);
xb=buscar_intervalos(f,0.1,H-0.1,10);
nb=size(xb);
for i=1:nb(1)
    bandas(i)=fzero(f,[xb(i,1) xb(i,2)]);
end

f=@(E) kp(a,b,H, 1.0, E);
xb=buscar_intervalos(f,0.1,H-0.1,10);
nc=size(xb);
for i=1:nc(1)
    bandas(i+nb(1))=fzero(f,[xb(i,1) xb(i,2)]);
end
bandas=sort(bandas);  %ordena de menor a mayor
if rem(length(bandas),2)==1
    bandas=[bandas H];
end
hold on
%representa las bandas de energía mediante rectángulos de color azul
for j=1:2:length(bandas)
    xx=[3 3 4 4];
    yy=[bandas(j), bandas(j+1), bandas(j+1), bandas(j)];
    fill(xx,yy,'b')
    text(4, bandas(j),num2str(bandas(j)),
'VerticalAlignment','top', 'HorizontalAlignment','right');
    text(4, bandas(j+1),num2str(bandas(j+1)),
'VerticalAlignment','bottom', 'HorizontalAlignment','right');
end
.....

La función buscar_intervalos busca los intervalos donde la función f(E) cambia de signo. Se ha definido en la página Atomo, molécula... sólido lineal

Zonas de Brillouin

Completamos el código del script para representar las zonas de Brillouin

....
f=@(E) kp(a,b,H, 0.0, E);

for j=1:2:length(bandas)
    E=bandas(j):0.01:bandas(j+1);
    k=zeros(length(E));
    k(1)=(j-1)/2;
    for i=2:length(E)          
        y=f(E(i));
        if(abs(y)<1)
            angulo=acos(y);
            if rem((j-1)/2,2)==1
                k(i)=(j-1)/2+1-angulo/pi;
            else
                k(i)=(j-1)/2+angulo/pi;
            end
        end
    end
    plot(k,E,'r') %simetría
    plot(-k,E,'r')
end
line([0,0],[0,H],'color','k')

hold off
xlabel('k')
ylabel('E')
title('Modelo de Kronig-Penney')