El pozo de potencial de forma rectangular

La cuantización de la energía es uno de los conceptos más importantes de la Mecánica Cuántica, ya que explica las propiedades de los átomos que constituyen los componentes básicos de la materia.

Para calcular los niveles de energía, es necesario resolver una ecuación diferencial de segundo orden, la ecuación de Schrödinger para la función potencial especificada, que en muchos casos carece de solución analítica sencilla. Por simplicidad, elegiremos como modelo de átomo un pozo de potencial.

La ecuación de Schrödinger independiente del tiempo en una región unidimensional cuya energía potencial viene descrita por la función Ep(x) es

2 2m d 2 Ψ d x 2 + E p (x)Ψ=EΨ

Donde E es la energía total de la partícula de masa m

E= p 2 2m + E p (x)

La solución de la ecuación de Schrödinger Ψ(x) se denomina función de onda. La función de onda se normaliza de modo que

+ | Ψ(x) | 2 dx=1

La simetría de la función potencial Ep(x) hace que las funciones de onda puedan ser:

Partícula libre

El caso más simple es el de una partícula libre, la energía potencial Ep(x)=0. La ecuación de Schrödinger se escribe

d 2 Ψ d x 2 + 2mE 2 Ψ=0

Su solución es

Ψ(x)=Aexp(iqx)+Bexp(iqx) q 2 = 2mE 2

Pozo de potencial de profundidad infinita

Consideremos una partícula obligada a moverse en una región entre x=-a y x=a, tal como una molécula en una caja, un electrón libre en un trozo de metal, etc. Si la energía cinética del electrón es pequeña comparada con la altura de la barrera de potencial, el electrón se podrá mover libremente a través del metal pero no podrá escapar de él.

Podemos representar estas situaciones físicas, por un potencial rectangular de altura infinita. Tenemos que Ep(x)=0 para -a<x<a, ya que la partícula se mueve libremente en esta región, fuera de esta región la energía potencial se hace infinita. Entonces, cualquiera que sea el valor de le energía E de la partícula, ésta no puede estar a la izquierda de x=-a, ni a la derecha de x=a. La función de onda en dichas regiones debe de ser nula.

La ecuación de Schrödinger en la región 0<x<a donde Ep(x)=0 se escribe

d 2 Ψ d x 2 + 2mE 2 Ψ=0

Su solución es

Ψ(x)=A e iqx +B e iqx q 2 = 2mE 2

Niveles de energía

Las condiciones de contorno requieren que Ψ(x)=0 en x=-a y también, que Ψ(x)=0 en x=a.

Ψ(a)=Aexp(iqa)+Bexp(iqa)=0 Ψ(a)=Aexp(iqa)+Bexp(iqa)=0 exp(2iqa)exp(iqa)=0 sin(2qa)=02qa=nπ

donde n es un número entero. La energía de la partícula será

E= n 2 π 2 2 8m a 2

Si E1 es la energía del primer nivel (n=1) la energía de los sucesivos niveles es 4E1, 9E1, 16E1... Concluimos que la partícula no puede tener una energía arbitraria, sino valores concretos, decimos que la energía de la partícula está cuantizada.

Funciones de onda

La simetría de la función potencial Ep(x) hace que los funciones de onda sean

En ambos casos el factor de escala, A se determina haciendo que

| Ψ(x) | 2 =1 a a | Ψ(x) | 2 =1 a a 4 A 2 cos 2 (qx)dx= 2 A 2 ( 2a+ sin(2qa) q ) a a 4 A 2 sin 2 (qx)dx= 2 A 2 ( 2a sin(2qa) q )

Las funciones de onda se parecen a los modos de vibración de una cuerda sujeta por ambos extremos o también denominadas ondas estacionarias.

Las funciones de onda dibujadas en color rojo, son simétricas. Las funciones de onda dibujadas en color azul, son antisimétricas

%caja de potencial
a=0.5;  %anchura de la caja 2a
Emax=18; %disminuir este valor al incrementar la anchura a
hold on
%potencial
xx=[-1 -1 -a -a a a 1 1];
yy=[0 Emax Emax 0 0 Emax Emax 0];
fill(xx,yy, [0.5 0.5 0.5])
%niveles de nergía y funciones de onda
for n=1:2:4
    %par, simétrico
    E=(n/(2*a))^2;
    line([-a a],[E E], 'color','k')
    q=n*pi/(2*a);
    suma=2*(2*a+sin(2*q*a)/q);
    x=-a:0.02:a;
    y=2*cos(q*x)/sqrt(suma);
    plot(x,E+y,'r');
    
    %impar, antisimétrico   
    E=((n+1)/(2*a))^2;
    line([-a a],[E E], 'color','k')
    q=(n+1)*pi/(2*a);
    suma=2*(2*a-sin(2*q*a)/q);
    x=-a:0.05:a;
    y=2*sin(q*x)/sqrt(suma);
    plot(x,E+y,'b');
end
hold off
xlim([-1,1])
xlabel('x');
ylabel('V')
title('Caja de potencial')

El pozo de potencial de forma rectangular

Sea un pozo de potencial de anchura 2a y de altura V0. Situamos el origen en el centro del pozo, tal como se muestra en la figura. Las ecuaciones de Schrödinger en las regiones (1) y (2) se escriben

2 2m d 2 Ψ d x 2 =EΨ0xa 2 2m d 2 Ψ d x 2 + V 0 Ψ=EΨx>a

Las soluciones de estas ecuaciones son, respectivamente.

Ψ 1 ( x )= A 1 exp( i q 1 x )+ B 1 exp( i q 1 x ) q 1 = 2mE 2 Ψ 2 ( x )= A 2 exp( i q 2 x )+ B 2 exp( i q 2 x ) q 2 = 2m(E V 0 ) 2

q1 es real y q2 es imaginario. E<V0 es la energía de un nivel, si la solución de le ecuación de Schrödinger, Ψ(x) tiende asintóticamente a cero para grandes valores de x.

lim x± Ψ(x)=0

Ψ2(x) debe tender a cero cuando x se hace grande, para ello B2 tiene que ser cero.

Las condiciones de continuidad de la función de onda Ψ(x) y su derivada primera en la frontera x=a entre las dos regiones de distinto potencial, constituyen un par de ecuaciones que relacionan A1 y B1 con A2.

x=a{ A 1 exp( i q 1 a )+ B 1 exp( i q 1 a )= A 2 exp( i q 2 a ) q 1 A 1 exp( i q 1 a ) q 1 B 1 exp( i q 1 a )= q 2 A 2 exp( i q 2 a ) A 1 = q 1 + q 2 2 q 1 A 2 exp( i q 2 a )exp( i q 1 a ) B 1 = q 1 q 2 2 q 1 A 2 exp( i q 2 a )exp( i q 1 a )

Este último parámetro, determina la escala vertical de la función de onda Ψ(x), y se puede obtener a partir de la condición de normalización

+ | Ψ(x) | 2 dx=1

Niveles de energía

La simetría de la función potencial Ep(x) hace que las funciones de onda sean:

Las raíces de las dos ecuaciones nos dan los niveles de energía de la partícula en el pozo de potencial.

Establecemos un sistema de unidades en el que ℏ=m=1. La anchura del potencial es 2a=2·2, la altura V0=5

function pozo
    a=2;   %la anchura del pozo de potencial es 2a
    V0=5;  %altura del pozo de potencial 
    f=@(E) sqrt(2*(V0-E)).*cos(sqrt(2*E)*a)-sqrt(2*E).*sin(sqrt(2*E)*a); %par
    g=@(E) sqrt(2*E).*cos(sqrt(2*E)*a)+sqrt(2*(V0-E)).*sin(sqrt(2*E)*a); %impar
    xx=linspace(0,V0-0.5,10);
    impar=raices(g,xx);
    par=raices(f,xx);
    disp(par)
    disp(impar)
    
     
    function r = raices(f, x)
        y=f(x);
        indices=find(y(1:end-1).*y(2:end)<0);
        r=zeros(1,length(indices));
        for m=1:length(indices)
            r(m)=fzero(f, [x(indices(m)), x(indices(m)+1)]);
        end
    end
end
0.2295    2.0230
0.9116    3.5005

La función raices calcula las raíces de las ecuaciones transcendentes empleando a función fzero de MATLAB.

Funciones de onda

Representamos las funciones de onda para cada uno de los niveles de energía

Ψ 1 ( x )= A 1 exp( i q 1 x )+ B 1 exp( i q 1 x ) Ψ 2 ( x )= A 2 exp( i q 2 x )

Subtituyendo q1=q y q2=ik

La función de onda simétrica, A1=B1 es

Ψ 1 ( x )=2 A 1 cos( qx )0xa Ψ 2 ( x )= A 2 exp( kx )x>a

La continuidad en x=a establece las relaciones entre A1 con A2 que actúa como factor de escala

A 1 = A 2 qcos(qa)+ksin(qa) 2q exp(ka)

La función de onda antisimétrica, A1=-B1 es

Ψ 1 ( x )=2i A 1 sin( qx )0xa Ψ 2 ( x )= A 2 exp( kx )x>a A 1 =i A 2 kcos(qa)qsin(qa) 2q exp(ka)

Determinamos el valor de A2 para que la integral

| Ψ(x) | 2 dx=1 0 | Ψ(x) | 2 dx= 0 a | Ψ 1 (x) | 2 dx+ a | Ψ 2 (x) | 2 dx= 1 2 0 a 4 A 1 2 cos 2 (qx)dx+ a A 2 2 exp(2kx)dx=2 A 1 2 ( a+ sin(2qa) 2q ) + A 2 2 2k exp(2ka) 0 a 4 A 1 2 sin 2 (qx)dx+ a A 2 2 exp(2kx)dx=2 A 1 2 ( a sin(2qa) 2q ) + A 2 2 2k exp(2ka)

Completamos el script anterior, representando el potencial V(x), los tres niveles de energía y las correspondientes funciones de onda, simétricas (color rojo), antisimétricas (azul)

function pozo_1
    a=2;   %la anchura del pozo de potencial es 2a
    V0=5;  %altura del pozo de potencial 
    f=@(E) sqrt(2*(V0-E)).*cos(sqrt(2*E)*a)-sqrt(2*E).*sin(sqrt(2*E)*a); %par
    g=@(E) sqrt(2*E).*cos(sqrt(2*E)*a)+sqrt(2*(V0-E)).*sin(sqrt(2*E)*a); %impar
    xx=linspace(0,V0-0.5,10);
    impar=raices(g,xx);
    par=raices(f,xx);
    disp(par)
    disp(impar)
    
    hold on
    %pozo de potencial
    xx=[-4 -4 -a -a a a 4 4];
    yy=[0 V0 V0 0 0 V0 V0 0];
    fill(xx,yy, [0.5 0.5 0.5])
    %niveles de nergía y funciones de onda
    for i=1:length(par)
       %nivel de energía
        line([-4 4],[par(i) par(i)], 'color','k')
        text(4, par(i),num2str(par(i)),'VerticalAlignment','bottom', 
'HorizontalAlignment','right');

        %función de onda (simétrica)
        q=sqrt(2*par(i));
        k=sqrt(2*(V0-par(i)));
        A1=(q*cos(q*a)+k*sin(q*a))*exp(-k*a)/(2*q);
        suma=2*A1^2*a+A1^2*sin(2*q*a)/q+exp(-2*k*a)/(2*k);
        x=linspace(0,4,50);
        y=((2*A1*cos(q*x)).*(x<=a)+(exp(-k*x)).*(x>a))/sqrt(suma);
        plot([fliplr(-x) x], [par(i)+fliplr(y) par(i)+y],'r');
    end
    for i=1:length(impar)
       %nivel de energía
        line([-4 4],[impar(i) impar(i)], 'color','k')
        text(4, impar(i),num2str(impar(i)), 'VerticalAlignment','bottom',
 'HorizontalAlignment','right');
        %función de onda (antisimétrica)
        q=sqrt(2*impar(i));
        k=sqrt(2*(V0-impar(i)));
        A1=(k*cos(q*a)-q*sin(q*a))*exp(-k*a)/(2*q);
        suma=2*A1^2*a-A1^2*sin(2*q*a)/q+exp(-2*k*a)/(2*k);
        x=linspace(0,4,50);
        y=((2*A1*sin(q*x)).*(x<=a)-(exp(-k*x)).*(x>a))/sqrt(suma);
        plot([fliplr(-x) x], [impar(i)+fliplr(-y) impar(i)+y],'b');
    end
    ylim([0 V0+0.1]);

    hold off
    xlabel('x');
    ylabel('V')
    title('Pozo de potencial')
     
    function r = raices(f, x)
        y=f(x);
        indices=find(y(1:end-1).*y(2:end)<0);
        r=zeros(1,length(indices));
        for m=1:length(indices)
            r(m)=fzero(f, [x(indices(m)), x(indices(m)+1)]);
        end
    end
end

Pozo de potencial de forma triangular

En este apartado, resolveremos la ecuación de Schrödinger unidimensional e independiente del tiempo en el potencial

V(x)={ V 0 ,x<a V 0 x a ,ax<0 V 0 x a ,0xa V 0 ,x>a

Para un potencial simétrico, resolveremos la ecuación de Schrödinger para E<V0 en las regiones I y II

d 2 ψ I (x) d x 2 2m 2 ( V 0 x a E ) ψ I (x)=0,0xa d 2 ψ II (x) d x 2 2m 2 ( V 0 E )· ψ II (x)=0,x>a

La solución de la segunda, en la región II es

ψ II (x)=Cexp( kx )+Dexp( kx ), k 2 = 2m 2 ( V 0 E )

El coeficiente D deberá ser nulo para evitar que la función de onda se haga muy grande cuando x→∞.

En la región I, hacemos el cambio de variable z=αx+β, dz=α·dx, de modo que la ecuación diferencial se transforme

d 2 ψ I (z) d z 2 z· ψ I (z)=0, ψ I (z)=A·Ai(z)+B·Bi(z) z= ( 2m 2 V 0 a 2 ) 1/3 ( x a E V 0 )=r( x a E V 0 )

cuya solución es una combinación lineal de las funciones de Airy, Ai(x) y Bi(x).

Las posiciones x=0 y x=a, corresponden a los valores de la variable z

z 0 =r E V 0 z a =r( 1 E V 0 )

Continuidad

La función de onda y su derivada primera son continuas en x=a.

ψ I ( z a )= ψ II (a) d ψ I (x) dx | x=a = d ψ II (x) dx | x=a

La segunda condición (continuidad de la derivada primera) se expresa en términos de la variable z

d ψ I (x) dx = d ψ I (z) dz dz dx = r a d ψ I (z) dz r a d ψ I (z) dz | z= z a = d ψ II (x) dx | x=a

La continuidad en x=a, es

x=a{ A·Ai( z a )+B·Bi( z a )=Cexp(ka) r( A·Ai'( z a )+B·Bi'( z a ) )=kaCexp(ka)

El símbolo ' significa la derivada de las función de Airy Ai(z) o Bi(z) respecto de la variable z. Esta derivada se evalúa en z=za

Niveles de energía

Establecemos un sistema de unidades en el que ℏ=m=1. La anchura del potencial es 2a=2·2, la altura V0=5

function airy_7
    a=2; %la anchura del pozo es 2a
    V0=5; %profundidad
    rr=(2*V0*a^2)^(1/3); %constante
    xx=linspace(0,V0-0.1,10);
    rI=raices(@fImpar,xx);
    rP=raices(@fPar,xx);
    En=sort([rI,rP]);
    disp(En)

    function res=fImpar(E)    
        z0=-rr*E/V0;
        za=rr*(1-E/V0);
        k=sqrt(2*(V0-E));
        res=(a*k.*airy(0,za)+rr*airy(1,za)).*airy(2,z0)-airy(0,z0).*
(a*k.*airy(2,za)+rr*airy(3,za));
    end

    function res=fPar(E)    
        z0=-rr*E/V0;
        za=rr*(1-E/V0);
        k=sqrt(2*(V0-E));
        res=(a*k.*airy(0,za)+rr*airy(1,za)).*airy(3,z0)-airy(1,z0).*
(a*k.*airy(2,za)+rr*airy(3,za));
    end
    function r = raices(f, x)
        y=f(x);
        indices=find(y(1:end-1).*y(2:end)<0);
        r=zeros(1,length(indices));
        for m=1:length(indices)
            r(m)=fzero(f, [x(indices(m)), x(indices(m)+1)]);
        end
    end
end
1.4893    3.4103    4.6449

Funciones de onda

{ ψ I (x)=A·Ai(z)+B·Bi(z),0x<a,z=r( x a E V 0 ) ψ II (z)=Cexp(kx),ax<

Los coeficientes B y C están relacionados con A. El coeficiente A se calcula de modo que

0 a | ψ I (x) | 2 dx+ a | ψ II (x) | 2 dx= 1 2

por la simetría de la función de onda. La segunda integral es inmediata. La primera integral se resuelve numéricamente utilizando integral de MATLAB para ello se hace el cambio de variable de x a z, dz=(r/adx, los nuevo límites de integración son z0 y za

a r z 0 z a ( A·Ai(z)+B·Bi(z) ) 2 dz+ 1 2k exp(2ka)= 1 2

Completamos el script anterior, representando el potencial V(x), los tres niveles de energía y las correspondientes funciones de onda, simétrica y antisimétrica

function airy_8
    a=2; %anchura 2a
    V0=5; %profundidad
    rr=(2*V0*a^2)^(1/3); %constante
    xx=linspace(0,V0-0.1,10);
    rI=raices(@fImpar,xx);
    rP=raices(@fPar,xx);
    En=sort([rI,rP]);
    disp(En)
    hold on
    %potencial V(x)
    line([-a,0],[V0,0], 'color','r', 'lineWidth',1.5)
    line([0,a],[0,V0], 'color','r', 'lineWidth',1.5)
%niveles de energíafunciones de onda
    x=linspace(0,2*a,100);
    for m=1:length(En)
        z0=-rr*En(m)/V0;
        za=rr*(1-En(m)/V0);
        k=sqrt(2*(V0-En(m)));
        if rem(m,2)==1 %si es par m=1,3,5...
             line([-En(m)*a/V0,En(m)*a/V0],[En(m),En(m)]) %nivel de energía
             B=-airy(1,z0)/airy(3,z0);
             C=exp(k*a)*(airy(0,za)+B*airy(2,za));
             f=@(z) airy(0,z)+B*airy(2,z);
             g=@(z) f(z).^2;
             A=sqrt(0.5/(a*integral(g,z0,za)/rr+C^2*exp(-2*k*a)/(2*k)));
             y=A*((x<=a).*f(rr*(x/a-En(m)/V0))+C*(x>a).*exp(-k*x));
             yy=[fliplr(y)+En(m),y+En(m)];
        else %si es impar m=2,4,6...
            line([-En(m)*a/V0,En(m)*a/V0],[En(m),En(m)]) %nivel de energía
            B=-airy(0,z0)/airy(2,z0);
            C=exp(k*a)*(airy(0,za)+B*airy(2,za));
            f=@(z) airy(0,z)+B*airy(2,z);
             g=@(z) f(z).^2;
             A=sqrt(0.5/(a*integral(g,z0,za)/rr+C^2*exp(-2*k*a)/(2*k)));
             y=A*((x<=a).*f(rr*(x/a-En(m)/V0))+C*(x>a).*exp(-k*x));
             yy=[-fliplr(y)+En(m),y+En(m)];
        end 
        plot([-fliplr(x),x],yy)
        text(4,En(m)+0.1,sprintf('%1.3f',En(m)),'VerticalAlignment','bottom', 
'HorizontalAlignment','right')
    end
    hold off
    ylim([0,V0])
    grid on
    xlabel('x')
    ylabel('V(x),\Psi(x)')
    title('Niveles de energía y funciones de onda')

    function res=fImpar(E)    
        z0=-rr*E/V0;
        za=rr*(1-E/V0);
        k=sqrt(2*(V0-E));
        res=(a*k.*airy(0,za)+rr*airy(1,za)).*airy(2,z0)-airy(0,z0).
*(a*k.*airy(2,za)+rr*airy(3,za));
    end
    function res=fPar(E)    
        z0=-rr*E/V0;
        za=rr*(1-E/V0);
        k=sqrt(2*(V0-E));
        res=(a*k.*airy(0,za)+rr*airy(1,za)).*airy(3,z0)-airy(1,z0).
*(a*k.*airy(2,za)+rr*airy(3,za));
    end

    function r = raices(f, x)
        y=f(x);
        indices=find(y(1:end-1).*y(2:end)<0);
        r=zeros(1,length(indices));
        for m=1:length(indices)
            r(m)=fzero(f, [x(indices(m)), x(indices(m)+1)]);
        end
    end
end