El oscilador armónico cuántico

La ecuación de Schrödinger unidimensional e independiente del tiempo es

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

La energía potencial de un oscilador armónico es Ep=kx2/2, donde k es la constante elástica y m la masa de la partícula.

Tomando una escala de energías y distancias de la forma

x= ( mω ) ξE=ωεω= k m

La ecuación de Schrödinger se transforma en otra más simple

1 2 ( d 2 d ξ 2 + ξ 2 )Ψ(ξ)=εΨ(ξ)

Los niveles de energía vienen dados por

ε n =n+ 1 2

Y las funciones de onda

Ψ n (ξ)= H n (ξ) 2 n π n! exp( ξ 2 2 )

Siendo H(ξ) los polinomios de Hermite.

Polinomios de Hermite

Los primeros seis polinomios de Hermite son los siguientes.

H0(x)=1

H1(x)=2x

H2(x)=4x2-2

H3(x)=8x3-12x

H4(x)=16x4-48x2+12

H5(x)=32x5-160x3+120x

Todos los polinomios de Hermite de orden n>2 se pueden expresar en términos de los dos primeros polinomios H0(x) y H1(x), de orden cero y uno respectivamente, mediante la siguiente relación de recurrencia.

H n (x)=2x H n1 (x)2(n1) H n2 (x)

El código de la función recursiva que calcula los valores del polinomio Hn(x) es

function res=hermite(n, x)
    if n==0
        res=1;
    elseif n==1
        res=2*x;
    else
        res=2*x.*hermite(n-1,x)-2*(n-1)*hermite(n-2,x);
    end
end

Mejoramos la función calculando los coeficientes del polinomio de Hermite de grado n y mediante la función polyval de MATLAB, se calcula el valor del polinomio para un valor de x o para un vector de datos x

function res=hermite_1(n, x)
    p=zeros(n+1,1);
    p1=zeros(n,1);
    p0=zeros(n-1,1);
%valores iniciales   
    p0=[1];
    p1=[2 0];
   if n==0
       p=p0;
   elseif n==1
        p=p1;
   else
        for k=2:n
            p=2*[p1 0]-2*(k-1)*[0 0 p0];
            p0=p1;
            p1=p;
        end
   end
   res=polyval(p,x);
end

Para representar las funciones de onda del oscilador armónico cuántico

Ψ n ( x ) = H n ( x ) 2 n π n ! exp ( x 2 2 )

nos hace falta el factorial de n.

>> prod(1:5)
>> 120

que es el factorial de 5!=1·2·3·4·5=120

Nuestra tarea va consistir en representar el potencial Ep(x)=x2/2 y la función de onda para un nivel dado de energía n=0, 1, 2...

n=3; %estado
hold on
x=-4:0.1:4;
y=x.^2/2;
xx=[-4 x 4];
yy=[0 y 0];
fill(xx,yy,[0.5 0.5 0.5])
plot(x,y,'b')
line([-4 4],[n+0.5 n+0.5], 'color','k')
y=n+0.5+hermite_1(n,x).*exp(-x.^2/2)/sqrt(2^n*sqrt(pi)*prod(1:n));
plot(x,y,'r') 
hold off
xlabel('x')
ylabel('\Phi(x)')
title('Oscilador armónico cuántico')