Una estrella bien simple

A simple star

Resumo

En este trabajo presentamos el problema de deducir e integrar numéricamente las ecuaciones de equilibrio hidrostático para una estrella con una ecuación de estado politrópica, que es el conocido modelo de Eddington. En particular brindamos los parámetros correspondientes para obtener un modelo del Sol. El tratamiento está orientado a estudiantes de los primeros años de la educación universitaria en física. Se incluye una implementación del algoritmo de integración.

Abstract

In this work we present the problem of obtaining and integrating numerically the hydrostatic equilibrium equations for a star with a polytropic equation of state, which is the well known Eddington model. In particular, we provide the parameters to obtain the corresponding model for the Sun. The treatment is geared towards first or second year physics students. An implementation of the integration algorithm is included.

Keywords:
stellar structure; hydrostatic equilibrium; polytrope

1. Introducción

En la enseñanza universitaria en física solemos buscar ejercicios o problemas que reúnan una serie no trivial de características. Por ejemplo, tienen que ser acordes con el nivel del estudiantado, ser representativos de los temas que se han discutido, y tienen que ayudar a perfeccionar las habilidades que se esperan de una física o de un físico. También tienen que ser suficientemente simples como para ser resueltos en los tiempos acotados de la enseñanza, pero no triviales, de forma de presentar desafíos que motiven a los alumnos y las alumnas. Además es conveniente que una vez resuelto un ejercicio en particular, este se pueda extender o complejizar, de forma de ahondar paulatinamente en los conocimientos relevantes. Como punto extra, está que el problema tenga carácter integrador, esto es, que combine conceptos que en la presentación de la teoría quizás no se ven obviamente conectados. Con estas ideas en mente es que en este artículo proponemos el problema de deducir e integrar las ecuaciones de equilibrio hidrostático con simetría esférica en mecánica newtoniana para un gas. Esto es, buscamos encontrar modelos de estrellas. Este problema en su formulación más básica requiere de fundamentos que incluso en algunos casos se encuentran en la física que se aprende en el secundario, a saber, la segunda ley de Newton y la ley de la gravitación universal. Quizás un poco más avanzado, pero sin pasar de los cursos elementales de física, se encuentra la descripción continua de la materia, en particular de un gas. Como contraparte matemática, se obtiene un sistema de ecuaciones en derivadas ordinarias, que se busca integrar en forma numérica. La parte más complicada del problema es justificar el uso de un politropo como ecuación de estado, sin embargo esto tampoco requiere mucho más allá de la ecuación de los gases ideales. Con estos elementos básicos se puede plantear el problema de encontrar un modelo estelar, el cual presenta conclusiones que no son triviales. Si bien el problema cuenta con una solución analítica, consideramos que es instructivo recurrir a un esquema numérico, que es también una capacidad que se intenta fortalecer dentro del currículo universitario, y que usualmente tiene dificultades para ser integrado en las materias tradicionales, dado el poco tiempo y la cantidad de material que se intenta brindar al alumnado.

Este problema, en su versión más simple, que es la que se presenta aquí, es suficientemente sencillo y corto como para ser parte de una guía de trabajos prácticos y ser resuelto en una semana, y a la vez contiene conceptos importantes y posibilidades de extensión como para formar la base de una discusión más detallada de modelos estelares. El problema plantea varios puntos interesantes, como la deducción de las ecuaciones diferenciales a partir de principios básicos y suposiciones, o la necesidad de una ecuación de estado que no proviene de la mecánica newtoniana. Asimismo se encuentran puntos a discutir respecto a la estructura del sistema de ecuaciones, su regularidad, y sus condiciones iniciales. En el algoritmo numérico también se plantean interrogantes, como el uso de variables adimensionales y la convergencia del método. Finalmente, una vez obtenida la solución numérica nos encontramos que podemos hacer y responder preguntas sobre la estructura estelar.

No es la primera vez que se propone la resolución de las ecuaciones de equilibrio hidrostático en modelos estelares como parte de la formación universitaria. Sin embargo, se han planteado para los últimos años de la carrera, o como parte de algún proyecto de investigación o tesis, como por ejemplo se propone en [15], [14] y [10]. En estas referencias el tratamiento es sustancialmente más avanzado, sin embargo los fundamentos básicos del problema son los mismos, lo cual a nuestro entender muestra la riqueza del problema. Más cercano al nivel de complejidad que hemos elegido se encuentran los trabajos [12] y [16]. Aquí las ecuaciones de equilibrio se resuelven dando el perfil de densidad de la estrella, con lo cual el resto de las variables se obtienen en forma analítica. Si bien es una aproximación absolutamente válida al problema, nos ha parecido superador el hecho de plantear la ecuación de estado como más fundamental que el perfil de densidad. Al imponer el perfil de densidad la ecuación de estado emerge como una consecuencia de la estructura estelar, y no al revés, que es lo que se espera de procesos físicos.

Del lado de la ecuación de estado, si bien los procesos politrópicos son simplificaciones, se han usado como casos especiales de utilidad en el contexto de modelos estelares. Sin duda el ejemplo más exitoso es la obtención del límite de Chandrasekhar para enanas blancas, al respecto recomendamos el tratamiénto didáctico en [5]. Para una discusión sobre un modelo del Sol con ecuaciones de estado politrópicas se puede recurrir a [8]. Algunas de las conclusiones que se desprenden del uso de modelos estelares con politropos relativos a la estabilidad y masas máximas se discuten en forma simplificada, aunque abordando otros temas de importancia, en [9], [1] y [2]. Los modelos de objetos compactos cuya ecuación de estado corresponde a un politropo se han usado también para ejemplificar otros conceptos, como es el tiempo de viaje entre dos puntos de un cuerpo cuando se toma en cuenta la fuerza gravitatoria [4] [6].

Adelantándonos un poco a la discusión de las próximas secciones, nos gustaría aclarar que el problema planteado, en su mínima expresión, consta de las ecuaciones diferenciales (1) y (2), de la ecuación de estado (5), y del algoritmo implementado en el ^{Apéndice A}. La ecuación de estado en particular que vamos a utilizar se puede presentar sin mayor justificación en el planteo del problema si se desea reducirlo a un mínimo. Sin embargo, con dedicarle un poco más de tiempo, y con la ayuda de la ecuación de estado de gases ideales, más una hipótesis sobre la presión de radiación, se puede justificar, y forma una parte interesante de la discusión del modelo. Asimismo, en las siguientes secciones, además de la descripción del problema, hemos incluído preguntas de corte conceptual, dado que el problema puede servir para una rica discusión sobre conceptos que suelen ser difíciles de abordar en los primeros cursos de física. Así, en las secciones que siguen vamos construyendo el problema junto con su solución, motivando los pasos que vamos dando, de forma de que sea posible comprender las motivaciones detrás de las suposiciones y el uso de la teoría.

El artículo está organizado de la siguiente manera. En la Sección 2 se discuten las ecuaciones de equilibrio hidrostático, haciendo especial hincapié en las suposiciones que se utilizan y en las propiedades del sistema de ecuaciones obtenido. En la Sección 3 se presenta y justifica la ecuación de estado. En la Sección 4 se realiza el planteo numérico del problema, la resolución y el análisis de la convergencia. Finalmente en la Sección 5 se presentan las conclusiones, donde también se discuten líneas para continuar profundizando en el estudio de la estructura estelar. Se incluyen dos apéndices. En al ^{Apéndice A} está la implementación del algoritmo numérico, lo que realizamos en SageMath [13], que es un software matemático libre de código abierto. En el ^{Apéndice B} se han dejado las ecuaciones de equilibrio hidrostático relativista y la ecuación de estado para un gas de electrones degenerados, con lo cual, junto con modificar el algoritmo de integración, se pueden obtener modelos de enanas blancas.

2. Ecuaciones de Equilibrio Hidrostático

El objetivo que nos planteamos es encontrar un modelo sencillo de estrella. Entonces antes de comenzar la pregunta fundamental es:

• •

  ¿Qué queremos decir con un “modelo de estrella”?

En general, en la educación secundaria e incluso primaria, cuando se presentan temas de astronomía, en particular cuando se habla de estrellas, se dice “qué son”. O sea, se da una lista de propiedades y características. Por ejemplo, que son esféricas, que están compuestas por gas (mayormente hidrógeno), que “viven” gracias a las reacciones nucleares (mayormente fusión de hidrógeno), que tienen una cierta composición química, temperatura, evolución, etc. Si bien esta descripción cuenta con ventajas pedagógicas obvias en esa etapa de la formación, deja afuera la pregunta de “por qué son” las estrellas. Claramente las estrellas no son objetos fundamentales en física y por lo tanto esperamos que la teoría nos permita explicar y predecir sus características. El hacer esto, el comenzar por principios más fundamentales y predecir cómo es una estrella es lo que llamamos un “modelo de estrella”.

El siguiente paso en nuestro modelo es decidir qué simplificaciones vamos a realizar, lo cual irá en forma de suposiciones. Las estrellas son objetos complicados y un modelo completo incluye muchos ingredientes, entre ellos distribución de densidad, presión, temperatura, composición química, generación y transporte de energía y estado dinámico. El incluir todos estos componentes hacen a un mejor modelo, pero también lo hacen mucho más complicado. En lo que sigue vamos a presentar los pasos para obtener el llamado “modelo de Eddington”, el cual incluye varias simplificaciones fuertes, pero que sin embargo da predicciones suficientemente buenas como para haber sido el modelo utilizado por la comunidad científica durante buena parte del siglo XX.

Comenzamos entonces con dos suposiciones fuertes, que la estrella sea esféricamente simétrica y que sea estática. Consideramos que este es un buen momento para discutir en el curso el valor de las suposiciones en física, con preguntas del tipo:

• •

  ¿Cuáles son las suposiciones que estamos realizando?

• •

  ¿Las suposiciones son realistas?

• •

  ¿Qué libertad nos queda luego de realizar las suposiciones?

• •

  ¿Podemos saber a priori si las suposiciones son razonables?

• •

  ¿Podemos ver a posteriori si las suposiciones fueron razonables?

• •

  ¿Qué pasa si no encontramos ninguna estrella que cumpla con las suposiciones?

• •

  ¿Qué pasa si encontramos una sola estrella que cumpla con las suposiciones?

• •

  ¿Qué pasa si encontramos muchas estrellas que cumplan con las suposiciones?

También es conveniente plantear que “estática” es la forma que usamos para decir que no cambia con el tiempo, y notar que hay una diferencia sutil pero fundamental con que sea “estacionaria”. Asimismo es interesante el concepto de “esféricamente simétrico”, lo que se puede abordar con las siguiente preguntas:

• •

  ¿Cómo decimos que algo es esféricamente simétrico?

• •

  ¿Si no decimos que algo es esféricamente simétrico, entonces podríamos saber que sí lo es?

• •

  ¿Qué pasa si usamos coordenadas que no ayudan, por ejemplo cartesianas o elípticas?

Por supuesto que las preguntas anteriores ayudan en caso de que haya tiempo para una discusión de conceptos, aunque no son estríctamente necesarias para el desarrollo del modelo. Suponemos que las y los estudiantes ya manejan las definiciones de densidad y presión, la segunda ley de Newton y la ley de la gravitación universal, junto con ser capaces de hacer una descripción continua de un gas. Con esto más las suposiciones de estaticidad y simetría esférica las y los estudiantes deberían ser capaces, guíandose con la Figura 1, de deducir las dos ecuaciones de equilibrio hidrostático,

Figura 1
Diagrama para deducir las ecuaciones de equilibrio hidrostático.

(1) $$\frac{dm}{dr}=4\pi r^2\rho ,$$

(2) $$\frac{dp}{dr}=-G\frac{m\rho }{r^2},$$

donde $m\left(r\right)$ es la masa contenida en una esfera de radio $r$, $\rho \left(r\right)$ es la densidad y $p\left(r\right)$ es la presión en el gas, y $G$ es la constante de gravitación universal. La ecuación (1) es la definición de densidad, mientras que (2) es la ley de la gravitación universal, junto con las suposiciones del caso. Justamente el hecho de la convergencia de estos conceptos es lo que consideramos hace especialmente interesante este problema, a saber, densidad, presión, leyes del movimiento de Newton, ley de la gravitación universal y descripción continua de la materia. Y a pesar de incluir todos estos ingredientes, el problema resultante es suficientemente accesible como para ser tratado con un algoritmo numérico sencillo, o a través de un tratamiento analítico si se prefiere seguir ese camino.

En este momento vale la pena detenerse y aprovechar para discutir el carácter del sistema formado por (1) y (2). Se trata de un sistema de ecuaciones diferenciales ordinarias, y es un sistema acoplado, a través de la función $\rho$. Para poder ser resuelto se necesita una ecuación de estado, esto es, una ecuación que relacione $\rho$ y $p$. También es importante notar que dicha ecuación tendría que venir de principios fundamentales o de experimentos respecto de la materia que compone la estrella, pero que dichos experimentos no tienen por qué tener relación con la situación gravitatoria del sistema. Es interesante notar que si uno tiene algún tipo de materia y realiza experimentos en una situación sin gravedad, y encuentra la relación entre $\rho$ y $p$, entonces esa relación también tiene que ser válida cuando se considera la influencia de la gravedad en el sistema. Esta idea de que la ecuación de estado es independiente de la situación gravitatoria del material es fundamental en relatividad general. Aquí también es un buen momento para plantear qué sucede en caso de que no haya una relación que sólo implique $\rho$ y $p$, o sea, si se tiene una ecuación de estado más general, donde por ejemplo $p$ depende de $\rho$ y la temperatura. Volveremos sobre este punto en la siguiente sección.

Habiendo obtenido el sistema de ecuaciones es un buen momento para discutir las propiedades que esperamos tengan las funciones involucradas, entre ellas los ordenes de magnitud esperados. Algunas consideraciones generales son las siguientes:

• •

  Esperamos que haya una región con materia, la estrella, y una región de vacío, el exterior.

• •

  El lugar donde se pasa de una situación a otra define el radio de la estrella, $R$.

• •

  La masa de la estrella está dada por $M=m\left(R\right)$.

• •

  Fuera de la estrella, esto es, para $r>R$, $\rho =0$, lo que también implica $m(r>R)=M$ y $p(r>R)=p\left(R\right)$.

• •

  En el vacío tenemos $p(r>R)=0$, con lo que $p\left(R\right)=0$, y esta va a ser la condición que determina el borde de la estrella.

• •

  En el interior debemos tener $\rho$ positiva y acotada. Entonces $m$ es monótonamente creciente y $p$ es monótonamente decreciente.

Debido a que el sistema de ecuaciones es un sistema de ecuaciones diferenciales ordinarias de primer orden necesitamos como condiciones iniciales el valor de las variables en un punto del intervalo de integración. Dado que hay que integrar el sistema desde $r=0$ hasta $r=R$, y todavía no sabemos cuál es el valor de $R$, es conveniente dar las condiciones iniciales en $r=0$. Esto si bien es conveniente, rápidamente nos damos cuenta que el sistema de ecuaciones es formalmente singular allí. Aquí es un buen momento para discutir al respecto:

• •

  ¿Qué significa que el sistema sea formalmente singular?

• •

  ¿Es el sistema realmente singular?

• •

  ¿Qué produce la singularidad?

• •

  ¿Se puede arreglar la singularidad con un sistema de coordenadas diferente?

• •

  ¿Es más conveniente utilizar otro sistema de coordenadas o buscar una manera de lidiar con la singularidad?

Entonces, como el sistema es formalmente singular, no es posible dar libremente las condiciones de contorno en $r=0$. Por razones físicas podemos ver que $m\left(0\right)=0$, pero también de la ecuación (2) tenemos que para que sea regular necesitamos que $m\left(0\right)=0$, y aún más, que $\frac{dm}{dr}\left(0\right)=0$. De la ecuación (1) tenemos que $\frac{d^{2}m}{d{r}^2}\left(0\right)=0$ y obtenemos la expansión de Taylor de $m$,

(3) $$m\left(r\right)=\frac{4}{3}\pi {\rho }_0{r}^3+O\left(r^4\right),$$

y usando (2),

(4) $$p\left(r\right)=p_0+O\left(r^2\right),$$

donde ${\rho }_0=\rho \left(0\right)$ y $p_0=p\left(0\right)$. Cuando tengamos la ecuación de estado veremos que el único parámetro libre es $p_0$. Por ahora nos quedamos con esta condición de contorno, que es la que vamos a utilizar para la integración, y en la Sección 5 discutiremos brevemente otras opciones.

3. Ecuación de Estado

En general, para construir un modelo estelar necesitaríamos ecuaciones para la generación y transporte de energía, para la variación de la temperatura y para la composición química dentro de la estrella. Esto complicaría enórmemente el problema. Sin embargo, de (1) y (2) vemos que si tuviéramos una ecuación de estado que relacione $\rho$ y $p$, entonces dichas ecuaciones podrían ser integradas. La forma de interpretar esto es que el resto de las variables de las que dependería la estrella se pueden deducir a posteriori habiendo integrado (1) y (2), a través de relaciones entre esas otras variables y $\rho$ y $p$. Vamos a utilizar un politropo con exponente $\frac{4}{3}$,

(5) $$p=K{\rho }^{\frac{4}{3}},$$

donde $K$ es una constante que consideramos puede ser diferente para distintas estrellas. En una primera aproximación al problema se puede dar esta ecuación de estado sin mayor justificación, sin embargo, como haremos a continuación, se puede justificar utilizando una suposición adicional. Así nos encontramos con un ejemplo interesante en sí mismo de cómo una condición extra puede fijar la relación funcional de la temperatura en la ecuación de estado. Para esto vamos a considerar que la materia está dada por un gas ideal, con ecuación de estado

(6) $$p_{gas}=\frac{R}{\mu }\rho T,$$

que es una ecuación a la vez útil y simple, que suponemos que ya se ha enseñado en la carrera. Aquí, $R$ es la constante universal de los gases y $\mu$ es la masa molar del gas en cuestión. Nuevamente es un buen momento para discutir si la suposición es válida o no, cuál será su rango de validez y qué tanto podemos confiar en las conclusiones que obtengamos. También se puede discutir qué sucede si simplemente reemplazamos (6) en (1) y (2) e intentamos integrarlas. Como camino alternativo podemos pensar en una estrella isotérmica, esto es, una estrella que tiene en todo su interior la misma temperatura. Sin embargo en la integración esto presenta otros inconvenientes que preferimos evitar. En su lugar vamos a considerar que también existe presión de radiación, cuya ecuación es

$$p_{rad}=\frac{4\sigma }{3c}{T}^4,$$

donde $\sigma$ es la constante de Stefan-Boltzmann y $c$ es la velocidad de la luz en el vacío. La presión total está dada por

$$p=p_{gas}+p_{rad}.$$

La condición extra que vamos a incluir es que la proporción de $p_{gas}$ y $p_{rad}$ permanece constante en todo el interior, para lo que definimos una constante, $\beta$, a través de

$$\beta =\frac{p_{gas}}{p}.$$

Eliminando $T$ obtenemos

$$p\left(\rho \right)={\left(\frac{3c{R}^4}{4\sigma {\mu }^4}\right)}^{\frac{1}{3}}{\left(\frac{1-\beta }{{\beta }^4}\right)}^{\frac{1}{3}}{\rho }^{\frac{4}{3}}.$$

Esta expresión es la misma que (5) y como $\beta$ es un parámetro en este modelo $K$ no está fijada por constantes fundamentales.

4. Integración Numérica

Finalmente, el sistema de ecuaciones es

$$\frac{dm}{dr}=\frac{4\pi }{K^{\frac{3}{4}}}{r}^2{p}^{\frac{3}{4}},$$

$$\frac{dp}{dr}=-\frac{G}{K^{\frac{3}{4}}}\frac{m{p}^{\frac{3}{4}}}{r^2},$$

con las condiciones iniciales

$$m\left(0\right)=0,p\left(0\right)=p_0.$$

Para integrarlo vamos a utilizar el método de Euler, que es sencillo y está directamente relacionado al concepto de derivada, por lo que permite también que las alumnas y los alumnos afiancen estos temas. Vamos a recapitular brevemente el método de Euler. Supongamos que tenemos una ecuación diferencial de la forma

$$\frac{dy}{dx}=f(x,y),$$

entonces aproximamos la derivada por la taza de cambio,

$$\frac{dy}{dx}\approx \frac{\Delta y}{\Delta x},$$

donde

$$\Delta x=x_{i+1}-x_i,\Delta y=y\left(x_{i+1}\right)-y\left(x_i\right),$$

siendo $x_i$ el punto $i$-ésimo de la grilla de integración. Despejando obtenemos

$$y\left(x_{i+1}\right)=y\left(x_i\right)+f(x_i,y\left(x_i\right))\Delta x.$$

El algoritmo de integración es directo, dado $y\left(x_0\right)$ cada paso de la integración calcula $y\left(x_{i+1}\right)$, hasta llegar a una condición que haga que se detenga la integración, que puede ser que se haya recorrido la grilla en la cual se quiere integrar, o bien en nuestro caso que una de las variables dependientes sea cero. Con esta breve descripción el estudiantado debería estar en condiciones de entender el algoritmo para el caso que estamos considerando.

Para la primera integración sugerimos utilizar parámetros correspondientes al Sol, los cuales han sido ajustados por fuera de lo que hemos presentado, dados por

$${\rho }_{\odot 0}\approx 8\times {10}^4\frac{kg}{m^3},p_{\odot 0}\approx 1,3\times {10}^{16}Pa,$$

con lo que se obtiene

$$K_{\odot }=3,8\times {10}^9\frac{Pa}{{(kg/m^3)}^{\frac{4}{3}}}.$$

Para una comparación gruesa pero instructiva con el Sol, tenemos que

$$M_{\odot }\approx 2\times {10}^{30}kg,R_{\odot }\approx 7\times {10}^{5}km.$$

Con estos parámetros se puede ir al integrador y obtener un modelo del Sol. Sin embargo, antes de recurrir al algoritmo sugerimos un paso más. Es importante que quede claro que al realizar la integración numérica se incurre en errores numéricos. En particular, si se tiene que trabajar con cantidades cuyos valores numéricos difieren en muchos órdenes de magnitud entonces se potencia la posibilidad de que se acumulen errores numéricos. Por esto es conveniente trabajar dentro de lo posible con cantidades adimensionales de órden de magnitud cercano a la unidad. Como estamos buscando un modelo de estrella para el Sol vamos a usar los parámetros que estimamos del mismo como unidad de medida para las cantidades involucradas. Por esto definimos

$$\widehat{r}=\frac{r}{R_{\odot }},\widehat{m}=\frac{m}{M_{\odot }},\widehat{\rho }=\frac{\rho }{{\rho }_{\odot 0}},\widehat{p}=\frac{p}{p_{\odot 0}},$$

con lo que las ecuaciones toman la forma

(7) $$\frac{d\widehat{m}}{d\widehat{r}}=A{\widehat{r}}^2\widehat{p},A=\frac{4\pi R_{\odot }^3{p}_{\odot 0}^{\frac{3}{4}}}{M_{\odot }{K}_{\odot }^{\frac{3}{4}}},$$

(8) $$\frac{d\widehat{p}}{d\widehat{r}}=-B\frac{\widehat{m}{\widehat{p}}^{\frac{3}{4}}}{{\widehat{r}}^2},B=\frac{G{M}_{\odot }}{R_{\odot }{K}_{\odot }^{\frac{3}{4}}{p}_{\odot 0}^{\frac{1}{4}}}.$$

Queda solo un detalle por discutir, y es que como el sistema es formalmente singular en $r=0$ el método de Euler no va a funcionar si intentamos integrar desde ese punto. En lugar de esto vamos a trasladar las condiciones iniciales al primer punto de la grilla. O sea, si $\Delta \widehat{r}$ es el espaciado en la grilla, y las condiciones iniciales son $\widehat{m}\left(0\right)=0$ y $\widehat{p}\left(0\right)={\widehat{p}}_0$, entonces vamos a comenzar a integrar en $\widehat{r}=\Delta \widehat{r}$ y las condiciones iniciales son

$$\widehat{m}\left(\Delta \widehat{r}\right)=0,\widehat{p}\left(\Delta \widehat{r}\right)={\widehat{p}}_0.$$

Aunque esto pueda parecer extremadamente burdo, está justificado por las expansiones (3) y (4), y en caso de requerirse mayor precisión se pueden utilizar órdenes superiores de la expansión de Taylor para dar condiciones iniciales más precisas en $\widehat{r}=\Delta \widehat{r}$.

Con esto cerramos el trabajo de preparación y en el ^{Apéndice A} presentamos una implementación del algoritmo en SageMath. El algoritmo es suficientemente sencillo como para ser implementado casi sin experiencia previa y puede ser implementado en múltiples lenguajes de programación, e incluso en una hoja de cálculo, como se ejemplifica para enanas blancas en [10]. Si bien SageMath cuenta con integradores más precisos ya implementados, e incluso también con una función implementada del método de Euler, preferimos hacer un algoritmo explícito, dado el carácter pedagógico del modelo. El programa está comentado y consideramos que es suficientemente transparente como para que con una mínima guía quienes aprenden lo entiendan y sean capaces de modificarlo a voluntad. En la Figura 2 se muestra el resultado de la integración por este método obtenida con el algoritmo presentado. En una primera versión del problema nos podríamos dar por satisfechos con llegar hasta aquí.

Figura 2
Resultados de la integración numérica.

Si se quiere aprovechar el trabajo realizado para explorar las posibilidades del modelo, el siguiente paso puede ser la construcción de una familia de soluciones. Para esto lo más directo es variar la presión central. En la Figura 3 se muestran cinco soluciones con distintas presiones centrales.

Figura 3
Soluciones con diferente ${\widehat{p}}_0$.

Una vez que se tiene un esquema numérico es importante saber que el mismo converge a la solución. En el caso que estamos tratando existe la solución analítica a través de las funciones de Lane-Emden, lo que permitiría saber exactamente cuánto error se comete en la solución numérica. Esto sería especialmente útil si luego el modelo se va a usar para situaciones donde no se tiene una solución analítica, de forma que este caso particular serviría como testeo del funcionamiento correcto del algoritmo. A pesar de esto, consideramos que en esta primera aproximación al problema, y como quizás también es para el estudiantado una primera aproximación a una solución numérica, es conveniente testear el método sin recurrir a la solución analítica. Para esto lo que hacemos es resolver el problema con diferentes grillas, y comparamos las soluciones con grillas más espaciadas con respecto a la solución con la grilla menos espaciada. En la Figura 4 se presentan los errores con respecto a la solución con $\Delta \widehat{r}={10}^{-6}$. El eje vertical está en escala logarítmica, lo que muestra que la relación entre el espaciado y el error sigue una ley de potencia, lo que es típico de los métodos numéricos convergentes, y que nos da confianza de que la solución es correcta y precisa, además de darnos una idea de la magnitud del error.

Figura 4
Errores en las soluciones con ${\widehat{p}}_0=1$ para distintas grillas con respecto a la solución con $\Delta \widehat{r}={10}^{-6}$.

5. Conclusiones

En las secciones anteriores hemos presentado el problema de encontrar un modelo estelar para el Sol, en una versión suficientemente simple como para que sea un problema que requiere conocimientos elementales de física y prácticamente ninguna experiencia numérica. Por lo tanto consideramos que es un buen primer paso en la resolución de problemas no triviales a través de métodos numéricos. El problema está planteado de forma que puede ser parte de una guía de ejercicios y que en su versión más elemental no requeriría más de una semana para ser resuelto, y pensamos que es un buen problema para integrar conceptos, empezar a tomar habilidad con los métodos numéricos, y profundizar en el entendimiento de cómo los métodos analíticos, la teoría y los métodos numéricos se complementan.

Por supuesto que un problema así puede ser el punto de partida para discusiones más elaboradas. Por el lado analítico, se puede continuar con el entendimiento de los politropos, la ecuación de Lane-Emden, sus propiedades y las propiedades de sus soluciones. Por el lado numérico, se pueden realizar mejoras sustanciales al algoritmo, comenzando por analizar la precisión con respecto al tamaño de la grilla y el tiempo de cómputo. También una adición sencilla pero interesante es implementar la serie de Taylor para las condiciones iniciales, siempre controlando cómo el cambio afecta la precisión del método. En esta misma línea se puede modificar la grilla de forma de tener más puntos cerca del borde singular. Esto se puede hacer poniendo una sección más densa cerca del borde singular y otra menos densa en el otro extremo, o haciendo que el espaciado de la grilla no sea uniforme. Se puede pasar a métodos numéricos más sofisticados, siendo el candidato natural el método de Runge-Kutta-Fehlberg, el cual se encuentra implementado en SageMath así como en otros softwares matemáticos. Finalmente, consideramos que otra forma interesante de modificar el algoritmo, y una que es utilizada en modelos estelares, es el uso de métodos de shooting. Esto se puede implementar para por ejemplo encontrar estrellas con una masa y radio fijo, o para ajustar el valor de $K$ en el modelo una vez que se han determinado otros parámetros. Por el lado de la física, también hay varios caminos a seguir. Uno es el pasar de las ecuaciones newtonianas a las ecuaciones relativistas, las ecuaciones de Tolman-Oppenheimer-Volkoff, que se presentan en el ^{Apéndice B}. Es interesante que una vez que se tiene el modelo newtoniano, el paso al modelo relativista no es particularmente difícil, mientras que deducir las ecuaciones de Tolman-Oppenheimer-Volkoff ciertamente requiere conocimientos más avanzados. Otro camino interesante es mejorar la ecuación de estado, por ejemplo para un gas de Fermi, que aplica a enanas blancas o a la interacción entre nucleones, lo cual implica una profundización importante en otras áreas por fuera de la mecánica newtoniana. Finalmente, se puede seguir el camino de hacer el modelo más completo, incluyendo diferentes elementos químicos, termodinámica, tazas de reacción y generación de energía, consultando por ejemplo los libros [3], [11] y [7]. Es llamativo que empezando con un problema humilde como el que se plantea aquí se puede llegar a temas de investigación de actualidad, pero de alguna manera esa ha sido la carta ganadora de la física desde que cristalizó como ciencia.

Reconocimientos

La cálculos numéricos y las figuras fueron realizadas en SageMath [13].

Este trabajo se realiza con el apoyo de los subsidios PIP 112-201301-00532 de CONICET, Argentina, y M060 de SIIP, Universidad Nacional de Cuyo, Argentina.

A. Implementación del algoritmo en SageMath

En este apéndice presentamos una implementación del método de Euler en SageMath para las ecuaciones (7) y (8) con los parámetros estimados para el Sol. Lo que sigue se puede copiar y pegar en una notebook de SageMath y se obtienen los gráficos de la Figura 2.

# Constantes utilizadas en la integracion.

G = 6.674*10^(-11)

Rs = 7*10^8

Ms = 2*10^30

rhos = 8*10^4

ps = 1.3*10^16 #

Constantes calculadas de las constantes anteriores.

Ks = float(ps/(rhos^(4/3))) ; print(Ks)

A = float(4*pi*Rs^3*ps^(3/4)/(Ms*Ks^(3/4))); print(A)

B = float(G*Ms/(Ks^(3/4)*ps^(1/4)*Rs)); print(B)

# Condiciones iniciales para la integracion. m0 = 0 p0 = 1

# Paso en la integracion. delta = 0.001

# Punto inicial en la grilla. r0 = delta

# Listas en las cuales se guarda la solucion.

m = [[r0,m0]]

p = [[r0,p0]]

# Integracion, utilizando la condicion de que

# la presion es positiva en el interior. while p0 > 0:

m0 = m0 + A * r0^2 * p0^(3/4)

* delta

p0 = p0 - B * m0 * p0^(3/4) / r0^2

* delta

r0 = r0 + delta

m.append([r0,m0])

p.append([r0,p0])

# Graficos de la solucion.

figM = list_plot(m,color='blue')

figP = list_plot(p,color='red')

show(figM + figP)

B. Ecuaciones relativistas para enanas blancas

En el caso de objetos compactos con altas densidades, el modelo de gas ideal deja de ser válido. El ejemplo más directo de esto son las enanas blancas. En este caso hay que pasar a la descripción relativista, y usar una ecuación de estado que tenga en cuenta el estado de la materia.

Las ecuaciones de equilibrio hidrostático para una configuración esféricamente simétrica en relatividad general son las llamadas ecuaciones de Tolman-Oppenheimer-Volkoff,

$$\frac{dm}{dr}=4\pi r^2\rho ,$$

$$\frac{dp}{dr}=-G\frac{m\rho }{r^2}\left(1+\frac{p}{c^2\rho }\right)\left(1+\frac{4\pi r^{3}p}{m{c}^2}\right){\left(1-\frac{2Gm}{c^{2}r}\right)}^{-1}.$$

La ecuación de estado para un gas de electrones degenerados en el límite relativistas es

$$p=\frac{\hslash c}{12{\pi }^2}{\left(\frac{3{\pi }^{2}Z}{m_{N}A}\right)}^{\frac{4}{3}}{\rho }^{\frac{4}{3}},$$

donde $A/Z$ es el número de nucleones por electrón y $m_N$ es la masa de cada nucleón.

Referências

• [1] M. Bandecchi, J. Horvath y P. Bretones, Rev. Bras. Ens. Fís. 41, e20180250 (2019).
• [2] M. Bandecchi, P. Bretones y J. Horvath, Rev. Bras. Ens. Fís. 41, e20190031 (2019).
• [3] S. Chandrasekhar, An Introduction to the Study of Stellar Structure (Dover, New York, 2010).
• [4] W. Dean Pesnell, Am. J. Phys. 84, 192 (2016).
• [5] E. Evangelista, Rev. Bras. Ens. Fís. 41, e20180167 (2019).
• [6] A. Gjerløv y W. Dean Pesnell, Am. J. Phys. 87, 452 (2019).
• [7] C. Hansen, S. Kawaler y V. Trimble, Stellar Interiors (Springer, New York, 2004), 2nd ed.
• [8] A. Hendry, Am. J. Phys. 61, 906 (1993).
• [9] F. Herrmann y H. Hauptmann, Am. J. Phys. 65, 292 (1997).
• [10] C. Jackson, J. Taruna, S. Pouliot, B. Ellison, D. Lee y J. Piekarewicz, Eur. J. Phys. 26, 695 (2005).
• [11] R. Kippenhahndolf, A. Weigert y A. Weiss, Stellar Structure and Evolution (Springer, Heidelberg, 2012), 2nd ed.
• [12] H. Rodrigues, Eur. J. Phys. 34, 667 (2013).
• [13] SageMath, the Sage Mathematics Software System (Version 9.1), accessed in 12/08/2020: https://www.sagemath.org
  » https://www.sagemath.org
• [14] I. Sagert, M. Hempel, C. Greiner y J. Schaffner-Bielich, Eur. J. Phys. 27, 577 (2006).
• [15] R. Silbar y S. Reddy, Am. J. Phys. 72, 892 (2004).
• [16] A. Smirnov y R. Oliveira, Rev. Bras. Ens. Fís. 37, 1308 (2015).

Fechas de Publicación

• Publicación en esta colección
  06 Nov 2020
• Fecha del número
  2020

Histórico

• Recibido
  12 Ago 2020
• rev-received
  15 Set 2020
• Acepto
  22 Set 2020