Numerical solution of the TOV equations for CGF stars

This Sagemath notebook does calculations for the paper Dark matter stars.

Preamble

Code to set floating point precision

mpf(x) converts sage floating point numbers to mp numbers
Reals(x) or myR(x) converts mp numbers to sage

setting mp.dps tells mpmath to use decimal precision mp.dps or a small amount more
Sagemath decimal precision is set to mpmath decimal precision + 3

Numerical parameters

Copied from A theory of the dark matter

CGF scale

Code

The TOV equations form an ode with independent variable $\hat r$ and dependent variables $\hat p$, $\hat m$, $\hat e$.

There are two odes:

TOV equations

Prepare the core ode

change variables

The density $\hat \rho$ decreases monotonically with the radius $\hat r$ so we can take $\hat \rho$ as the independent variable in the ode.
This allows integrating the ode from the initial condition $\hat \rho(0)$ down to $\hat \rho = \hat \rho_{EW}$ which is the density at the outer surface of the core.
The ode solver wants the independent variable to increase, so we use $$ u = \frac{\hat \rho_{EW}}{\hat \rho} \qquad u_0 \le u \le 1 $$ as independent variable. To simplify the ode, we change dependent variables to $$ w = \frac{\hat m}{\hat\rho_{EW}\hat r^3} \qquad x=\hat r^2 \qquad z = \frac{\hat e}{\hat\rho_{EW}} $$

make ode derivatives into functions callable by the ode solver

$\frac{dw}{du}$, $\frac{dx}{du}$, $\frac{dz}{du}$

prepare the initial condition

The ode is singular at $\hat r=0$ so start the ode at $u= u_0(1+\epsilon)$.

define a function to convert ode-solver output to star data

The ode-solver outputs y=[w,x,z].
The star data is ($\hat r$, $\hat m$, $\hat e$, $\hat \rho$, $\hat p$)

(rhohat_norm_0_actual, core_ode) = make_core_ode(rhohat_norm_0)

rhohat_norm_0 $ = \hat \rho_{norm}(0)$ is the initial condition at $\hat r=0$.
$$ \hat \rho_{norm} = {\frac{\hat \rho}{\hat \rho_{EW}}} = \frac1{u} $$

rhohat_norm_0_actual = $(1 - \epsilon)\hat \rho_{norm}(0)$ is the actual initial value given to the ode-solver.

core_ode $\colon \hat \rho_{norm} \mapsto (\hat r,\hat m,\hat e,\hat \rho, \hat p)$ is the solution of the ode for $$ \hat \rho_{norm} \le (1 - \epsilon)\hat \rho_{norm}(0) $$

The actual ode-solver is   [w,x,z] = the_odefun(u) for $$ (1+\epsilon) u_0 \le u $$

Prepare the shell ode

The low density equation of state is expressed implicitly $$ \hat \rho, \hat p = \hat \rho(k^2), \hat p(k^2) \qquad \frac12 \ge k^2 \ge0 \qquad \hat \rho(1/2) = \hat\rho_{EW} \qquad \hat \rho(0) = 0 $$ The functions $\hat \rho(k^2)$, $\hat p(k^2)$ are expressed by a series of formulas listed in Table 2 of A theory of the dark matter.
The formulas involve the the ratio $K/E$ where $K(k)$ and $E(k)$ are the complete elliptic integrals. $k^2$ is the elliptic parameter.
The usual notation is $m=k^2$.

The TOV equations are solved with $$ \text{independent variable}\quad \sigma = 1-2m \qquad \text{dependent variables}\quad \hat r,\,\hat m,\,\hat e $$ During the preparation of the ode it is convenient to use $$ s = -\ln m \qquad \frac{d}{ds} = - m \frac{d}{dm} = \frac{1-\sigma}{2}\frac{d}{d\sigma} $$

code to define functions of $m$

The following code enters the Table 2 formulas as SageMath symbolic expressions.

Next, substitutions are made in the formulas to define functions $$ m, \, E/K \mapsto \hat\rho, \,\hat p,\, \frac{d \hat p}{d\sigma} $$ For efficiency, the elliptic integrals $E,\,K$ will be calculated separately.

define functions $$ m \mapsto (\hat \rho,\hat p, d\hat p/ds) \qquad m \mapsto (\hat \rho,\hat p) $$

define the function that provides the derivatives of the dependent variables to the ode solver. $$ \sigma, y \mapsto \left[ \frac{d\hat r}{d\sigma}, \,\frac{d\hat m}{d\sigma},\,\frac{d\hat e}{d\sigma} \right ] $$

How stars are represented

A star is generated from an initial condition $\hat \rho_{norm}(0)$.

For $\hat \rho_{norm}(0)>1$, first the core_ode is constructed and solved to the endpoint $u=1$ to get the core radius, mass, and energy. These are then used as initial conditions for the shell_ode which is solved to the endpoint $\sigma=1$ to get the star radius, mass, and energy.

For $\hat \rho_{norm}(0)\le 1$, only the shell_ode is constructed and run.

The results are stored in a star which is a dictionary with two keys, 'data' and 'ode'.

Numerical values are mpf floating point numbers. 'core profile' and 'shell profile' are not used here. Their values would be the intermediate values of the TOV solution for the star. Here we only calculate at the ode endpoints: the data at the boundary of the core and at the boundary of the star.

How stars are generated and stored

The object is to generate a set of stars more or less evenly spaced on the star curve in the $M$, $R$ plane.

The even space is more or less accomplished by an algorithm.

The data is collected in subsets to allow for twiddling the decimal precision and the spacing parameters.

A subset of stars is called a star_set. This is a dictionary of stars. The dictionary key is the star's initial value converted to a string. The value is a star.

The star_sets are stored in two lists, one for the stars with cores, the other for the coreless stars.

code to make a store with core

star = make_star( $\hat\rho_{norm}(0)$)

code to adaptively sample parameter space

The last 3 points are used to calculate the radius of curvature $A/B$ (radius of the circle passing through the points).

At the beginning of a run, 3 stars have to be generated to initialize the adaptive algorithm.

The spacing in the $\hat M$-$\hat R$ plane to the next point is $$ \frac{A}{A+c B} \epsilon' $$ If the radius of curvature is large, the spacing is $\epsilon'$. If the radius of curvature is small, the spacing is $\frac{A}{cB} \epsilon'$.

Then linear extrapolation gives the parameter step that will give this spacing.

code for coreless stars

For $\hat \rho_{norm}(0) \le 1$ the star has no core. Only the shell ode is run.

Again $\hat r =0$ is a singular point of the ode. The initial condition is again shifted by $\epsilon$.

$$ \begin{gathered} \sigma = \sigma_0+\epsilon \\ c_1 = \hat r \frac{d\hat r}{d\sigma}_{/\sigma_0} \qquad \hat r^2 = 2 c_1 \epsilon \\ \frac{d\hat p}{d \hat r} = -(\hat \rho_0+\hat p_0) \left( \frac13 \hat\rho_0 +\hat p_0\right) \hat r \qquad \hat p = \hat p_0 - \hat \rho_0+\hat p_0) \left( \frac13 \hat\rho_0 +\hat p_0\right) \frac12 \hat r^2 \\ \hat m = \hat e = \frac13 \hat \rho_0 \hat r^3 \end{gathered} $$

code to initialize the adaptive algorithm from the last 3 stars in a star_set

code for printing and plotting

Plot equation of state

generate stars with cores

generate coreless stars

The star at $\rho(0)=\rho_{EW}$

The star at $\hat\rho(0)= 10^8$

Plots

collect star_data in two dictionaries, cored_stars and coreless_stars

stars is the union of the two, a dictionary of all the stars

The star's key in these dictionaries is central_density (number, not string), for sorting.

Special stars

The star of minumum radius.

3.20219250e+00 r=1.204027e+00 m=1.883777e-01 e=2.477540e-01 BE=5.937627e-02 BE/M=3.151979e-01 t=6

The star of largest mass

3.47709480e-01 r=1.896051e+00 m=3.108239e-01 e=3.670901e-01 BE=5.626628e-02 BE/M=1.810230e-01 t=12
3.41803251e-01 r=1.918823e+00 m=3.109591e-01 e=3.661058e-01 BE=5.514674e-02 BE/M=1.773441e-01 t=12
3.35874144e-01 r=1.941588e+00 m=3.108993e-01 e=3.648964e-01 BE=5.399716e-02 BE/M=1.736805e-01 t=12

The mass where transitions become possible

9.23643135e+00 r=1.271494e+00 m=1.733260e-01 e=2.278757e-01 BE=5.454964e-02 BE/M=3.147227e-01 t=6
1.02705121e+01 r=1.283822e+00 m=1.731541e-01 e=2.274269e-01 BE=5.427287e-02 BE/M=3.134369e-01 t=6
1.14397938e+01 r=1.296846e+00 m=1.732126e-01 e=2.272499e-01 BE=5.403726e-02 BE/M=3.119707e-01 t=6

The star of largest $\mathrm{BE}/M$

Not used in the paper

radius / Schwarzschild radius

Core radius / radius and core mass / mass

mass profile of a cloud

mass profiles of massive stars

3 stars of same mass, different radius and binding energy

1.19892328e-01 r=2.721241e+00 m=1.812195e-01 e=1.913617e-01 BE=1.014224e-02 BE/M=5.596662e-02 t=12
2.86784664e+01 r=1.403815e+00 m=1.812014e-01 e=2.350428e-01 BE=5.384138e-02 BE/M=2.971356e-01 t=6
4.47166984e+00 r=1.211533e+00 m=1.809907e-01 e=2.385063e-01 BE=5.751560e-02 BE/M=3.177822e-01 t=6