ObjetosCompactos101404.html

MapleV, Equilibrio y Estabilidad de

Distribuciones Esféricas Autogravitantes:

Una Propuesta Pedagógica de Simulación de Situaciones en Astrofísica

Adriana Bisceglia
Centro Nacional de Cálculo Científico, Universidad de Los Andes (CeCalCULA).
Corporación Parque Tecnológico de Mérida (CPTM), Mérida, 5101. Venezuela
Correo-e: abisceg@ula.ve

Ángel G. Muñoz S.
Facultad de Ingeniería. Universidad Rafael Urdaneta.Maracaibo, 4004. Venezuela. y
Grupo de Investigaciones de Física Teórica (GIFT).
Depto. de Física. Facultad de Ciencias. La Universidad del Zulia, Maracaibo, 4004. Venezuela.
Correo-e: agmunoz@luz.edu.ve

Luis A. Núñez
Centro de Física Fundamental (CFF).
Departamento de Física, Facultad de Ciencias. Universidad de Los Andes, Mérida, 5101. Venezuela y
Centro Nacional de Cálculo Científico, Universidad de Los Andes (CeCalCULA).
Corporación Parque Tecnológico de Mérida (CPTM), Mérida, 5101. Venezuela
Correo-e: nunez@ula.ve
http://webdelprofesor.ula.ve/ciencias/nunez

Versión 1.1 Octubre de 2004

Resumen

Se presentan en esta hoja de Maple v.9 una propuesta para simulación del equilibrio y estabilidad de objetos compactos en Relatividad General. Para ello de utilizan los casos más conocidos de libros de texto en Astrofísica Relativista. Se ilustra con ellos modelos de estrellas polítropas, enanas blancas y estrellas de neutrones. La intención es ofrecer a los estudiantes la posibilidad de asomarse a los conceptos de Astrofísica Relativista y, al mismo tiempo, iniciarse en los conceptos de simulación de situaciones Físicas. Este documento evolucionará y sus versiones más recientes podrán ser encontradas en http://webdelprofesor.ula.ve/ciencias/nunez/

> restart:with(plots):

Warning, the name changecoords has been redefined

Introducción. Consideraciones Generales.

Antes de empezar es necesario introducir algunas nociones fundamentales para el estudio de equilibrio, estabilidad e inestabilidad de objetos autogravitantes. Nótese que emplearemos en lo que sigue c=1 , y consideraremos configuraciones isoentrópicas y con potencial electro-químico constante.

Como se sabe, idea básica es definir la energía total del sistema y a partir de ella encontrar los casos de equilibro estable e inestable por medio de los criterios usuales variacionales (equilibrio: primera variación de la energía igual a cero; estabilidad: segunda variación de la energía positiva). Un procedimiento equivalente sería encontrar la ecuación de Sturm-Liouville de perturbaciones radiales del sistema y ubicar los casos en los que los eigenvalores correspondientes son positivos.

Vale resaltar que aunque, en términos de la determinación de la estructura, la mayoría de las estrellas pueden ser descritas por medio de la teoría newtoniana de gravitación, surgen diferencias importantes cuando se estudia equilibrio y estabilidad por medio de la teoría de Newton y la de Einstein.

A continuación se presenta, para los casos newtoniano y relativista, la condición necesaria y suficiente para que la energía total del sistema autogravitante sea un extremal (ecuaciones de equilibrio hidrostático). Por simplicidad consideraremos el caso esféricamente simétrico para ambos estudios.

Caso Newtoniano

En este caso, la energía total del sistema se escribirá como (se supone que no existen movimientos macroscópicos)

> E:=Ein+Eg;

donde E[i] corresponde a la energía interna total y E[g] a la gravitacional. Denotemos por e=e(rho,s) la energía interna por unidad de masa. Ahora bien, si contamos con una ecuación de estado barótropa, es posible escribir e=e(rho(p(r)),s), que será útil a continuación para lo que nos proponemos. Vale recordar que, consideramos configuraciones isoentrópicas, de modo que podemos obviar la dependencia con la entropía s . La energía interna será, entonces:

> alias(e=e(rho(p(r)))):Ein:=Int(e,m=0..M);

y la gravitacional,

> Eg:=-G*Int(m/r,m=0..M);

La variación de E para perturbaciones radiales puede escribirse como

> delta(E)=diff(E,r)*delta(r);

pero, ya que (termodinámica)

> D(e)(rho(p(r)))=p*rho**(-2);

y (cálculo diferencial)

> D(rho)(p(r))=rho/p;

obtenemos

> delta(Energía)=Int((1/rho*diff(p(r),r)+G*m/r²)*delta.r,m=0..M);

Y si la variación de la energía es nula (extremal), entonces:

> diff(p(r),r)=-rho*G*m/r**2;

la cual se constituye en la ecuación newtoniana de equilibrio hidrostático.

Caso Relativista

Para nuestro propósito, de gran importancia serán la energía interna total de la estrella y la densidad propia de energía interna de materia.

La energía interna total de la estrella será:

> restart:E:=M-M[0];

donde

> M0:=m[N]*N;

es la energía que la materia de la estrella tendría si se dispersara hasta el infinito. m[N]=1.66e-24 gramos, es la masa en reposo de un nucleón y N es el número de nucleones en la estrella, y puede calcularse, si es necesario, mediante:

> N=Int(Int(Int(sqrt(g)*J°[N],r),theta),phi);

(donde J°[N] es la corriente conservada de número de nucleones y g es el determinante de la métrica), o bien mediante los parámetros métricos:

> N=Int(4*Pi*r**2*(J°N)*exp((lambda+nu)/2),r=0..R);

La masa total de la estrella, definida como su energía total más la de su campo gravitacional, se escribe como la función masa usual:

> M=int(4*Pi*r**2*rho(r),r=0..R);

Por otra parte, la densidad propia de energía interna de materia se escribe como

> e(r)=rho(r)-M0/N*n(r);

donde n(r) es la densidad propia de número por nucleón,

> n(r)=exp(nu/2)*J°[N];

Nótese que en virtud de e(r) es posible reescribir la energía interna total E como

> Eit:=T+G;

donde

> T:=int(4*Pi*r**2*exp(lambda(r)/2)*e(r),r=0..R);

corresponde a la energía térmica , mientras que

> G:=int(4*Pi*r**2*(1-exp(lambda(r)/2))*rho(r),r=0..R);

es la energía gravitacional . Estas expresiones nos serán de utilidad más adelante.

Ahora bien, en el caso relativista, una configuración estelar con entropía por nucleón y composición química uniformes estará en equilibrio si la energía interna total presentada arriba es un extremal, o equivalentemente, si la masa total M es estacionaria con respecto a todas la variaciones de rho(r) que dejen invariante el número de nucleones N. En términos de variaciones (empleando el método de los multiplicadores de Lagrange, sigma) esto se traduce a:

> restart:
v1:=delta*M-sigma*delta*N=4*Pi*int(r**2*delta(rho(r)),r=0..infinity)-sigma*4*Pi*int(r**2*(1-2*m(r)/r)**(-1/2)*delta(n(r)),r=0..infinity)-sigma*4*Pi*int(r*(1-2*m(r)/r)**(-3/2)*n(r)*delta(m(r)),r=0..infinity);

donde se ha hecho uso de la relación

> exp(-lambda)=1-2*m(r)/r;

Ahora bien,

> delta(m(r)):=4*Pi*int(r**2*delta(rho(r)),r=0..r);

y ya que las variaciones no cambian la entropía por nucleón, es cierto que

> delta(rho/n)+p*delta(1/n)=0;

de donde

> delta(n(r)):=delta(rho(r))*(n(r)/(p(r)+(rho(r))));

Así pues, introduciendo estos resultados y factorizando, se obtiene

> v2:=delta*M-sigma*delta*N=Int(4*Pi*r**2*(1-sigma*n(r)/(rho(r)+p(r))*(1-2*m(r)/r)**(-1/2)-sigma*Int(4*Pi*r*n(r)*(1-2*m(r)/r)**(-3/2),r=r..infinity)),r=0..infinity);

que se anulará si y sólo si la cantidad entre paréntesis es nula, de donde:

> Paréntesis:=1-sigma*n(r)/((p(r)+rho(r))*(1-2*m(r)/r)**(1/2))-sigma*int(4*Pi*r*n(r)/(1-2*m(r)/r)**(3/2),r =r..infinity)=0;

> sigma:=factor(solve(Paréntesis,sigma));

Este será el caso para algún multiplicador de Lagrange sigma si y sólo si el lado derecho es independiente de r, esto es, si su r-derivada es nula

> simplify(diff(sigma,r)=0);

De esta expresión se sigue que

> presquetov:=-diff(p(r),r)=collect(collect(collect(-solve(collect(simplify(diff(sigma,r)=0),diff(p(r),r)),diff(p(r),r)),m(r)),p(r)),rho(r));

Ahora bien, la condición de entropía por nucleón uniforme reza

> Diff(rho(r)/n(r),r)+p(r)*Diff(1/n(r),r)=0;

de manera que

> nd(r):=solve(diff(rho(r)/n(r),r)+p(r)*diff(1/n(r),r)=0,diff(n(r),r));

y, tras sustituir, obtenemos finalmente

> TOV:=r**2*((collect(simplify(subs(diff(n(r),r)=nd(r),diff(m(r),r)=4*Pi*r**2*rho(r),presquetov)),m(r))));

que es la ecuación de equilibrio hidrostático para distribuciones relativistas isótropas esféricamente simétricas, o, equivalentemente, toda configuración de materia que satisfaga esta ecuación poseerá una M estacionaria con respecto a variaciones de la densidad que dejen estacionaria N. Esta ecuación se conoce bajo el nombre de ecuación de Tolman-Oppenheimer-Volkov para fluidos pascalianos o isótropos .

Por otra parte, y en términos de auto-consistencia, nótese que si G = 1, p << rho , 2*G*m << r , la ecuación TOV se reduce a

> tov2:=subs(p(r)=p,rho(r)=rho,m(r)=m,TOV);
newt:=-r**2*Diff(p(r),r)=convert(convert(series(series(rhs(tov2),p=0,1),m=0,2),polynom),polynom);

la cual, es la ecuación newtoniana de equilibrio hidrostático. Siguiendo la misma lógica es posible construir una ecuación de equilibrio hidrostático para fluidos no pascalianos (materia anisótropa). Esto es una ecuación de equilibrio para fluidos con presiones distintas en distintas direcciones. La ecuación de equilibrio de Tolman-Oppenheimer-Volkov para fluidos no-pascalianos o anisótropos puede escribirse como

> Diff(P[rad](r),r)=-(rho(r) + P[rad](r))*((m(r)-4*Pi*r^3*P[rad](r))/(r*(r-2*m(r)))) -(2/r)*(P[tan](r)-P[rad](r));

la cual claramente si se reduce a la ecuación TOV para fluidos pascalianos, esto es.

> Diff(P[rad](r),r)=-(rho(r) + P[rad](r))*((m(r)-4*Pi*r^3*P[rad](r))/(r*(r-2*m(r))));

La expresión para la masa viene descrita por como

> m(r)=Int(4*Pi*r^2*rho(r),r=0..R);

o equivalentemente

> diff(m(r),r)=4*pi*r^2*rho(r);

Es claro que para todos los caso (Newtoniano, Relativista Isótropo y Relativista Anisótropo) se requiere de una ecuación de estado. Vale decir una ecuación que ligue la presión con la densidad, para el caso Isótropo y para el caso anisótropo. Estas ecuaciones se denominan ecuaciones de estado y, como ligan la presión con la densidad, se llaman ecuaciones de estado barótropas. Una de las ecuaciones de estado más utilizadas se conoce como ecuación de estado Polítropa. Ella relaciona la presión con una potencia de la densidad: , En este caso tendremos para el caso Newtoniano

> NewtonHidroEstatica1:= m(r)= (r^2/rho(r))*diff(P(r),r);
diff(%,r);
NewtonHidroEstatica2:=simplify(4*pi*r^2*rho(r)=2*r/rho(r)*diff(P(r),r)-r^2/rho(r)^2*diff(P(r),r)*diff(rho(r),r)+r^2/rho(r)*diff(P(r),`$`(r,2)));

Si suponemos la ecuación de estado Polítropa tendremos que

> P(r):=K*rho(r)**gamma1;
simplify(NewtonHidroEstatica2);

Es claro que en este caso (que es uno de los más simples) la ecuación diferencial que emerge no es nada sencilla. Es una ecuación diferencial ordinaria, NO lineal de segundo orden. En la próxima sección estudiaremos, en detalle las configuaciones que emegen de esta ecuación.

Estrellas Newtonianas polítropas

En varios escenarios astrofísicos, sobre todo newtonianos, ocurre que la densidad de energía interna es proporcional a la presión:

> restart:with(plots):
e(r):=rho(r)-m[N]*n(r)=(gamma1-1)**(-1)*P(r);

Warning, the name changecoords has been redefined

La idea de expresar la presión como función de una potencia de la densidad surge de la condición de entropía uniforme por nucleón la cual nos dice que

> Diff(rho(r)/n(r),r)+P(r)*Diff(1/n(r),r)=0;
Diff(rhs(e(r))/n(r),r)+P(r)*Diff(1/n(r),r)=diff(rhs(e(r))/n(r),r)+P(r)*diff(1/n(r),r);

de donde, se puede resolver el lado derecho y se obtiene

> dsolve(rhs(%),P(r));

con _C1 una constante de integración.

En astrofísica newtoniana, la energía interna (y la presión) son mucho menores que la densidad de masa en reposo, de manera que puede escribirse

> n(r):=rho(r)/m[N];

y así

> P(r)=K*rho(r)**gamma1;

Cualquier estrella con una ecuación de estado como la de la ecuación anterior se denomina polítropa . En general, un proceso politrópo es uno cuasi-estático en el cual el calor específico permanece constante (dQ/dT=cte). Estos modelos son muy versátiles puesto que permiten describir una gran variedad de escenarios utilizando una misma ecuación de estado y variando el único parámtro libre: el exponente . Como veremos, esencialmente las enanas blancas se modelan con ecuaciones de estado polítropas, así como también casos límite de estrellas de neutrones y también estrellas supermasivas. Más adelante extenderemos este concepto, , para como una ecuación de estado para describir cualquier configuración material, relativista o no.

Para determinar la estructura del objeto autogravitante tenemos, como única fuerza la gravedad por unidad de volumen

> F[grav]:=-G*m(r)*rho(r)/r^2;

es decir, un volumen diferencial de materia estará sometida a esta fuera por la masa contenida desde el origen ( r=0 ) hasta el punto donde se encuentra el volumen ( ) de ella derivamos la ecuación de ecuación de equilibrio hidrostático

> EcHidroEstNew:=diff(P(r),r)=F[grav];

y la definición de función masa contenida hasta ese volumen diferencial.

> ec2:=diff(m(r),r)=4*pi*r^2*rho(r);

Combinando las ecuaiones: de estado polítropa, la expresión para la Fuerza de Gravitación y la expresión para la masa, es posible obtener la siguiente ecuación diferencial de segundo orden:

> 4*Pi*r^2*rho(r) = r*K*gamma1*(2*rho(r)^gamma1*diff(rho(r),r)-2*r*rho(r)^(gamma1-1)*diff(rho(r),r)^2+r*rho(r)^(gamma1-1)*gamma1*diff(rho(r),r)^2+r*rho(r)^gamma1*diff(rho(r),`$`(r,2)))/rho(r)^2;

que no es más que la correspondiente ecuación de Poisson.

Para manipular una ecuación diferencial resulta conveniente introducir una variable adimensional de distancia tal que :

> r=A[gamma1]*xi;

con

> A[gamma1]:=sqrt( (K*gamma1/(4*Pi*(gamma1-1)))) *rho0^((gamma1-2)/2);

donde es la densidad central de la configuración autogravitante. Adicionalmente, con la utilización de la regla de la cadena podemos traducir el operador diferencial

> diff(f(r),r)=diff(f(r(xi)),xi);

el radio, la densidad y la presión vienen dado por:

> rho(r):=rho0*Theta(xi)^(1/(gamma1 -1));
P(r):=K*(rho0^gamma1)*(Theta(xi)^(gamma1/(gamma1-1)));

Con estas nuevas variables se puede transformar la ecuación diferencial de hasta llegar a la ecuación de Lane-Emden la cual, puede ser escrita:

> LaneEmden:=2*diff(Theta(xi),xi)+xi*diff(Theta(xi),`$`(xi,2))+Theta(xi)^(1/(gamma1-1))*xi=0;

con las condiciones iniciales y

Históricamente, la manera (Chandrasekhar 1939) de resolver esta ecuación fue mediante el método de series de potencias y la serie solución (como era de suponerse) se el conoce con el nombre de funciones de Lane-Emden . Resolviendo la ecuación mediante el método de serie de potencias, tendremos que:

> Order:=10:
dsolve({LaneEmden,Theta(0)=1, D(Theta)(0)=0},Theta(xi),series);
ThetaGammaSerie:= rhs(%):
ThetaGammaPoly:= convert(%,polynom):

El enfoque tradicional (Chandrasekhar 1939) es estudiar esta serie para distintos valores de .

No es nada trivial darse cuenta que para valores se encuentran valores , tales que , con lo cual encontramos el radio de la distribución.

explorando el comportamiento de la serie para distintos valores de no resulta intuitivo.

> Thetagamma65:=subs(gamma1=6/5,ThetaGammaSerie);
Thetagamma65poly:= convert(%,polynom):
Thetagamma119:=subs(gamma1=11/9,ThetaGammaSerie);
Thetagamma119poly:= convert(%,polynom):
Thetagamma54:=subs(gamma1=5/4,ThetaGammaSerie);
Thetagamma54poly:= convert(%,polynom):
Thetagamma97:=subs(gamma1=9/7,ThetaGammaSerie);
Thetagamma97poly:= convert(%,polynom):
plot([Thetagamma65poly,Thetagamma119poly,Thetagamma54poly,Thetagamma97poly],xi=1..2,
labels=["xi","Theta"],labeldirections=[HORIZONTAL,VERTICAL], title="Theta vs xi",
legend=["gamma = 6/5","gamma = 11/9","gamma = 5/4","gamma = 9/7"]);

Aún extendiendo la expansión en serie

> Order:=20:
dsolve({LaneEmden,Theta(0)=1, D(Theta)(0)=0},Theta(xi),series);
ThetaGammaSerie:= rhs(%):
ThetaGammaPoly:= convert(%,polynom):

> Thetagamma65:=subs(gamma1=6/5,ThetaGammaSerie);
Thetagamma65poly:= convert(%,polynom):
Thetagamma119:=subs(gamma1=11/9,ThetaGammaSerie);
Thetagamma119poly:= convert(%,polynom):
Thetagamma54:=subs(gamma1=5/4,ThetaGammaSerie);
Thetagamma54poly:= convert(%,polynom):
Thetagamma97:=subs(gamma1=9/7,ThetaGammaSerie);
Thetagamma97poly:= convert(%,polynom):
plot([Thetagamma65poly,Thetagamma119poly,Thetagamma54poly,Thetagamma97poly],xi=0.1..5,
y=-1..1, labels=["xi","Theta"],labeldirections=[HORIZONTAL,VERTICAL], title="Theta vs xi",
legend=["gamma = 6/5","gamma = 11/9","gamma = 5/4","gamma = 9/7"]);

Los comportamientos de las series son muy distintos al aumentar el número de términos. Algo muy distinto sucede si resolvemos numéricamente la ecuación de Lane-Emden para distintos valores del parámetro . Vale decir, = 6/5, 11/9, 5/4 y 9/7

> LaneEmden65:=subs(gamma1=6/5,LaneEmden);
solnumerica65:= dsolve({LaneEmden65,Theta(0)=1,D(Theta)(0)=0},Theta(xi),numeric,range=0..50):
plotsolnumerica65:= odeplot(solnumerica65,color=blue):
LaneEmden119:=subs(gamma1=11/9,LaneEmden);
solnumerica119:=dsolve({LaneEmden119,Theta(0)=1,D(Theta)(0)=0},Theta(xi),numeric,range=0..30):
plotsolnumerica119:=odeplot(solnumerica119,color=red):
LaneEmden54:=subs(gamma1=5/4,LaneEmden);
solnumerica54:=dsolve({LaneEmden54,Theta(0)=1,D(Theta)(0)=0},Theta(xi),numeric,range=0..50):
plotsolnumerica54:=odeplot(solnumerica54,color=green):
LaneEmden97:=subs(gamma1=9/7,LaneEmden);
solnumerica97:=dsolve({LaneEmden97,Theta(0)=1,D(Theta)(0)=0},Theta(xi),numeric,range=0..9.5):
plotsolnumerica97:=odeplot(solnumerica97, color=yellow):
display([plotsolnumerica65,plotsolnumerica119,plotsolnumerica54,plotsolnumerica97],view=[0..50,0..1]);

podemos explorar esos valores de para los cuales se anula . Para el caso = 11/9, estudiamos los valores de alrededor de = 30

> solnumerica119(31);
solnumerica119(31.8);
solnumerica119(31.9);

igualmente, procedemos explorar valores de para los cuales se anula . Para el caso = 5/4, estudiamos los valores de alrededor de = 15

> solnumerica54(14.97);
solnumerica54(14.98);

Lo importante de estos puntos de corte es que constituyen puntos en los cuales se anulan también las variables físicas (presión, o densidad) y el valor se relaciona con la frontera de la distribución. Con ello hemos encontrado los radios correspondientes a distribuciónes autogravitantes estables con ecuaciones de estado politrópas del la forma ,

Para calcular la masa total de la distribució expresamos la densidad en términos de las funciones de a

> M=Int(4*Pi*r^2*rho(r),r=0..R);
M=4*Pi*(1/2/Pi^(1/2)*(K*gamma1/(gamma1-1))^(1/2)*rho0^(1/2*gamma1-1))^2*rho0*((K*gamma1/(gamma1-1))^(1/2)*rho0^(1/2*gamma1-1))*Int(xi^2*Theta(xi)^(1/(gamma1-1)),xi=0..xi[1]) ;

Utilizando la ecuación de Lane Emden para sustituir la expresión para tendremos

> -(1/xi^2)*Diff(xi^2*Diff(Theta(xi),xi),xi)=Theta(xi)^(1/(gamma1-1));

> M=4*Pi*(1/2/Pi^(1/2)*(K*gamma1/(gamma1-1))^(1/2)*rho0^(1/2*gamma1-1))^2*rho0*((K*gamma1/(gamma1-1))^(1/2)*rho0^(1/2*gamma1-1))*Int(xi^2*((-1/xi^2)*Diff(xi^2*Diff(Theta(xi),xi),xi)),xi=0..xi[1]) ;

> M=(gamma1/(gamma1-1)*(rho0^(1/2*gamma1-1))^3*rho0*(K*gamma1/(gamma1-1))^(1/2))*(abs(xi^2*Diff(Theta(xi),xi))[xi[1]]);

A continuación se presenta, a modo ilustrativo, un procedimiento para resolverla numéricamente.

Lane-Emden

> restart;
with(plots):

Warning, the name changecoords has been redefined

Proporciónese el valor del índice polítropo:

> n:=3;
alpha:=0.1:

Este "procedimiento" de Maple informa detalladamente cómo hacer el cambio de variable para producir el sistema de ecs. de primer orden:

>    solproc := proc(N, x, Y, YP)
    YP[1] := Y[2];
    YP[2] := -Y[1]^n-2/x*Y[2];
  end:

Ahora se introducen las condiciones iniciales en forma de arreglo y se especifican las variables

> ci := array([1,0]):
dvars := [y(x),z(x)]:

El cálculo se asigna en memoria interna a "LE". Nótese que, debido a que en x=0 la ecuación es singular, comienzo con un número bastante cercano, pero distinto de 0. OJO : Ha de variarse la cota superior del rango de acuerdo al índice polítropo escogido, de manera que justo coincida con la visualización del punto del primer corte con el eje radial. Esto es importante debido a que Maple almacena en memoria todas las soluciones (con interpolación de orden 4 entre los puntos) y un exceso en el rango consume recursos innecesariamente.

Escriba la cota superior del rango de integración a utilizar:

> cota:=7:

> LE:=dsolve(numeric,number=2,method=rkf45,output=listprocedure,procedure=solproc,range=1e-16..cota,initial=ci,procvars=dvars);

> odeplot(LE,[[x,y(x)],[x,z(x)]],0..cota,axes=boxed,labels=[" x "," Y , D Y "],labeldirections=[HORIZONTAL,VERTICAL], title="Soluc. Lane-Emden (Newtoniana)",legend=[" Y ","D Y "]);

Warning, cannot evaluate the solution further left of .19894130e-70, probably a singularity

Escribimos las variables de interés como "procedimientos" maple.

> psi:=eval(y(x),LE);xi:=eval(x,LE);

Determinamos el radio xi1 de la polítropa (aquí "err" es nuestro "0 numérico"):

> err:=1e-4:for i from cota-1 by 1e-3 while psi(i)>=err do end do;xi1:=i;

Recuperamos las variables físicas:

> rhoc:=alpha*A^2*(n+1)/4/Pi;

> rho:=simplify(rhoc*psi^n);

> K:=alpha/rhoc^(1/n);

> P:=simplify(K*rho^(1+1/n));

> A:=2;

> plot([P,rho],0..cota,0..0.16,title="Variables Físicas.\n Polítropa con n=3",axes=boxed, labels=[ x ,("P, rho")],labeldirections=[HORIZONTAL,VERTICAL],legend=["P","rho"]);

Por último, incluiremos aquí que, a partir de la función masa introducida en la primera sección, la masa total de la estrella queda escrita según

> restart:M:=4*Pi*Int(r**2*rho,r)=-4*Pi*((n+1)*K/4/Pi)**(3/2)*rhoc**((3-n)/2/n)*(xi**2*Diff(Psi,xi))[xi1];

Equilibrio y Estabilidad.

Ahora obtengamos los casos de equilibrio.

Las energías térmica y gravitacional, tras emplear la ecuación de equilibrio hidrostático newtoniana e integrar por partes [Weinberg], son (ver primera sección)

> restart:T:=1*M**2/R*1/(5*gamma-6);V:=-3*(gamma-1)/((5*gamma-6))*M**2/R;

y la energía interna total es

> E:=simplify(T+V);

Empleando el criterio de la primera derivada, notamos que el punto crítico corresponde a un valor de gamma que satisfaga

> dE:=diff(E,R)=0;de:=diff(subs(gamma=gam,E),R):

o bien

> gamma=solve(de=0,gam);

que es el valor también en el que la energía total se hace exactamente nula, independientemente del radio de la estrella. Se aprecia también inmediatamente que esta expresión presenta una singularidad para

> gamma=6/5;

> M:=0.25*R:R:=9e6:plot([E,0],gamma=0..5,-1e6..1e5,axes=boxed,labels=["gamma","E"],labeldirections=[HORIZONTAL,VERTICAL], title="Energía Interna Total de una configuración polítropa con M/R=0.25 y R=90km");

En aras de determinar realmente el equilibrio y estabilidad de estas configuraciones vale la pena recordar que la condición necesaria y suficiente es que E sea un mínimo con respecto a todas las variaciones de densidad que dejen N invariable. De las ecuaciones presentadas en la primera sección se obtiene que:

> restart:ene:=N=4*Pi/3/mn*rho*R**3;T:=4*Pi/3*(gamma-1)**(-1)*K*rho**gamma*R**3;V:=-16*Pi**2/15*rho**2*R**5;

De manera que la energía interna total es

> E:=simplify(T+V);

Eliminando R, y teniendo en cuenta que para estrellas newtonianas la masa M está dominada por la masa en reposo total Nmn,

> R:=(3*M/4/Pi/rho)**(1/3);

obtenemos, finalmente para este caso,

> E:=collect(E,rho);

Tras derivar esta expresión con respecto a la densidad y anularla, nos percatamos de que el valor de gamma que satisface el estado de equilibrio ha de ser

> dE:=diff(subs(gamma=gam,E),rho)=0:g:=factor(solve(dE,gam))=4/3;

que, según vimos antes, debe coincidir con el valor 4/3. Inmediatamente podemos despejar la masa, obteniendo (presentamos la solución positiva)

> Ms:=solve(g,M):M[c]=Ms[2];

que es la masa que ha de poseer la configuración para que E sea estacionaria. Nótese que por lo recién visto, una estrella con esta masa posee una energía interna total E nula. Esta masa corresponde al caso marginal de paso de configuración estable a inestable (gamma=4/3), esto es, es una masa crítica.

Calculemos finalmente a qué valores de gamma corresponden los casos de estabilidad e inestabilidad. Ya sabemos que el equilibrio se determina, para M constante, anulando la primera derivada de la enrgía interna total con respecto a la densidad:

> dE:=diff(E,rho)=0;e1:=solve(dE,M):

que corresponde a una masa igual a

> M:=simplify(e1[1]);

El estado de equilibrio estable corresponderá a los mínimos de tal masa con respecto variaciones de la densidad, mientras que el equilibrio inestable corresponderá a los máximos. Así pues, notamos que

> f1:=factor(diff(M,rho)):Diff(Mas,rho)=f1;

es proporcional a

> -4+3*gamma;

De modo que los casos de estabilidad corresponden a

> es:=f1>0;

o a un valor de gamma en el intervalo:

> assume(gam>0):solve(subs(gamma=gam,es),gam):%[2];

Análogamente, los casos de inestabilidad corresponden a

> ines:=f1<0;

o a un valor de gamma en el intervalo:

> a1:=subs(gamma=gam,ines):solve(a1,gam):%[1];

Enanas Blancas

Cuando las estrellas han consumido su combustible nuclear y comienzan a enfriarse y contraerse, algunas se convierten en lo que genéricamente se denominan "enanas blancas": estrellas constituidas fundamentalmente por electrones y cuya estructura se sostiene por medio de presión de degeneración electrónica, impuesta por el Principio de Exlusión de Pauli.

Efectivamente, para temperaturas suficientemente bajas, los electrones se ubican en los menores niveles de energía disponibles, y el principio de Pauli entonces nos dice que podrá haber sólo dos electrones con espines orientados en distintos sentidos por cada nivel.

Ya que existen

> restart:4*Pi*k**2*(2*Pi*hb)**(-3)*dk;

niveles por unidad de volumen con momenta entre k y k+dk, el número de electrones por unidad de volumen será (recordar que hay que multiplicar por 2 la expresión para tomar en cuenta los 2 electrones por nivel):

> n=8*Pi*(2*Pi*hb)**(-3)*Int(k**2,k=0..k[F]);n:=8*Pi*(2*Pi*hb)**(-3)*int(k**2,k=0..k[F]);

La densidad de masa será

> e1:=rho=n*m[N]*mu;

> e2:=solve(e1,k[F]):k[F]:=factor(e2[1]);

La condición, mencionada brevemente arriba, que indica que la temperatura es despreciable y puede escribirse como (el signo < corresponde a <<)

> k*T<(k[F]**2+m[e]**2)**(1/2)-m[e];

La energía cinética y la presión de los electrones son, respectivamente,

> e:=8*Pi/(2*Pi*hb)**3*Int(((k**2+me**2)**(1/2)-me)*k**2,k=0..k[F]);

> p:=8*Pi/3/(2*Pi*hb)**3*Int(1/(k[F]**2+me**2)**(1/2)*k**4,k=0..k[F]);

Y la ecuación de estado general correspondiente a este caso se puede obtener integrando directamente la expresión anterior.

A pesar de que la ecuación de estado resulta un poco complicada de tratar, existen dos casos particulares y extremos en los que es posible obtener expresiones simplificadas. Los casos corresponden a densidad estelar (a) mucho menor y (b) mucho mayor que una cierta densidad crítica , tal que el momentum máximo k[F] se hace igual a me , y de valor numérico:

> rho[c]:=1/3*me^3*m[N]*mu/(Pi^2*hb^3)=0.97*1e6*mu;

donde las unidades de esta última ecuación corresponden a g/cm**3.

CASO 1. Densidad<<Densidad crítica.

> restart:p:=8*Pi/3/(2*Pi*hb)**3*Int(k**4/me,k=0..k[F])=8*Pi/3/(2*Pi*hb)**3*int(k**4/me,k=0..k[F]);

> k[F] := 3^(1/3)*(rho*Pi^2*m[N]^2*mu^2)^(1/3)*hb/(m[N]*mu);

de manera que

> p:=rhs(p);

que es una ecuación polítropa con gamma=5/3 y

> K:=simplify(1/5*hb^2*3^(2/3)*(Pi^2*m[N]^2*mu^2)^(5/3)/(Pi^2*me*m[N]^5*mu^5),symbolic);

Y al calcular la función masa correspondiente (ver final de la sección de las polítropas) se obtiene

> rhoc:=1/3*me^3*m[N]*mu/(Pi^2*hb^3):M:=1/2*(3*Pi/8)**(1/2)*(2.71)*(hb**3/(m[N]**2*mu**2))*(rho(0)/rhoc)**(1/2)=2.79*mu**(-2)*(rho(0)/rhoc)**(1/2)*Ms;

donde Ms es la masa del Sol.

El radio será

> R:=2.0*1e4*mu**(-1)*(rho(0)/rhoc)**(-1/6);

el cual viene ya dado en km.

CASO 2. Densidad>>Densidad crítica.

En este caso tenemos que k[F]>>m[e], y así

> p:=8*Pi/3/(2*Pi*hb)**3*Int(k**3,k=0..k[F])=8*Pi/3/(2*Pi*hb)**3*int(k**3,k=0..k[F]);

de manera que tenemos una polítropa con gamma=4/3 y

> K:=1/4*hb*3^(1/3)*(Pi^2*m[N]^2*mu^2)^(4/3)/(Pi^2*m[N]^4*mu^4);

Es importante notar que debido al valor del índice polítropo, se obtiene una masa única de valor

> M:=1/2*(3*Pi)**(1/2)*(2.018)*(hb**(3/2)/(m[N]**2*mu**2))=5.87*mu**(-2)*Ms;

donde, nuevamente, Ms es la masa del Sol.

El radio será

> R:=1/2*(3*Pi)**(1/2)*(6.896)*(hb**(3/2)/(m[e]*m[N]*mu))*(rho(0)/rhoc)**(1/2)=5.3*1e4*mu**(-1)*(rhoc/rho(0))**(1/3);

dado en km.

Conclusiones

Vale la pena detenerse un instante a interpretar estos resultados.

Para Densidad central<<Densidad crítica , nos percatamos de que gamma>4/3 (ver CASO 1), de modo que las enanas blancas menos masivas son ESTABLES.
Nos percatamos de que la masa M de la estrella tiende a crecer monótonamente a medida que aumenta la densidad central.
Para Densidad central>>Densidad crítica , o bien Densidad central==>infinito, la configuración posee una masa máxima, y la estructura continúa siendo estable. Esta masa máxima se conocer como LÍMITE DE CHANDRASEKHAR. Si se escoge mu=56/26 [Weinberg], el valor es de:

> mu:=56/26;M[Chandrasekhar]:=M;

NO HAY un valor de densidad central para el cual las estrellas puedan ser INESTABLES, al menos hasta la Masa de Chandrasekhar.
Cuando k[F] es del orden de 5m[e], es energéticamente favorable para los electrones unirse a protones y anti-neutrinos y convertirse en neutrones. De manera que llegada una cierta densidad, se incrementará el número de nucleones por electrón y es así como a partir de aproximadamente

> rho(0)=5**3*rhoc;

la masa M alcanza un máximo y comienza a decrecer. Así pues, el resultado más importante es que:

Las estrellas de neutrones sólo existen de manera estable para

> M<M[Chandrasekhar];

Esferas Relativistas Isótropas

> restart;with(plots):

Warning, the name changecoords has been redefined

Partimos de la ecuación para de TOV para esferas relativistas isótropas.

> EcuacTOVIso:=Diff(P(r),r)=-(rho(r) + P(r))*((m(r)-4*Pi*r^3*P(r))/(r*(r-2*m(r))));

Incluimos la definición de masa

> m(r):=int(4*Pi*r^2*rho(r),r);eval(EcuacTOVIso);

Adimensionalizando con y se obtiene la versión de la ecuación de TOV adimensional

>

> EcuacTOVIsoAdim:= diff(P(zeta),zeta)/R = -(rho1(zeta)*rho0+P(zeta))*(int(4*Pi*(R*zeta)^2*rho1(zeta)*rho0*R,zeta)-4*Pi*(R*zeta)^3*P(zeta))/(R*zeta)/((R*zeta)-2*int(4*Pi*(R*zeta)^2*rho1(zeta)*rho0*R,zeta));

La cual constituye la versión adimensional de la ecuación TOV para fluidos isótropos. Ahora consideraremos varias ecuaciones de estado

Esferas Relativistas de densidad constante

Consideremos el caso más simple desde el punto de vista matemático pero que no está exento de interés físico. Las configuraciones autogravitantes de densidad uniforme constituyen el caso de esferas incompresibles y es el caso límite para este tipo de objetos. Para el caso adimensional que nos compete la expresión a utilizar será

> rho1(zeta):=1;simplify(eval(EcuacTOVIsoAdim));

Simplificando aún más hacemos

> rho0:=M/((4/3)*Pi*R^3);
M:=eta*R;
P(zeta):=P1(zeta)/(4*Pi*R^2);
simplify(eval(EcuacTOVIsoAdim));

y resolvemos

> dsolve(%):P1(zeta):=rhs(%);

donde es una constante de integración que se determina a partir considerando el punto, frontera de la distribución donde la presión se anula.

> subs(zeta=1,rhs(%%)):
_C1:= solve(%,_C1);

> plotP1(zeta):=simplify(eval(P1(zeta)));

> animate(plot,[plotP1(zeta),zeta=0..1],eta=0.1..0.49, frames=10);

Esferas Relativistas con densidad Modelo Tipo Tolman VI

Consideremos las siguiente densidad

> rho1(zeta):=Kappa/zeta^2;

con lo cual la ecuación de equilibrio hidróstático queda como

> simplify(eval(EcuacTOVIsoAdim));

Simplificando aún más hacemos

> rho0:=M/((4/3)*Pi*R^3);
M:=eta*R;
P(zeta):=P1(zeta)/(4*Pi*R^2);
simplify(eval(EcuacTOVIsoAdim));

y resolvemos

> dsolve(%):P1(zeta):=rhs(%);

Error, (in dsolve) expecting an ODE or a set or list of ODEs. Received 3*zeta*eta^2*(-1+2*eta)/(-2+eta+3*zeta^2*eta)^2/Pi/R^3 = -9/4*eta^2*(-2*Kappa+Kappa*eta+3*Kappa*zeta^2*eta-zeta^2*eta+zeta^4*eta)*(2*Kappa-Kappa*eta-3*Kappa*zeta^2*eta-zeta^2*eta+zeta^4*eta)/zeta^3/R^3/Pi/(-2+eta+3*zeta^2*eta)^2/(-1+6*Kappa*eta)

donde es una constante de integración que se determina a partir considerando el punto, frontera de la distribución donde la presión se anula.

> subs(zeta=1,rhs(%%)):
_C1:= solve(%,_C1);

luego lo sustituimos

> simplify(eval(P1(zeta)));

En este caso hay proveer el parámetro K

> plotP11:=subs(Kappa=1,eval(P1(zeta)));

> animate(plot,[plotP11(zeta),zeta=0..1],eta=0.1..0.49, frames=10);

Plotting error, empty plot

Esferas Relativistas Anisótropas

Iniciamos esta sección de un modo distinto. Partimos de una métrica del espacio tiempo

de la métria del espacio tiempo se derivan, utilizando las ecuaciones de Einstein las variables físicas

, y

Nos damos cuenta que este tipo de métrica genera una ecuación de estado (relación funcional entre la presión y densidad) no local, en el sentido que la presión y la densidad no sólo están relacionadas por propiedades de la materia local sino que heredan en parte el comportamiento de todo el interior de la configuración material.

Cláculo de las Variables Físicas

> restart;grtw();

> makeg(nolocal);

Makeg 2.0: GRTensor metric/basis entry utility

To quit makeg, type 'exit' at any prompt.

Do you wish to enter a 1) metric [g(dn,dn)],

                       2) line element [ds],

                       3) non-holonomic basis [e(1)...e(n)], or

                       4) NP tetrad [l,n,m,mbar]?

makeg> 2;

Enter coordinates as a LIST (eg. [t,r,theta,phi]):

makeg> [t,r,theta,phi];

Enter the line element using d[coord] to indicate differentials.

(for example,  r^2*(d[theta]^2 + sin(theta)^2*d[phi]^2)

[Type 'exit' to quit makeg]

 ds^2 =

makeg> (C*d[t]^2 -d[r]^2)/h(r) -r^2*(d[theta]^2 + sin(theta)^2*d[phi]^2);

If there are any complex valued coordinates, constants or functions

for this spacetime, please enter them as a SET ( eg. { z, psi } ).

Complex quantities [default={}]:

makeg> ;

{}

You may choose to 0) Use the metric WITHOUT saving it,

                  1) Save the metric as it is,

                  2) Correct an element of the metric,

                  3) Re-enter the metric,

                  4) Add/change constraint equations,

                  5) Add a text description, or

                  6) Abandon this metric and return to Maple.

makeg> 1;

Information written to: `C:/Grtii(6)/Metrics/nolocal.mpl`

Do you wish to use this spacetime in the current session?

(1=yes [default], other=no):

makeg> 1;

Initializing: nolocal

makeg() completed.

> grcalcalter(G(up,dn),factor,simplify);grcalc(R(up,dn,dn,dn),factor,simplify);

Created definition for G(up,dn)

Simplification will be applied during calculation.

Applying routine factor to object g(up,up)

Applying routine factor to object g(dn,dn,pdn)

Applying routine factor to object Chr(dn,dn,dn)

Applying routine factor to object Chr(dn,dn,up)

Applying routine factor to object R(dn,dn)

Applying routine factor to object tRicciscalar

Applying routine factor to object G(dn,dn)

Applying routine factor to object G(up,dn)

Applying routine simplify to object g(up,up)

Applying routine simplify to object g(dn,dn,pdn)

Applying routine simplify to object Chr(dn,dn,dn)

Applying routine simplify to object Chr(dn,dn,up)

Applying routine simplify to object R(dn,dn)

Applying routine simplify to object tRicciscalar

Applying routine simplify to object G(dn,dn)

Applying routine simplify to object G(up,dn)

Created definition for R(up,dn,dn,dn)

> grdisplay(G(up,dn));

> m(r):=r*simplify(grcomponent(R(up,dn,dn,dn),[phi,theta,phi,theta]))/2;

> rho(r):=simplify(grcomponent(G(up,dn),[t,t])/(8*pi));

> Prad(r):=-grcomponent(G(up,dn),[r,r])/(8*pi);

> Ptan(r):=-grcomponent(G(up,dn),[theta,theta])/(8*pi);

La Ecuación de Equilibrio Hidrostático

> EcuacHidroAni:=diff(Prad(r),r)=-(rho(r) + Prad(r))*((m(r)-4*pi*r^3*Prad(r))/(r*(r-2*m(r)))) -(2/r)*(Ptan(r)-Prad(r));

> simplify(EcuacHidroAni);

> dsolve(%);

$h(r) = `&where`(_a,[{diff(_b(_a),_a) = (-2+2*_a)*_b(_a)^3-(4*_a-1)/_a*_b(_a)^2}, {_a = h(r), _b(_a) = 1/(r*diff(h(r),r))}, {h(r) = _a, r = exp(Int(_b(_a),_a)+_C1)}])$
$h(r) = `&where`(_a,[{diff(_b(_a),_a) = (-2+2*_a)*_b(_a)^3-(4*_a-1)/_a*_b(_a)^2}, {_a = h(r), _b(_a) = 1/(r*diff(h(r),r))}, {h(r) = _a, r = exp(Int(_b(_a),_a)+_C1)}])$

> EcuacHidroAniAdim:=-1/8*((diff(h(zeta),`$`(zeta,2))/R^2)*(zeta*R)^2-2*(zeta*R)*diff(h(zeta),zeta)/R-2+2*h(zeta))/(zeta*R)^3/pi = 1/8*(2*(zeta*R)*diff(h(zeta),zeta)/R-4*h(zeta)*(zeta*R)*diff(h(zeta),zeta)/R+(diff(h(zeta),`$`(zeta,2))/R^2)*h(zeta)*(zeta*R)^2-2*h(zeta)+2*h(zeta)^2)/(zeta*R)^3/pi/h(zeta);

>

> EcuacHidroAniAdim1:=(diff(h(zeta),`$`(zeta,2))*zeta^2-2*zeta*diff(h(zeta),zeta)-2+2*h(zeta))*h(zeta)+ 2*zeta*diff(h(zeta),zeta)-4*h(zeta)*zeta*diff(h(zeta),zeta)+diff(h(zeta),`$`(zeta,2))*h(zeta)*zeta^2-2*h(zeta)+2*h(zeta)^2=0;

> dsolve({EcuacHidroAniAdim1,h(0)=1, D(h)(0)=0},h(zeta),series);

> solnumerica:= dsolve({EcuacHidroAniAdim1,h(0)=1,D(h)(0)=0},h(zeta),numeric,range=0..1):

Error, (in dsolve/numeric/make_proc) ode system is singular at the initial point

>

Otras estrategias

Ejemplo 1

Ejemplo 2

Ejemplo 3

Ejemplo 4

Ejemplo 5

Ejemplo 6

>

Referencias