Sistema de potenciales delta de Dirac con un defecto
En la página titulada Sistema de potenciales delta de Dirac calculamos los niveles de energía y representamos las funciones de onda, para potenciales delta de Dirac negativos. El problema se divide en dos partes: E<0 y E>0. En esta página, los potenciales delta son positivos, por lo que solamente tenemos que estudiar el caso E>0
Sustitución
Resolveremos la ecuación de Schrödinger, en un sistema formado por cuatro potenciales delta, igualmente espaciados a, el parámetro del tercero es β y de los restantes α. A partir de este ejemplo, el lector puede generalizar para cualquier número de potenciales iguales e igualmente espaciados con un único defecto, tal como se puede apreciar en el código MATLAB
El potencial V(x) es
Siendo ε→0, una cantidad que tiende a cero
Las soluciones de la ecuación de Schrödinger en las regiones I, I, III, IV y V con potencial V(x)=0 son
-
Condiciones en los extremos x=0, y x=5a donde V(x)=∞
-
En en x=a, x=2a, x=3a y x=4a, la función de onda es continua
La derivada primera de la función de onda NO es continua en x=a, x=2a, x=3a y x=4a
Integramos la ecuación de Schrödinger en el pequeño intervalo de a-ε a a+ε
De modo similar, establecemos las otras tres ecuaciones que relacionan los coeficientes en x=2a (parámetro α), en x=3a (parámetro β) y x=4a (parámetro α)
En forma matricial, la relación entre los coeficientes es
Niveles de energía
Relacionamos A1 y B1 con A5 y B5
Teniendo en cuenta, la relaciones en los extremos, en x=0 (B1=0) y en x=5a, (A5sin(ka)+B5cos(ka)=0), obtenemos la ecuación transcendente
Representamos los niveles de energía para un sistama de N=4 potenciales delta de Dirac, fijamos la separación constante a=1, el parámetro α=2, de todos los potenciales excepto el tercero, n=3, cuyo parámetro β va cambiando entre 0 y 10
function delta_dirac_30 a=1; alfa=2; N=4; n=3; %posición defecto hold on for beta=0:0.5:10 ff=@(x) niveles(x); xb=buscar_intervalos(ff,0.1,4,50); nb=size(xb); for m=1:nb(1) raiz=fzero(ff,[xb(m,1) xb(m,2)]); plot(beta,raiz^2,'ro','markersize',3,'markerfacecolor','r') end end line([alfa,alfa],[0,11],'lineStyle','--') ylim([0,11]) hold off grid on xlabel('\beta'); ylabel('E') title('Defecto en el sólido lineal') function res=niveles(k) T=eye(2); Ma=[0,1; 1,-alfa/k]; Mb=[0,1; 1,-beta/k]; MI=[sin(k*a),cos(k*a);cos(k*a),-sin(k*a)]; for s=N:-1:n+1 T=(MI*Ma)*T; end T=(MI*Mb)*T; for s=n-1:-1:1 T=(MI*Ma)*T; end res=T(2,1)*cos(k*a)-T(2,2)*sin(k*a); end function xb = buscar_intervalos(f,a,b,n) x = a:(b-a)/n:b; j = 0; y1=f(x(1)); xb=[]; for i = 1:length(x)-1 y2=f(x(i+1)); if sign(y1) ~= sign(y2) j = j + 1; xb(j,1) = x(i); xb(j,2) = x(i+1); end y1=y2; end if isempty(xb) disp('no se han encontrado cambios de signo') else disp(['número de intervalos:' int2str(j)]) end end end
Cuando el parámetro β=α=2, línea a trazos, no hay defecto.
Funciones de onda
Dados los coeficientes A1 y B1=0, obtenemos los coeficientes A2 y B2, A3 y B3, A4 y B4, A5 y B5. En forma matricial
El coeficiente A1 se determina de modo que la suma
Representamos las funciones de onda correspondientes a los dos primeros niveles de energía para un sistama de N=4 potenciales delta de Dirac, fijamos la separación constante a=1, el parámetro α=2, de todos los potenciales excepto el tercero, n=3, cuyo parámetro β=0.5
function delta_dirac_27 a=1; alfa=2; beta=0.5; N=4; n=3; %posición defecto ff=@(x) niveles(x); xb=buscar_intervalos(ff,0.1,4,50); nb=size(xb); kn=zeros(1,nb(1)); for m=1:nb(1) kn(m)=fzero(ff,[xb(m,1) xb(m,2)]); end disp(kn) A=zeros(1,N+1); B=zeros(1,N+1); A(1)=1; B(1)=0; for k=kn(1:2) NaI=[alfa/k,1;1,0]; NbI=[beta/k,1;1,0]; M=[sin(k*a),cos(k*a);cos(k*a),-sin(k*a)]; suma=A(1)^2*(a-sin(2*k*a)/(2*k))/2+B(1)^2*(a+sin(2*k*a)/(2*k))/2- A(1)*B(1)*(cos(2*k*a)-1)/(2*k); for m=1:n-1 Z=(NaI*M)*[A(m);B(m)]; A(m+1)=Z(1); B(m+1)=Z(2); suma=suma+A(m+1)^2*(a-sin(2*k*a)/(2*k))/2+B(m+1)^2*(a+sin(2*k*a)/(2*k)) /2-A(m+1)*B(m+1)*(cos(2*k*a)-1)/(2*k); end m=n; Z=(NbI*M)*[A(m);B(m)]; A(m+1)=Z(1); B(m+1)=Z(2); suma=suma+A(m+1)^2*(a-sin(2*k*a)/(2*k))/2+B(m+1)^2*(a+sin(2*k*a)/(2*k)) /2-A(m+1)*B(m+1)*(cos(2*k*a)-1)/(2*k); for m=n+1:N Z=(NaI*M)*[A(m);B(m)]; A(m+1)=Z(1); B(m+1)=Z(2); suma=suma+A(m+1)^2*(a-sin(2*k*a)/(2*k))/2+B(m+1)^2*(a+sin(2*k*a)/(2*k)) /2-A(m+1)*B(m+1)*(cos(2*k*a)-1)/(2*k); end hold on for m=0:N fplot(@(x) (A(m+1)*sin(k*(x-m*a))+B(m+1)*cos(k*(x-m*a)))/sqrt(suma), [m*a,(m+1)*a]) end end hold off grid on xlabel('x') ylabel('\Phi(x)') title('Funciones de onda') function res=niveles(k) T=eye(2); Ma=[0,1; 1,-alfa/k]; Mb=[0,1; 1,-beta/k]; MI=[sin(k*a),cos(k*a);cos(k*a),-sin(k*a)]; for s=N:-1:n+1 T=(MI*Ma)*T; end T=(MI*Mb)*T; for s=n-1:-1:1 T=(MI*Ma)*T; end res=T(2,1)*cos(k*a)-T(2,2)*sin(k*a); end function xb = buscar_intervalos(f,a,b,n) x = a:(b-a)/n:b; j = 0; y1=f(x(1)); xb=[]; for i = 1:length(x)-1 y2=f(x(i+1)); if sign(y1) ~= sign(y2) j = j + 1; xb(j,1) = x(i); xb(j,2) = x(i+1); end y1=y2; end if isempty(xb) disp('no se han encontrado cambios de signo') else disp(['número de intervalos:' int2str(j)]) end end end
Las primeras raíces k de la ecuación transcendente son
1.0145 1.8762 3.1416 3.5450
Desplazamiento
Resolveremos la ecuación de Schrödinger, en un sistema formado por cuatro potenciales delta, igualmente espaciados a, el parámetro de todos los potenciales es α. El tercer potencial delta se desplaza δ<a. A partir de este ejemplo, el lector puede generalizar para cualquier número de potenciales iguales con un único defecto, tal como se puede apreciar en el código MATLAB
El potencial V(x) es
Siendo ε→0, una cantidad que tiende a cero
Las soluciones de la ecuación de Schrödinger en las regiones I, I, III, IV y V con potencial V(x)=0 son
-
Condiciones en los extremos x=0, y x=5a donde V(x)=∞
-
En en x=a, x=2a, x=3a+δ y x=4a, la función de onda es continua
La derivada primera de la función de onda NO es continua en x=a, x=2a, x=3a+δ y x=4a. En la sección anterior hemos visto cómo se relacionan los coeficientes. En ésta, examinamos las posiciones x=3a+δ y x=4a
Integramos la ecuación de Schrödinger en el pequeño intervalo de 3a+δ-ε a 3a+δ+ε
Integramos la ecuación de Schrödinger en el pequeño intervalo de 4a-ε a 4a+ε
Las cuatro ecuaciones que relacionan los coeficientes en x=a, x=2a, x=3a+δ y x=4a son
En forma matricial, la relación entre los coeficientes es
Niveles de energía
Relacionamos A1 y B1 con A5 y B5
Representamos los niveles de energía para un sistama de N=4 potenciales delta de Dirac, fijamos la separación constante a=1, el parámetro α=2, de todos los potenciales. El tercero, n=3, se desplaza δ<a que va cambiando entre 0 y 0.9
function delta_dirac_29 a=1; alfa=2; N=4; n=3; %posición defecto hold on for delta=linspace(0,0.9,20) ff=@(x) niveles(x); xb=buscar_intervalos(ff,0.1,4,50); nb=size(xb); for m=1:nb(1) raiz=fzero(ff,[xb(m,1) xb(m,2)]); plot(delta,raiz^2,'ro','markersize',3,'markerfacecolor','r') end end hold off grid on ylim([0,11]) xlabel('\delta'); ylabel('E') title('Defecto en el sólido lineal') function res=niveles(k) T=eye(2); M=[0,1; 1,-alfa/k]; MI=[sin(k*a),cos(k*a);cos(k*a),-sin(k*a)]; MIa=[sin(k*(a+delta)),cos(k*(a+delta));cos(k*(a+delta)),-sin(k*(a+delta))]; MIb=[sin(k*(a-delta)),cos(k*(a-delta));cos(k*(a-delta)),-sin(k*(a-delta))]; for s=N:-1:n+2 T=(MI*M)*T; end T=(MIb*M)*T; T=(MIa*M)*T; for s=n-1:-1:1 T=(MI*M)*T; end res=T(2,1)*cos(k*a)-T(2,2)*sin(k*a); end function xb = buscar_intervalos(f,a,b,n) x = a:(b-a)/n:b; j = 0; y1=f(x(1)); xb=[]; for i = 1:length(x)-1 y2=f(x(i+1)); if sign(y1) ~= sign(y2) j = j + 1; xb(j,1) = x(i); xb(j,2) = x(i+1); end y1=y2; end if isempty(xb) disp('no se han encontrado cambios de signo') else disp(['número de intervalos:' int2str(j)]) end end end
Para δ=0, no hay defecto
Funciones de onda
Dados los coeficientes A1 y B1=0, obtenemos los coeficientes A2 y B2, A3 y B3, A4 y B4, A5 y B5. En forma matricial
El coeficiente A1 se determina de modo que la suma
El resultado es
Representamos las funciones de onda correspondientes a los dos primeros niveles de energía para un sistama de N=4 potenciales delta de Dirac, fijamos la separación constante a=1, el parámetro α=2, de todos los potenciales. El tercero, n=3, se desplaza δ=0.75
function delta_dirac_28 a=1; alfa=2; delta=0.75; N=4; n=3; %posición defecto ff=@(x) niveles(x); xb=buscar_intervalos(ff,0.1,4,50); nb=size(xb); kn=zeros(1,nb(1)); for m=1:nb(1) kn(m)=fzero(ff,[xb(m,1) xb(m,2)]); end disp(kn) A=zeros(1,N+1); B=zeros(1,N+1); A(1)=1; B(1)=0; for k=kn(1:2) NI=[alfa/k,1;1,0]; M=[sin(k*a),cos(k*a);cos(k*a),-sin(k*a)]; Ma=[sin(k*(a+delta)),cos(k*(a+delta));cos(k*(a+delta)), -sin(k*(a+delta))]; Mb=[sin(k*(a-delta)),cos(k*(a-delta));cos(k*(a-delta)), -sin(k*(a-delta))]; suma=A(1)^2*(a-sin(2*k*a)/(2*k))/2+B(1)^2*(a+sin(2*k*a)/(2*k))/2 -A(1)*B(1)*(cos(2*k*a)-1)/(2*k); for m=1:n-1 Z=(NI*M)*[A(m);B(m)]; A(m+1)=Z(1); B(m+1)=Z(2); suma=suma+A(m+1)^2*(a-sin(2*k*a)/(2*k))/2+B(m+1)^2*(a+sin(2*k*a) /(2*k))/2-A(m+1)*B(m+1)*(cos(2*k*a)-1)/(2*k); end suma=suma+A(m+1)^2*(a+delta-sin(2*k*(a+delta))/(2*k))/2+B(m+1)^2* (a+delta+sin(2*k*(a+delta))/(2*k))/2-A(m+1)*B(m+1)*(cos(2*k*(a+delta))-1) /(2*k); m=n; Z=(NI*Ma)*[A(m);B(m)]; A(m+1)=Z(1); B(m+1)=Z(2); suma=suma+A(m+1)^2*(a-delta-sin(2*k*(a-delta))/(2*k))/2+B(m+1)^2* (a-delta+sin(2*k*(a-delta))/(2*k))/2-A(m+1)*B(m+1)*(cos(2*k*(a-delta))-1) /(2*k); m=n+1; Z=(NI*Mb)*[A(m);B(m)]; A(m+1)=Z(1); B(m+1)=Z(2); suma=suma+A(m+1)^2*(a-sin(2*k*a)/(2*k))/2+B(m+1)^2*(a+sin(2*k*a) /(2*k))/2-A(m+1)*B(m+1)*(cos(2*k*a)-1)/(2*k); for m=n+2:N Z=(NI*M)*[A(m);B(m)]; A(m+1)=Z(1); B(m+1)=Z(2); suma=suma+A(m+1)^2*(a-sin(2*k*a)/(2*k))/2+B(m+1)^2*(a+sin(2*k*a) /(2*k))/2-A(m+1)*B(m+1)*(cos(2*k*a)-1)/(2*k); end hold on for m=0:n-2 fplot(@(x) (A(m+1)*sin(k*(x-m*a))+B(m+1)*cos(k*(x-m*a)))/ sqrt(suma),[m*a,(m+1)*a]) end m=n-1; fplot(@(x) (A(m+1)*sin(k*(x-m*a))+B(m+1)*cos(k*(x-m*a)))/ sqrt(suma),[m*a,(m+1)*a+delta]) m=n; fplot(@(x) (A(m+1)*sin(k*(x-m*a-delta))+B(m+1)*cos(k*(x-m*a-delta))) /sqrt(suma),[m*a+delta,(m+1)*a]) for m=n+1:N fplot(@(x) (A(m+1)*sin(k*(x-m*a))+B(m+1)*cos(k*(x-m*a))) /sqrt(suma),[m*a,(m+1)*a]) end end hold off grid on xlabel('x') ylabel('\Phi(x)') title('Funciones de onda') function res=niveles(k) T=eye(2); M=[0,1; 1,-alfa/k]; MI=[sin(k*a),cos(k*a);cos(k*a),-sin(k*a)]; MIa=[sin(k*(a+delta)),cos(k*(a+delta));cos(k*(a+delta)),-sin(k*(a+delta))]; MIb=[sin(k*(a-delta)),cos(k*(a-delta));cos(k*(a-delta)),-sin(k*(a-delta))]; for s=N:-1:n+2 T=(MI*M)*T; end T=(MIb*M)*T; T=(MIa*M)*T; for s=n-1:-1:1 T=(MI*M)*T; end res=T(2,1)*cos(k*a)-T(2,2)*sin(k*a); end function xb = buscar_intervalos(f,a,b,n) x = a:(b-a)/n:b; j = 0; y1=f(x(1)); xb=[]; for i = 1:length(x)-1 y2=f(x(i+1)); if sign(y1) ~= sign(y2) j = j + 1; xb(j,1) = x(i); xb(j,2) = x(i+1); end y1=y2; end if isempty(xb) disp('no se han encontrado cambios de signo') else disp(['número de intervalos:' int2str(j)]) end end end
Las primeras raíces k de la ecuación transcendente son
1.2432 1.8356 2.3518 2.7170 3.1965
Referencias
Todd K. Timberlake, Neilson Woodfield. Band formation and defects in a finite periodic quantum potential. Am. J. Phys. 90 (2), February 2022. pp. 93-102