Resonancias

Sistema simple

Resolveremos la ecuación de Schrödinger para una barrera de potencial delta de Dirac, centrada en la posición x=a, y una pared de potencial de altura infinita situda en el origen x≤0

2 2m d 2 ψ(x) d x 2 +V(x)·ψ(x)=E·ψ(x) V(x)={ ,x0 2 2m α·δ(xa),x>0

Representamos el potencial delta de Dirac por la función que se muestra en la figura, una barrera de anchura 2ε y altura que tiende a infinito, siendo ε→0

En x=a

Vamos a calcular el cuadrado del módulo de la amplitud |A|2 en función de el coeficiente D, resolviendo el sistema de dos ecuaciones

{ 2iAsin(ka)=Cexp(ika)+Dexp(ika) 2kAcos( ka )+2αAsin(ka)=kCexp( ika )kDexp( ika ) A=D k( cos(ka)icos(ka) ) kcos( ka )αsin(ka)+iksin(ka) A=Dk ( cos(ka)icos(ka) )( kcos( ka )αsin(ka)iksin(ka) ) ( kcos( ka )αsin(ka)+iksin(ka) )( kcos( ka )αsin(ka)iksin(ka) ) = A=Dk ( αsin(ka)cos( ka )+k )iα sin 2 (ka) ( kcos( ka )+αsin(ka) ) 2 + k 2 sin 2 (ka) A 2 = D 2 k 2 ( αsin(ka)cos( ka )+k ) 2 + α 2 sin 4 (ka) ( k 2 +2kαsin(ka)cos( ka )+ α 2 sin 2 (ka) ) 2 = D 2 k 2 k 2 +2kαsin(ka)cos( ka )+ α 2 sin 2 (ka) A 2 = D 2 k 2 k 2 +kαsin(2ka)+ α 2 2 ( 1cos(2ka) )

Derivamos respecto de k, para calcular los extremos (máximos y mínimos) de |A|2

d A 2 dk =0 k( 1αa )sin(2ka)( α+2 k 2 a )cos(2ka)+α=0

Establecemos un sistema de unidades en el que ℏ=m=1. Representamos la amplitud |A|2 en función de la energía E, k= 2E . Resolvemos la ecuacion transcendente en k, calculando las energías E de los máximos y minimos de la amplitud

La posición de la barrera es a=3, el parámetro α=10

function resonancias
    ee=linspace(0,5,100); %energías
    a=3; %posición barrera
    alfa=20; %parámetro
    hold on
    g=@(k) (1-alfa*a)*k.*sin(2*k*a)-(alfa+2*k.^2*a).*cos(2*k*a)+alfa;
    k=sqrt(2*ee);
    rr=raices(g,k);
    disp(rr.^2/2)
    plot(rr.^2/2,0,'ro','markersize',3,'markerfacecolor','r')
    f=@(k) k./sqrt(k.^2+alfa*k.*sin(2*k*a)+alfa^2*(1-cos(2*k*a))/2);
    plot(ee,f(k))
    hold off
    xlabel('E')
    ylabel('A')
    grid on
    title('Resonancias')
    
    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

Los extremos se señalan mediante puntos de color rojo en la figura. Los mínimos corresponden al índice par y los máximos al impar

    0.5305    1.0866    2.1229    3.2131    4.7790

Caso particular

Cuando el parámetro α es muy grande, tenemos un pozo de potencial de anchura a y altura infinita.

k( 1 α a )sin(2ka)( 1+ 2 k 2 a α )cos(2ka)+1=0 kasin(2ka)cos(2ka)+10 2kasin(ka)cos(ka)+2 sin 2 (ka)=0,{ sin(ka)=0,ka=nπ tan(ka)=ka

Los primeros corresponden a los niveles de energía de un pozo de potencial de altura infinita y de anchura a. Se advierte que en el apartado 'Pozo de potencial de profundidad infinita', la anchura del pozo de potencial es 2a

nπ a = 2m 2 E , E n = n 2 π 2 2 2m a 2 ,n=1,2,3...

function resonancia_1
    ee=linspace(0,5,100);
    a=3; %posición barrrera
    alfa=50; %parámetro
    %exacto
    g=@(k) (1-alfa*a)*k.*sin(2*k*a)-(alfa+2*k.^2*a).*cos(2*k*a)+alfa; %exacto
    k=sqrt(2*ee);
    rr=raices(g,k);
    disp(rr.^2/2)
    %aproximado (alfa grande)
    g=@(k) sin(k*a)-a*k.*cos(k*a); 
    rr=raices(g,k);
    aprox=sort([(1:length(rr))*pi/a,rr]);
    disp(aprox.^2/2)

    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

Los mínimos corresponden al índice par y los máximos al impar

    0.5411    1.1071    2.1644    3.2725    4.8700 %exacto
    0.5483    1.1217    2.1932    3.3155           %aproximado

Funciones de onda

{ ψ I (x)=2iAsin(kx),0<x<a A=Dk ( αsin(ka)cos( ka )+k )iα sin 2 (ka) ( kcos( ka )+αsin(ka) ) 2 + k 2 sin 2 (ka) { ψ II (x)=Cexp(ikx)+Dexp(ikx),x>a Cexp(ika)=2iAsin(ka)Dexp(ika)

Sea a=3, α=20, tomamos D=1. Representamos la función de onda, por ejemplo, para el segundo máximo de amplitud |A|2, E=2.1229

a=3; %posición barrera
alfa=20; %parámetro
hold on
e=2.1229; %energía (máximo)
k=sqrt(2*e);
D=1;
A=-D*k*(alfa*sin(k*a)*cos(k*a)+k-1i*alfa*sin(k*a)^2)/
((k*cos(k*a)+alfa*sin(k*a))^2+k^2*sin(k*a)^2);
 C=(2*1i*A*sin(k*a)-D*exp(-1i*k*a))*exp(-1i*k*a);
disp(C*C') %comprobación |C|^2=1
fplot(@(x) real(2*1i*A*sin(k*x)),[0,a])
fplot(@(x) real(C*exp(1i*k*x)+D*exp(-1i*k*x)),[a,3*a]) 
hold off
xlabel('x')
ylabel('\Psi(x)')
grid on
title('Función de onda')

Cambiamos la energía al segundo mínimo, E=3.2131

Comprobamos que el cuadrado de la amplitud de la función de onda de las partículas incidentes |D|2 es igual a la de las reflejadas |C|2=1

   1.0000

Desfase

La representación gráfica de la función de onda en los intervalos 0<x<a y x>a, nos sugiere que la función de onda se puede simplificar cambiando de escala

{ ψ I (x)=Asin(kx),0<x<a ψ II (x)=sin( k(xa)+δ ),x>a

La amplitud A y la fase δ se calculan a partir de las condiciones de continuidad

{ ψ I (a)= ψ II (a) 2 2m ( d ψ II (x) dx | a d ψ I (x) dx | a )+α 2 2m ψ I (a)=0 { Asin(ka)=sinδ ( kcosδkAcos(ka) )+αAsin(ka)=0

Dividiendo las dos ecuaciones, eliminamos la amplitud A y despejamos la fase δ

tanδ= sin(ka) cos(ka)+ α k sin(ka) = tan(ka) 1+ α k tan(ka)

La fase se anula, cuando sin(ka)=0, la energía es

ka=nπ,E= n 2 π 2 2 a 2 ,n=1,2,3...

Tomamos la primera ecuación, utilizando la relación entre sinδ y tanδ, despejamos la amplitud A

sinδ= tanδ 1+ tan 2 δ Asin(ka)= sin(ka) sin 2 (ka)+ ( cos(ka)+ α k sin(ka ) 2 A= k k 2 +2kαsin(ka)cos( ka )+ α 2 sin 2 (ka)

El mismo resultado que hemos obtenido anteriormente

Con los datos, a=3, y α=20. Representamos la función de onda, para tres valores de la energía:

a=3; %posición barrera
alfa=20; %parámetro
hold on
colores=['k','b','r'];
j=1;
for e=[2.1932,1.0866,2.1229]
    k=sqrt(2*e);
    A=k/sqrt(k^2+2*alfa*k*sin(k*a)*cos(k*a)+(alfa*sin(k*a))^2);
    delta=atan2(sin(k*a), cos(k*a)+alfa*sin(k*a)/k);
    fplot(@(x) A*sin(k*x),[0,a], 'color',colores(j))
    fplot(@(x) sin(k*(x-a)+delta),[a,3*a],'color',colores(j))
    j=j+1;
end
hold off
xlabel('x')
ylabel('\Psi(x)')
grid on
title('Función de onda')

Pozo de potencial asimétrico con barrera

En la página titulada Otros potenciales rectangulares de potencial, estudiamos el potencial asimétrico de profundidad V0 y anchura a, calculamos los niveles de energía y representamos las funciones de onda

En este apartado, le añadimos una barrera de potencial de altura V1 y de anchura b, tal como se muestra en la figura

El interés de este ejemplo, estriba en el estudio de la amplitud de la función de onda en la cavidad (región I) cuando las partículas de masa m y energía V0<E<V0+V1 inciden sobre la barrera

Estudiaremos las soluciones de la ecuación de Schrödinger para el potencial V(x) en los intervalos de energías E<V0 y V0<E<V0+V1

2 2m d 2 ψ(x) d x 2 +V(x)·ψ(x)=E·ψ(x) V(x)={ ,x0 0,0<xa V 0 + V 1 ,a<xa+b V 0 ,x>a+b

E<V0

Las condiciones de continuidad de la función de onda y de su derivada primera son

Niveles de energía

Eliminamos F en este sistema de dos ecuaciones

( α 3 α 2 )Cexp( α 2 (a+b) )+( α 3 + α 2 )Dexp( α 2 (a+b) )=0

Junto con las expresiones de C y D obtenemos la ecuación transcendente de la energía E

( α 3 α 2 1 )( α 2 sin( k 1 a) k 1 cos( k 1 a) )exp( α 2 b )+( α 3 α 2 +1 )( α 2 sin( k 1 a)+ k 1 cos( k 1 a) )exp( α 2 b )=0

Calculamos y representamos los niveles de energía para un pozo de potencial de anchura a=3, profundidad V0=5 y una barrera de anchura b=1 y altura V1=2

function resonancia_3
    V0=5; %profundidad pozo
    a=3; %anchura pozo
    V1=2; %altura barrera
    b=1; %anchura barrera
    
    ee=linspace(0,V0, 20);
    rr=raices(@ecuacion,ee);
    disp(rr)
    %potencial
    xx=[-0.1, 0, 0, a, a, a+b, a+b,2*a, 2*a, -0.1, -0.1];
    yy=[V0+V1+1, V0+V1+1, 0, 0, V0+V1, V0+V1 V0, V0, -0.1, -0.1, V0+1];
    fill(xx,yy, [0.5 0.5 0.5])
    %niveles de energía
    for E=rr
        line([0,a+1],[E,E],'color','b')
        text(a, E,num2str(E),'VerticalAlignment','bottom', 
'HorizontalAlignment','right')
   end
    xlabel('')
    ylabel('E')
    title('Niveles de energía')
    
    function res=ecuacion(E)
        alfa_2=sqrt(2*(V1+V0-E));
        alfa_3=sqrt(2*(V0-E));
        k1=sqrt(2*E);
        res=(alfa_3./alfa_2-1).*(alfa_2.*sin(k1*a)-k1.*cos(k1*a)).*
exp(-alfa_2*b)+(alfa_3./alfa_2+1).*(alfa_2.*sin(k1*a)+k1.*cos(k1*a)).
*exp(alfa_2*b);
    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 k=1:length(indices)
            r(k)=fzero(f, [x(indices(k)), x(indices(k)+1)]);
        end
    end
end

Los niveles de energía que hemos obtenido son muy parecidos a los de un potencial asimétrico para los mismo valores de la anchura a=3 y profundidad V0=5

Funciones de onda

{ ψ I (x)=2iAsin( k 1 x),0<xa ψ II (x)=Cexp( α 2 x)+Dexp( α 2 x),a<xa+b ψ III (x)=Fexp( α 3 x),x>a+b C=iA α 2 sin( k 1 a) k 1 cos( k 1 a) α 2 exp( α 2 a) D=iA α 2 sin( k 1 a)+ k 1 cos( k 1 a) α 2 exp( α 2 a) F=( Cexp( α 2 (a+b) )+Dexp( α 2 (a+b) ) )exp( α 3 (a+b) )

El coeficiente A o iA se calcula de modo que la integral

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

El resultado es

2 A 2 ( a sin(2 k 1 a) 2 k 1 )+2CD·b C 2 2 α 2 ( exp( 2 α 2 (a+b) )exp( 2 α 2 a ) )+ D 2 2 α 2 ( exp( 2 α 2 (a+b) )exp( 2 α 2 a ) )+ F 2 2 α 2 exp( 2 α 3 (a+b) )=1

Representamos las funciones de onda, correspondientes a los tres niveles de energía. Vemos que son similares a las representadas para un pozo de potencial asimétrico, salvo que tienden a cero más rápidamente, al estar por encima, en la región II, una barrera de potencial de altura V1.

function resonancia_4
    V0=5; %profundidad pozo
    a=3; %anchura pozo
    V1=2; %altura barrera
    b=1; %anchura barrera 
    ee=linspace(0,V0, 20);
    rr=raices(@ecuacion,ee);
    disp(rr)
    %funciones de onda
    hold on
    for E=rr
        alfa_2=sqrt(2*(V1+V0-E));
        alfa_3=sqrt(2*(V0-E));
        k1=sqrt(2*E);
        C=(alfa_2*sin(k1*a)-k1*cos(k1*a))*exp(alfa_2*a)/alfa_2;
        D=(alfa_2*sin(k1*a)+k1*cos(k1*a))*exp(-alfa_2*a)/alfa_2;
        F=(C*exp(-alfa_2*(a+b))+D*exp(alfa_2*(a+b)))*exp(alfa_3*(a+b));
       suma=2*(a-sin(2*k1*a)/(2*k1))+2*C*D*b-C^2*(exp(-2*alfa_2*(a+b))-
exp(-2*alfa_2*a))/(2*alfa_2)+D^2*(exp(2*alfa_2*(a+b))-exp(2*alfa_2*a))
/(2*alfa_2)+F^2*exp(-2*alfa_3*(a+b))/(2*alfa_3);
       A=1/sqrt(suma);
       f1=@(x) 2*A*sin(k1*x);
       f2=@(x) (C*exp(-alfa_2*x)+D*exp(alfa_2*x))*A;
       f3=@(x) A*F*exp(-alfa_3*x);
       fplot(f1,[0,a])
       fplot(f2,[a,a+b])
       fplot(f3,[a+b,2*a])
    end
    hold off
    grid on
    xlabel('x')
    ylabel('\Phi(x)')
    title('Funciones de onda')

     function res=ecuacion(E)
        alfa_2=sqrt(2*(V1+V0-E));
        alfa_3=sqrt(2*(V0-E));
        k1=sqrt(2*E);
        res=(alfa_3./alfa_2-1).*(alfa_2.*sin(k1*a)-k1.*cos(k1*a)).*
exp(-alfa_2*b)+(alfa_3./alfa_2+1).*(alfa_2.*sin(k1*a)+k1.*cos(k1*a)).
*exp(alfa_2*b);
    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 k=1:length(indices)
            r(k)=fzero(f, [x(indices(k)), x(indices(k)+1)]);
        end
    end
end

V0<E<V1

Las condiciones de continuidad de la función de onda y de su derivada primera son

Despejamos F y G en este sistema de dos ecuaciones

{ ( i k 3 + α 2 )Cexp( α 2 (a+b) )+( i k 3 α 2 )Dexp( α 2 (a+b) )=2i k 3 Fexp( i k 3 (a+b) ) ( i k 3 α 2 )Cexp( α 2 (a+b) )+( i k 3 + α 2 )Dexp( α 2 (a+b) )=2i k 3 Gexp( i k 3 (a+b) )

Sustituyendo C y D en la primera ecuación, relacionamos A y F

F= A α 2 k 3 exp( i k 3 (a+b) ){ α 2 ( α 2 sin( k 1 a)sinh( α 2 b)+ k 1 cos( k 1 a)cosh( α 2 b) ) +i k 3 ( α 2 sin( k 1 a)cosh( α 2 b)+ k 1 cos( k 1 a)sinh( α 2 b) ) }

Sustituyendo C y D en la segunda ecuación, relacionamos A y G

G= A α 2 k 3 exp( i k 3 (a+b) ){ α 2 ( α 2 sin( k 1 a)sinh( α 2 b)+ k 1 cos( k 1 a)cosh( α 2 b) )+ i k 3 ( α 2 sin( k 1 a)cosh( α 2 b)+ k 1 cos( k 1 a)sinh( α 2 b) ) }

Comprobamos que los módulos de F y G son iguales

Despejamos el módulo de la amplitud |A| tomando el de |F| como unidad

| A |= α 2 k 3 α 2 2 ( α 2 sin( k 1 a)sinh( α 2 b)+ k 1 cos( k 1 a)cosh( α 2 b) ) 2 + k 3 2 ( α 2 sin( k 1 a)cosh( α 2 b)+ k 1 cos( k 1 a)sinh( α 2 b) ) 2 | F |

Representamos el módulo de |A| en función de la energía V0<E<V0+V1, para la anchura del pozo a=3, profundidad V0=5, anchura de la barrera b=0.5, altura V1=10. Observamos que para ciertos valores de la energía E el módulo de |A| presenta máximos, similares al ejemplo estudiado en el apartado anterior

V0=5; %profundidad pozo
a=3; %anchura pozo
V1=10; %altura barrera
b=0.5; %anchura barrera
ee=linspace(V0,V0+V1, 200);
j=1;
A=zeros(1, length(ee));
for E=ee
    alfa_2=sqrt(2*(V1+V0-E));
    k3=sqrt(2*(E-V0));
    k1=sqrt(2*E);
    A(j)=alfa_2*k3/sqrt(alfa_2^2*(alfa_2*sin(k1*a)*sinh(alfa_2*b)
+k1*cos(k1*a)*cosh(alfa_2*b))^2+k3^2*(alfa_2*sin(k1*a)*cosh(alfa_2*b)
+k1*cos(k1*a)*sinh(alfa_2*b))^2);
    j=j+1;
end
plot(ee,A)
grid on
xlabel('E')
ylabel('A')
title('Resonancias')

Representamos la función de onda para un valor de la energía E cercano a un máximo y para un valor cercano a un mínimo.

V0=5; %profundidad pozo
a=3; %anchura pozo
V1=10; %altura barrera
b=0.5; %anchura barrera
%funciones de onda
hold on
for E=[7.47,9.12]
    alfa_2=sqrt(2*(V1+V0-E));
    k3=sqrt(2*(E-V0));
    k1=sqrt(2*E);
    F=1;
    A=alfa_2*k3*exp(-1i*k3*(a+b))/(-alfa_2*(alfa_2*sin(k1*a)*
sinh(alfa_2*b)+k1*cos(k1*a)*cosh(alfa_2*b))+1i*k3*(alfa_2*sin(k1*a)*
cosh(alfa_2*b)+k1*cos(k1*a)*sinh(alfa_2*b)));
    C=(alfa_2*sin(k1*a)-k1*cos(k1*a))*exp(alfa_2*a)*1i*A/alfa_2;
    D=(alfa_2*sin(k1*a)+k1*cos(k1*a))*exp(-alfa_2*a)*1i*A/alfa_2;
    G=(C*exp(-alfa_2*(a+b))+D*exp(alfa_2*(a+b))-F*exp(-1i*k3*(a+b)))*
exp(-1i*k3*(a+b));
    f1=@(x) real(2*A*1i*sin(k1*x));
    f2=@(x) real(C*exp(-alfa_2*x)+D*exp(alfa_2*x));
    f3=@(x) real(F*exp(-1i*k3*x)+G*exp(1i*k3*x));
    fplot(f1,[0,a])
    fplot(f2,[a,a+b])
    fplot(f3,[a+b,2*a])
end
hold off
grid on
xlabel('x')
ylabel('\Psi(x)')
title('Funciones de onda')

Referencias

Herbert Massmann. Illustration of resonances and the law of exponential decay in a simple quantum-mechanical problem. Am. J. Phys. 53 (7) July 1985, pp. 679-683

One Dimensional Scattering.