La ecuación de Schrödinger en coordenadas esféricas (II)

La ecuación de Schrödinger en la dirección radial es

1 r 2 d dr ( r 2 dR dr )+{ 2m 2 ( EV(r) ) l(l+1) r 2 }R=0

donde V(r) es el potencial

V(r)={ V 0 ,0ra 0,r>a

Tenemos que resolver dos ecuaciones diferenciales, una para el interior de la esfera de radio a y otra, para el exterior. En r=a la componente radial de la función de onda R(r) deberá ser continua y también su derivada primera. Lo que nos permitirá calcular los niveles de energía En dado el valor del número l

{ 1 r 2 d dr ( r 2 dR dr )+{ 2m 2 ( E+ V 0 ) l(l+1) r 2 }R=0,0ra 1 r 2 d dr ( r 2 dR dr )+{ 2m 2 E l(l+1) r 2 }R=0,r>a

Haciendo el cambio de variable, u=rR(r),

R= u r dR dr =( r du dr u ) 1 r 2 , d dr ( r 2 dR dr )=r d 2 u d r 2

Las ecuaciones diferenciales se transforman en

{ d 2 u d r 2 +{ q 2 l(l+1) r 2 }u=0, q 2 = 2m 2 ( E+ V 0 ),0ra d 2 u d r 2 +{ k 2 l(l+1) r 2 }u=0, k 2 = 2m 2 E,r>a

k2>0 y q2>0 ya que los niveles de energía En se encuentran en el intervalo -V0<En<0

Caso particular: l=0

El caso más sencillo se presenta cuando l=0

{ d 2 u d r 2 + q 2 u=0, q 2 = 2m 2 ( E+ V 0 ),0ra d 2 u d r 2 k 2 u=0, k 2 = 2m 2 E,r>a q 2 = 2m 2 ( E+ V 0 )= k 2 + 2m 2 V 0

Niveles de energía

En r=a, la solución u(r) o R(r) ha de ser continua y también su derivada primera

{ Asin( qa )=Cexp( ka ) Aqcos( qa )=Ckexp( ka ) tan( qa )= q k

Como hemos visto, los coeficientes q y k están relacionados

( qa ) 2 + ( ka ) 2 = r 0 2 , r 0 2 = 2m 2 V 0 a 2

Estas dos últimas ecuaciones determinan los posibles valores de k o los niveles En<0 de la energía. Llamamos x=qa, y=ka

{ y= x tan(x) x 2 + y 2 = r 0 2

Representamos circunferencias de radio r0=π/2, 3π/2, 5π/2, 7π/2 y la función y=-x/tan(x)

hold on
for r0=pi/2:pi:7*pi/2
    fplot(@(x) r0*cos(x), @(x) r0*sin(x),[0,pi/2])
end
for n=0:3
    fplot(@(x) -x./tan(x),[pi/2+n*pi, pi+n*pi])
end
hold off
ylim([0,12])
xlabel('x')
ylabel('y')
grid on
set(gca,'XTick',0:pi/2:4*pi)
set(gca,'XTickLabel',{'0','\pi/2','\pi','3\pi/2','2\pi','5\pi/2','3\pi',
'7\pi/2','4\pi'})
title('Niveles de energía')

Dado el valor de r0 o de la energía potencial -V0. Las ordenadas y=ka de los puntos de intersección, son los posibles valores de k y por tanto, de la energía E

La primera ecuación se escribe, alternativamente, ysin(x)+xcos(x)=0

Representamos la función

f(y)=ysin( r 0 2 y 2 )+ r 0 2 y 2 cos( r 0 2 y 2 )

para r0=π/2, 3π/2, 5π/2, 7π/2. Calculamos las raíces múltiples de la ecuación transcendente f(y)=0, mediante la función raices. Las señalamos mediante puntos de color rojo

function pozo_potencial_1
    hold on
    for r0=pi/2:pi:7*pi/2
        f=@(y) y.*sin(sqrt(r0^2-y.^2))+sqrt(r0^2-y.^2).*cos(sqrt(r0^2-y.^2));
        xx=linspace(0,r0,50);
        rr=raices(f,xx);
        for n=1:length(rr)
            plot(rr(n),0,'ro','markersize',4,'markerfacecolor','r')
        end
        disp(rr)
        fplot(f,[0,r0])
    end
    hold off
    xlabel('y')
    ylabel('f(y)')
    grid on
    title('Raíces')

    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

r0/y=ka
π/2
3π/23.9526
5π/25.60057.3456
7π/26.93119.381810.6126

Función de onda radial R(r)

Una vez que tenemos los niveles de energía En o los posibles valores de k para un valor dado de V0 o r0, determinamos la parte radial de las funciones de onda

R(r)={ A sin( qr ) r ,0ra C exp( kr ) r ,r>a Asin( qa )=Cexp( ka ) ( qa ) 2 + ( ka ) 2 = r 0 2

Las dos últimas ecuaciones, relacionan los coeficientes A y C, q y k

El coeficiente A se determina de modo que la integral

0 a A 2 sin 2 ( qr ) r 2 r 2 dr + a C 2 exp( 2kr ) r 2 r 2 dr=1

El resultado es

A 2 { 0 a sin 2 ( qr )dr+ sin 2 ( qa ) exp( 2ka ) a exp( 2kr )dr }=1 A 2 { 1 2 ( a sin( 2qa ) 2q )+ sin 2 ( qa ) 2k }=1

Representamos las funciones de onda correspondientes a los tres niveles de energía para r0=7π/2. El nivel más bajo de energía (estado fundamental) es el valor más grande de k=10.6126

function pozo_potencial_2
    a=1;
    r0=7*pi/2;
    f=@(y) y.*sin(sqrt(r0^2-y.^2))+sqrt(r0^2-y.^2).*cos(sqrt(r0^2-y.^2));
    xx=linspace(0,r0,50);
    rr=raices(f,xx);
    disp(rr)
    hold on
    for k=rr
        q=sqrt(r0^2-k^2);
        A=1/sqrt((a/2-sin(2*q*a)/(4*q)+sin(q*a)^2/(2*k)));
        C=A*sin(q*a)/exp(-k*a);
        R=@(r) (A*(r<=a).*sin(q*r))./r+(C*(r>a).*exp(-k*r))./r;
        fplot(R,[0,2*a])
    end
    hold off
    xlabel('r')
    legend('2','1','fundamental','Location', 'best')
    ylabel('R(r)')
    grid on
    title('Función de onda, l=0')

    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 posibles valores de ka son

6.9311    9.3818   10.6126

Solución general, l≠0

Niveles de energía

La componente radial R(r) de la función de onda ha de ser continua y también su derivada primera en r=a

{ A j l ( qa )=C h l (1) ( ika ) A d j l ( qr ) dr | r=a =C d h l (1) ( ikr ) dr | r=a

Teniendo en cuenta las propiedades de estas funciones

d j l (z) dz = j l1 (z) l+1 z j l (z) d h l (1) (z) dz = h l1 (1) (z) l+1 z h l (1) (z)

El sistema de dos ecuaciones se convierte en una ecuación transcendente

{ A j l ( qa )=C h l (1) ( ika ) A k 0 { j l1 ( qa ) l+1 k 0 a j l ( qa ) }=Cik{ h l1 (1) ( ika ) l+1 ika h l (1) ( ika ) } k 0 ( j l1 ( qa ) l+1 k 0 a j l ( qa ) ) j l ( qa ) = ik( h l1 (1) ( ika ) l+1 ika h l (1) ( ika ) ) h l (1) ( ika ) q· j l1 ( qa ) h l (1) ( ika )ik h l1 (1) ( ika ) j l ( qa )=0

Los coeficientes q y k están relacionados

( qa ) 2 + ( ka ) 2 = r 0 2

Calculamos los ceros de la ecuación transcendente en k

f(k)= r 0 2 / a 2 k 2 · j l1 ( r 0 2 k 2 a 2 ) h l (1) ( ika )ik h l1 (1) ( ika ) j l ( r 0 2 k 2 a 2 )

La función f(k) es real para l par e imaginaria para l impar. Para calcular los ceros la función f(k) tiene que ser real por lo que se multiplica por il. El signo de la función carece de importancia

En la página titulada Funciones de Bessel el apartado, Funciones esféricas de Bessel), proporcionamos las expresiones de las primeras funciones esféricas de Bessel jl(x) y de Hankel. Vamos a trabajar con estas funciones antes de utilizar las propias de MATLAB

Establecemos el valor de r0=7π/2. Comprobamos que para l=0, obtenemos las mismos valores de k

function pozo_potencial_6
   j_1=@(x) cos(x)./x;
   j0=@(x) sin(x)./x;
   j1=@(x) (sin(x)-x.*cos(x))./(x.^2);
   j2=@(x) ((3-x.^2).*sin(x)-3*x.*cos(x))./x.^3;
   j3=@(x) (3*(5-2*x.^2).*sin(x)-x.*(15-x.^2).*cos(x))./x.^4;
   Jfunc={j_1,j0,j1,j2,j3};
   
   h_1=@(x) exp(1i*x)./x;
   h0=@(x) -1i*exp(1i*x)./x;
   h1=@(x) -((x+1i).*exp(1i*x))./(x.^2);
   h2=@(x) (1i*(x.^2+3*1i*x-3).*exp(1i*x))./x.^3;    
   h3=@(x) ((x.^3+6*1i*x.^2-15*x-15*1i).*exp(1i*x))./x.^4;
   Hfunc={h_1,h0,h1,h2,h3};
      
    l=0; %cambiar a 1,2, 3 como máximo
    a=1;
    r0=7*pi/2;
 
    q=@(y) sqrt(r0^2-y.^2);
    f=@(y) ((q(y).*Jfunc{l+1}(q(y)*a)).*Hfunc{l+2}(1i*y*a)-
1i*(y.*Hfunc{l+1}(1i*y*a)).*Jfunc{l+2}(q(y)*a))*(1i)^l;
    xx=linspace(0,r0,50);
    rr=raices(f,xx);
    disp(rr)

    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
6.9311    9.3818   10.6126

Calculamos los niveles de energía En o los posibles valores de k dado r0 o V0. Para los valores l=0, 1, 2, 3

l/k
06.93119.381810.6126
14.96868.451810.1981
21.13657.26289.6515
35.67338.9555

Los mismos resultados se obtienen utilizando las funciones besselj, besselh de MATLAB. Hay que tener en cuenta que las funciones esféricas se definen

j l (x)= π 2x J l+ 1 2 (x), h l (1) (x)= π 2x H l+ 1 2 (1) (x)

function pozo_potencial_4
    l=0; %cambiar a 1,2,3...
    a=1;
    r0=7*pi/2;
    %funciones esféricas de Bessel
    q=@(y) sqrt(r0^2-y.^2);
    h0=@(y) sqrt(pi./(2*y)).*besselh(l+1/2,1,y);
    h1=@(y) sqrt(pi./(2*y)).*besselh(l-1/2,1,y);
    j0=@(y) sqrt(pi./(2*y)).*besselj(l+1/2,y);
    j1=@(y) sqrt(pi./(2*y)).*besselj(l-1/2,y);
    
    f=@(y) ((q(y).*j1(q(y))).*h0(1i*y)-(1i*y.*h1(1i*y)).*j0(q(y)))*(1i)^l;
    xx=linspace(0,r0,50);
    rr=raices(f,xx);
    disp(rr)

    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

Función de onda radial R(r)

Una vez que tenemos los niveles de energía En o los posibles valores de k para un valor dado de V0 o r0. Determinamos la parte radial de la función de onda

R(r)={ A j l ( k 0 r ),0ra C h l (1) ( ikr ),r>a A j l ( qa )=C h l (1) ( ika ) ( qa ) 2 + ( ka ) 2 = r 0 2

Las dos últimas ecuaciones, relacionan los coeficientes A y C, q y k

El coeficiente A se determina de modo que la integral

0 a ( A j l ( qr ) ) 2 r 2 dr + a ( C h l (1) ( ikr ) ) 2 r 2 dr=1 A 2 { 0 a ( r· j l ( qr ) ) 2 dr+ ( j l ( qa ) h l (1) ( ika ) ) 2 a ( r· h l (1) ( ikr ) ) 2 dr }=1

Representamos las funciones de onda correspondientes a los tres niveles de energía para r0=7π/2 y para l=0. Comprobamos que obtenemos la misma representación gráfica que en el apartado anterior. Cambiamos a l=1.

function pozo_potencial_5
    l=1;
    a=1;
    r0=7*pi/2;
    %funciones esféricas de Bessel
    q=@(y) sqrt(r0^2-y.^2);
    h0=@(y) sqrt(pi./(2*y)).*besselh(l+1/2,1,y);
    h1=@(y) sqrt(pi./(2*y)).*besselh(l-1/2,1,y);
    j0=@(y) sqrt(pi./(2*y)).*besselj(l+1/2,y);
    j1=@(y) sqrt(pi./(2*y)).*besselj(l-1/2,y);
    f=@(y) ((q(y).*j1(q(y))).*h0(1i*y)-(1i*y.*h1(1i*y)).*j0(q(y)))*(1i)^l;
    xx=linspace(0,r0,50);
    rr=raices(f,xx);
    disp(rr)
    hold on
    for k=rr
        q=sqrt(r0^2-k^2);
        C=j0(q*a)/h0(1i*k*a);
        f_1=@(x) (x.*j0(q*x)).^2;
        f_2=@(x) (x.*h0(1i*k*x)).^2;
        A=1/sqrt(integral(f_1,0,a)+C^2*integral(f_2,a,10*a));
        C=C*A;
        R=@(r) A*(r<=a).*j0(q*r)+C*(r>a).*h0(1i*k*r);
        fplot(R,[0,2*a])
    end
    hold off
    xlabel('r')
    ylabel('R(r)')
    legend('2','1','fundamental','Location', 'best')
    grid on
    title('Función de onda, l=1')

    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

El nivel más bajo de energía (estado fundamental) es el valor más grande de k=10.1981

4.9686    8.4518   10.1981

Referencias

Masatsugu Sei Suzuki, Itsuko S. Suzuki. Finite spherical square well potential: deuteron, with the use of Mathematica https://www.researchgate.net/publication/348819348_Finite_spherical_square_ well_potential_deuteron_with_the_use_of_Mathematica