Energía potencial efectiva

La energía de una partícula en un campo de fuerzas centrales descrito por un potencial V(r) se expresan en coordenadas polares

E= 1 2 m v 2 +V(r) E= 1 2 m ( dr dt ) 2 + 1 2 m r 2 ( dθ dt ) 2 +V(r)

El momento angular se expresa

L=m r 2 dθ dt

Combinando ambas expresiones

E= 1 2 m ( dr dt ) 2 + L 2 2m r 2 +V(r)

Se denomina potencial efectivo a la suma de los dos últimos términos que dependen de la distancia radial r

V ef (r)= L 2 2m r 2 +V(r)

Cuando V(r)=-GMm/r encontramos que

Potencial de Yukawa

El potencial de Yukawa depende de la distancia radial r de la forma

V(r)= k r exp( r λ )

El potencial efectivo es

V ef (r)= k r exp( r λ )+ L 2 2m r 2 ,r0

Definimos la distancia adimensional, x=r/λ

V ef (x)= k λ ( e x x + 1 2α x 2 ),α= mk L 2 λ

Estudiamos la función dependiente de la variable x definida entre paréntesis

y= e x x + 1 2α x 2

Parámetro α=0.9

Representamos la función y(x) para el valor del parámetro α=0.9

a=0.9; %parámetro alfa
f=@(x) -exp(-x)./x+1./(2*a*x.^2);
fplot(f,[0.5,5])
E=0.2; %energía
g=@(x) f(x)-E;
xm=fzero(g,[0.5,2]); %distancia mínima
line([xm,5],[E,E],'color','k')
line([xm,xm],[0,E],'lineStyle','--','color','k')
ylim([0,0.5])
grid on
xlabel('x')
ylabel('y')
title('Potencial de Yukawa')

La función Vef(r) es similar a la de un campo de fuerzas repulsivo inversamente proporcional al cuadrado de la distancia

Cuando la energía de la partícula es E=0.2, la intesección de la recta horizontal y la función energía potencial efectiva, define la distancia de máximo acercamiento xm al centro de fuerzas

La ecuación E=y(x) es transcendente por lo que utilizamos la función fzero de MATLAB para calcularla

...
f=@(x) -exp(-x)./x+1./(2*a*x.^2);
E=0.2; %energía
g=@(x) f(x)-E;
xm=fzero(g,[0.5,2]); %distancia mínima
...
>> xm
xm =    0.9692

Parámetro α=1.3

Representamos la función y(x) para el valor del parámetro α=1.3

a=1.3; %parámetro
f=@(x) -exp(-x)./x+1./(2*a*x.^2);
hold on
fplot(f,[0.5,5])
g=@(x) x*(x+1)*a-exp(x);
x_min=fzero(g,1);
plot(x_min,f(x_min),'ro','markersize',3,'markeredgecolor','r',
'markerfacecolor','r')
x_max=fzero(g,2);
plot(x_max,f(x_max),'bo','markersize',3,'markeredgecolor','b',
'markerfacecolor','b')

E=0.04; %energía
h=@(x) f(x)-E;
x0=fzero(h,[0.5,x_min]); %distancia mínima
line([x0,5],[E,E],'color','r')
line([x0,x0],[0,f(x0)],'lineStyle','--','color','k')

E=0.02; %energía
h=@(x) f(x)-E;
x1=fzero(h,[0.5,x_min]); %distancia mínima
x2=fzero(h,[x_min,x_max]); %distancia máxima
x3=fzero(h,[x_max,5]); %distancia mínima
line([x1,x2],[E,E],'color','r')
line([x3,5],[E,E],'color','r')
line([x1,x1],[0,f(x1)],'lineStyle','--','color','k')
line([x2,x2],[0,f(x2)],'lineStyle','--','color','k')
line([x3,x3],[0,f(x3)],'lineStyle','--','color','k')
ylim([0,0.05])
hold off
grid on
xlabel('x')
ylabel('y')
title('Potencial de Yukawa')

Observamos que la función y(x) presenta un mínimo y un máximo, puntos de color rojo y azul

El extremo de una función se calcula haciendo que su derivada primera sea nula

dy dx =0, x e x + e x x 2 1 α x 3 =0,x(x+1)α e x =0

El mínimo está cerca de x=1 y el máximo está cerca de x=2, ambos los calculamos con fzero de MATLAB

...
g=@(x) x*(x+1)*a-exp(x);
x_min=fzero(g,1);
x_max=fzero(g,2);
...
>> x_min,x_max
x_min =    1.1010
x_max =    2.2578
>> f(x_min)
ans =    0.0153
>> f(x_max)
ans =    0.0291

La derivada segunda de y(x) es

d 2 y d x 2 = x 2 +2x+2 x 3 e x +3 1 α x 4

Comprobamos que en un mínimo d2y/dx2>0 y que en un máximo d2y/dx2<0

>> h=@(x) -(x^2+2+x*2)*exp(-x)/x^3+3/(a*x^4);
>> h(x_min)
ans =    0.2214
>> h(x_max)
ans =   -0.0167

Para la energía E1=0.04, la recta horizontal E1=0.04 y la función y(x) se cortan en un punto que define la distancia mínima x0 al centro de fuerzas. La partícula se puede mover en el intervalo (x0, ∞)

...
f=@(x) -exp(-x)./x+1./(2*a*x.^2);
E=0.04; %energía
h=@(x) f(x)-E;
x0=fzero(h,[0.5,x_min]); %distancia mínima
...
>> x0
x0 =    0.7971

El caso más interesante se presenta cuando la energía E2 está entre la altura del mínimo, 0.0153 y la del máximo, 0.0291, por ejemplo, E2=0.02

Para esta energía, la recta horizontal E2=0.02 y la función y(x) se cortan en tres puntos: x1 situado en el intervalo (0, xmín), x2 situado en el intervalo (xmín, xmáx) y x3 situado en el intervalo (xmáx, ∞). La partícula se puede mover en el intervalo (x1, x2) o en el intervalo (x3,∞)

...
f=@(x) -exp(-x)./x+1./(2*a*x.^2);
E=0.02; %energía
h=@(x) f(x)-E;
x1=fzero(h,[0.5,x_min]); %distancia mínima
x2=fzero(h,[x_min,x_max]); %distancia máxima
x3=fzero(h,[x_max,5]); %distancia mínima
...
>> x1,x2,x3
x1 =    0.9360
x2 =    1.3958
x3 =    3.9150

Podríamos estar interesados en que el mínimo de y(x) fuese cero, para calcular el valor del parámetro α, hacemos que y(x)=0 y su derivada dy/dx=0 sean nulas simultáneamente

{ y(x)=0, e x x + 1 2α x 2 =0, e x 2αx=0 dy dx =0,x(x+1)α e x =0

La solución es x=1, a=e/2=1.3591

Para comprobarlo, se sugiere al lector sustituir en la primera línea de código del script, el valor de α=1.3 por α=1.3591

Punto de inflexión

Hemos comprobado que para α=0.9, no hay máximos y mínimos, para α=1.3 la curva presenta un máximo y un mínimo. Habrá un valor de α para el cual la función y(x) presenta un punto de inflexión, aquél en el que la derivada primera dy/dx=0 y la derivada segunda, d2y/dx2=0, se anulan simultáneamente

{ dy dx =0,x(x+1)α e x =0 d 2 y d x 2 =0,x( x 2 +2x+2 ) 3 α e x =0

Obtenemos la ecuación x2-x-1=0, cuya raíz positiva es

x= 1+ 5 2 =1.6180 α c = e x x(x+1) =1.1905

hold on
x0=(1+sqrt(5))/2; %punto de inflexión
a=exp(x0)/(x0*(x0+1)); %parámetro
f=@(x) -exp(-x)./x+1./(2*a*x.^2);
fplot(f,[0.5,5])
line([x0,x0],[0,f(x0)],'lineStyle','--','color','k')
line([0,x0],[f(x0),f(x0)],'lineStyle','--','color','k')
plot(x0,f(x0),'ro','markersize',3,'markeredgecolor','r','markerfacecolor','r')
ylim([0,0.1])
hold off
grid on
xlabel('x')
ylabel('y')
title('Potencial de Yukawa')

El potencial de Yukawa da lugar a trayectorias cerradas cuando

λ>1.1905 L 2 mk

Ecuación de la trayectoria

La lagrangiana es

L= 1 2 m ( dr dt ) 2 + 1 2 m r 2 ( dθ dt ) 2 + k r exp( r λ )

Como la lagrangina no depende de θ, hay una constante del movimiento (el momento angular)

d dt ( L θ ˙ ) L θ =0 d dt ( m r 2 θ ˙ )=0 m r 2 dθ dt =cte

La ecuación del movimiento en la dirección radial es

d dt ( L r ˙ ) L r =0 m d 2 r d t 2 mr ( dθ dt ) 2 + k r 2 exp( r λ )+ k λr exp( r λ )=0 d 2 r d t 2 = L 2 m 2 r 3 k mλ ( λ r +1 ) exp( r λ ) r

Donde L es el momento angular que se mantiene constante. Para calcular la ecuación de la trayectoria, r=r(θ), eliminamos el tiempo t del siguiente modo

dr dt = dr dθ dθ dt = L m r 2 dr dθ d 2 r d t 2 = d dt ( L m r 2 dr dθ )= d dθ ( L m r 2 dr dθ ) dθ dt = L m r 2 d dθ ( L m r 2 dr dθ ) d 2 r d t 2 = L 2 m 2 r 4 ( d 2 r d θ 2 2 r ( dr dθ ) 2 )

La ecuación diferencial de la trayectoria se escribe

L 2 m 2 r 4 ( d 2 r d θ 2 2 r ( dr dθ ) 2 )= L 2 m 2 r 3 k mλ ( λ r +1 ) exp( r λ ) r

Hacemos el cambio de variable u=1/r

dr dθ = 1 u 2 du dθ d 2 r d θ 2 = 2 u 3 ( du dθ ) 2 1 u 2 d 2 u d θ 2

La ecuación diferencial de la trayectoria se escribe

L 2 m 2 u 4 ( 2 u 3 ( du dθ ) 2 1 u 2 d 2 u d θ 2 2u ( 1 u 2 du dθ ) 2 )= L 2 m 2 u 3 k mλ u( λu+1 )exp( 1 λu ) d 2 (λu) d θ 2 +λu=α( λu+1 λu )exp( 1 λu ),α= mkλ L 2

Energía de la partícula

E= 1 2 m{ ( dr dt ) 2 + r 2 ( dθ dt ) 2 } k r exp( r λ ) E= 1 2 L 2 m r 4 { ( dr dθ ) 2 + r 2 } k r exp( r λ )

Denominamos x=r/λ

Eλ k = 1 2 L 2 mkλ { 1 x 4 ( dx dθ ) 2 + 1 x 2 } 1 x exp( x ) e= 1 2α { 1 x 4 ( dx dθ ) 2 + 1 x 2 } 1 x exp( x )

Dada la energía e calculamos la distancia mínima al centro de fuerzas. Por ejemplo, para el valor del parámetro α=1.3, si la energía de la partícula e=0.02, calculamos x1=0.9360, raíz de la ecuación transcendente con dx/dθ=0

Supongamos que la partícula parte de esa posición x1=r1/λ=1/(λu1)=0.9360. Resolvemos la ecuación diferencial por procedimientos numéricos con las siguientes condiciones iniciales, en la posición (λu1, θ=0), d(λu)/=0. En el código, denominamos x(1) a (λu) y x(2) a d(λu)/

a=1.3; %parámetro
e=0.02; %energía
h=@(x) 1./(2*a*x.^2)-exp(-x)./x-e;
x1=fzero(h,[0.5,1]); %mímima distancia
f=@(t,x) [x(2); -x(1)+a*(x(1)+1)*exp(-1/x(1))/x(1)]; 
[t,x]=ode45(f,[0,8*pi],[1/x1,0]);
plot(cos(t)./x(:,1), sin(t)./x(:,1))
grid on
axis equal
xlabel('x')
ylabel('y')
title('Potencial de Yukawa')

Comprobamos que la mínima distancia es 0.9360 y la máxima 1.3958

>> min(1./x(:,1))
ans =    0.9359
>> max(1./x(:,1))
ans =    1.3956    

Comprobamos que la energía e se mantiene constante e igual a 0.02

>> e=(x(:,2).^2+x(:,1).^2)/(2*a)-x(:,1).*exp(-1./x(:,1))
   e =
    0.0200
    0.0200
    ....
    0.0200
    0.0200

Para el valor del parámetro α=1.3, si la energía de la partícula e=0.04, calculamos la mínima distancia al centro de fuerzas, x0=0.7971. Resolvemos la ecuación diferencial por procedimientos numéricos con las siguientes condiciones iniciales, en la posición (λu0=1/x0, θ=0), d(λu)/=0

a=1.3; %parámetro
e=0.04; %energía
h=@(x) 1./(2*a*x.^2)-exp(-x)./x-e;
x0=fzero(h,[0.5,1]); %mímima distancia
f=@(t,x) [x(2); -x(1)+a*(x(1)+1)*exp(-1/x(1))/x(1)]; 
[t,x]=ode45(f,[0,2*pi-pi/4],[1/x0,0]);
plot(cos(t)./x(:,1), sin(t)./x(:,1))
grid on
axis equal
xlabel('x')
ylabel('y')
title('Potencial de Yukawa')

Referencias

O. L. de Lange, J. Pierrus. Solved Problems in Classical Mechanics. Analytical and numerical solutions with comments. Oxford University Press (2010). Question 8.6, pp. 224-228