Numerical Construction of Neutron Star Models in GR

Specify Parameters of the integration.

K = 5.38 * 10^9;
Γ = 5/3;
(Subscript[ρ, 0  c]) = 10^15;
G  =  6.67*10^(-8);
c = 3*10^(10);
(Subscript[ρ, 0  c]) = (Subscript[ρ, 0  c])*
          G/c^2*(100*10^3)^2;
K  =  
        K*(G c)^(-2/3)*(100*10^3)^(-4/
                3);
\(\(Print[\*""(Subscript[ρ, 0  c]) = "", ρ\_\(0  c\), \*"" (km^(-2))""];\)\)
\(\(Print["K = ", K, \*"" \!\(km\^\(4/3\)\)""];\)\)

Specify TOV Equations.

ρ[
          r_]  =  ρ0[
            r] + ((K (((ρ0[r]))^Γ))/(Γ - 1));
P[r_]  =  K (((ρ0[r]))^Γ);
MassEq  =  (m'[r]  ==  
            4  π ρ[r] (r^2));
PressureEq  =  (P'[
              r]  == -(((ρ[r] + P[r]) (m[r] + 
                        4  π (r^3) P[r]))/(r (r - 
                        2 m[r]))));

Integrate and find \(TraditionalForm\`\(M\_*\)\ and\ \(R\_*\)\).

Choose small r to avoid singularity at 0 (or could linearize the system near 0...).

rSmall  =  0.1;
rEnd = 30;

Integrate. Find \(TraditionalForm\`\(R\_*\)\ and\ \(M\_*\)\ and\ \(\(ρ\_c\)\(.\)\)\)

state = 
        First[NDSolve`ProcessEquations[{MassEq, PressureEq, 
              m[rSmall] == 
                4  π ((Subscript[ρ, 0  c]) + ((K ((((Subscript[ρ, 0  c])))^Γ))/(Γ - 1))) (rSmall^3)/3 , ρ[
                  rSmall]  ==  (Subscript[ρ, 0  c]) + ((K ((((Subscript[ρ, 0  c])))^Γ))/(Γ - 1))}, {m, ρ0}, {r, 
              rSmall, rEnd}]];
Off[NDSolve::"mxst"];
NDSolve`Iterate[state,  rEnd]
solution  =  
        NDSolve`ProcessSolutions[state];
rStar  =  
        Last[state@"CurrentTime"[]];
mStar  =  
        Re[m[r] /. solution /. 
            r -> rStar];
\(\(Print[\*""(Subscript[ρ, c]) = "", ρ\_\(0  c\) + \(K \((ρ\_\(0  c\))\)\^Γ\)\/\(Γ - 1\), \*"" (km^(-2))""];\)\)
\(\(Print[\*""\!\(\(R\_*\)\) = "", 
        rStar, " km"];\)\)
\(\(Print[\*""\!\(\(M\_*\)\) = "", 
        mStar, \ " km"];\)\)
On[NDSolve::"mxst"];

Integrate proper mass; find Δ.

mProper = 
        Re[NIntegrate[
            4 π  ((ρ0[
                        r] + ((K (((ρ0[r]))^Γ))/(Γ - 1))) /. 
                  solution) (r^2) (((1 - ((2 (m[r] /. solution))/r)))^(-1/2)), {r, rSmall, 
              rStar}]];
Delta  =  ((mProper - mStar)/mStar);
\(\(Print[\*""\!\(M\_p\) = "", 
        mProper, \ " km"];\)\)
Print["Δ = ", Delta,  " or ", 
        Delta*100, "%."];

Plots.

Off[Plot::"plnr"];
Off[InterpolatingFunction::"dmval"];
Plot[Evaluate[m[r] /. solution], {r, rSmall, 6.  rStar}, 
        AxesLabel -> {"r (km)", "mass (km)"}];
\(\(Plot[Evaluate[ρ0[r] /. solution], {r, rSmall, 6.  rStar}, 
        AxesLabel -> {"r (km)", \*""(Subscript[ρ, 0]) (km)""}];\)\)
On[Plot::"plnr"];
On[InterpolatingFunction::"dmval"];