# Numerically solving the TOV equations: how to find the max radius and mass?

I'm solving the TOV (Tolman-Oppenheimer-Volkoff) equations using Mathematica 13.2, and I'm having some issues at finding the outer radius $$r_{\text{max}}$$ and mass $$M_{\text{max}}$$ of the spherical star. There's something fishy related to the physics and not to Mathematica, that I don't understand. That question have been asked here a few times, but the answers aren't satisfying to me. So here it is:

The equilibrium equation for a spherical mass is well known. From General Relativity: $$$$\tag{1} \frac{dp}{dr} = - G \frac{(M_r + 4 \pi p r^3)(\rho + p)}{r (r - 2 G M_r)},$$$$ where $$M_r(r)$$ is the mass enclosed inside a sphere of radius $$r$$: $$$$\tag{2} M_r(r) = \int_0^r \rho(r) 4 \pi r^2 dr \qquad \Rightarrow \qquad \frac{dM_r}{dr} = \rho 4 \pi r^2,$$$$ $$\rho$$ is the energy density (not the mass density!). Since there are three unknown: $$p(r)$$, $$\rho(r)$$ and $$M_r(r)$$, we need a third equation to close the system (the state equation). For a polytropic fluid: $$$$\tag{3} p(\rho_{\text{mass}}) = \kappa \rho_{\text{mass}}^{\gamma},$$$$ where $$\gamma > 1$$. Since $$\rho = \rho_{\text{mass}} + p/(\gamma - 1)$$, we may write the third equation as this: $$$$\tag{4} \rho = \Bigl( \frac{p}{\kappa} \Bigr)^{1/\gamma} + \frac{p}{\gamma - 1}.$$$$ A simple change of scale let us integrate the three equations (1), (2) and (4) with Mathematica: $$\tilde{p} = \alpha p$$, $$\tilde{\rho} = \alpha \rho$$, $$\tilde{r} = \beta r$$, $$\tilde{M}_r = \alpha \beta^3 M_r$$. The arbitrary constants $$\alpha$$ and $$\beta$$ are selected such that $$\tilde{G} = G/\alpha \beta^2 \equiv 1$$ and $$\kappa$$ disappears from equation (4). Thus: \begin{align}\tag{5} \tilde{G} &\equiv 1, &\tilde{\kappa} &\equiv 1. \end{align} The initial conditions are $$\tilde{p}(0) = q$$ (the central pressure, where $$q$$ is an arbitrary number) and $$\tilde{M}_r(0) = 0$$. All this is well known and I don't have any issue yet.

The integration with Mathematica works very well and I can plot the pressure $$\tilde{p}$$ as a function of $$\tilde{r}$$, for any values of $$q$$ (the central pressure) and $$\gamma$$ (the adiabatic index). But then, I need to find the outer radius $$R$$ and total mass $$M_{\text{tot}}$$ of the star. Usually, we use $$p(R) = 0$$ and $$M_{\text{tot}} = M_r(R)$$. This is where things get fishy. Here's a typical plot of the internal pressure $$\tilde{p}$$ as a function of $$\tilde{r}$$ (unitless variables). The blue curve is the relativistic solution of the TOV equations, in the case $$\gamma = \frac{4}{3}$$, while the dashed curve is the solution from newton's theory:

Plot of $$\mathrm{Log}_{10} \tilde{p}$$, for the same $$\gamma = \frac{4}{3}$$:

Apparently, the pressure never get to 0! Unless I made a mistake with my Mathematica programming (I have strong confidence there isn't any), I can't find the outer radius with $$\tilde{p} = 0$$! How can I find $$R = r_{\text{max}}$$?? Should I use some pressure cutoff at an arbitrary value? Using a cutoff like $$\tilde{p}_{\min} = 0.0001$$ is very arbitrary and can give any radius, so what's going on here?

EDIT: About the Mathematica coding of this numerical integration, you can see some of the details here:

https://mathematica.stackexchange.com/questions/288315/how-to-remove-an-apparent-singularity-in-ndsolve

EDIT 2: Here's an example of the Mass-Radius relation curve:

• Comments have been moved to chat; please do not continue the discussion here. Before posting a comment below this one, please review the purposes of comments. Comments that do not request clarification or suggest improvements usually belong as an answer, on Physics Meta, or in Physics Chat. Comments continuing discussion may be removed.
– Buzz
Commented Aug 6, 2023 at 1:05

The TOV equations in their original form $$\{p'(r),m'(r)\}$$ have some problems reagarding their numerical implementation as discussed in this question: The integration domain is not known at the beginning of numerical ODE integration and determining the stellar radius $$R$$ as $$p(R)=0$$ is for some EoS problematic. Additionally the system in $$\{p'(r),m'(r)\}$$ is very stiff and not very well conditioned. This has been known in literature for a long while and there is an -- in my opinion perfect solution -- a reformulation of the equations put forward in

L. Lindblom, Phase transitions and the mass-radius curves of relativistic stars, Phys. Rev. D 58.2 (1998) 024008, arXiv: gr-qc/9802072 [gr-qc],

and discussed in for example

S. Postnikov, Topics in the Physics and Astrophysics of Neutron Stars, Electronic Dissertation, Ohio University, 2010,

M. J. Steil, Structure of slowly rotating magnetized neutron stars in a perturbative approach, Masters Thesis, TU Darmstadt, 2017.

The last of the aforementioned references is my Master's thesis and this answer is basically based on Sec. 3.2.3 and App. B.3.

The idea is to reformulate the TOV equations in a thermodynamic variable to fix the integration range. For practical reasons the chemical potential or even better (due to Lindblom) the logarithm of the enthalpy per baryon $$h=\ln\left(\frac{\rho+p}{\bar{m}\, n}\right), \tag{1}$$ are good choices. In eq. (1) $$\rho$$ is the energy density, $$p$$ the pressure, $$\bar{m}$$ the mean particle mass, $$n$$ the particle number density and thus $$\bar{m}\,n$$ is the rest mass density. One can reformulate the TOV equation in terms of $$h$$ as \begin{align} \frac{\mathrm{d}r(h)^2}{\mathrm{d}h}&=- \frac{2 r(h)^2\left(1-2z(h)\right)}{4\pi r(h)^2 P(h)+z(h)}, \tag{2b}\\ \frac{\mathrm{d}z(h)}{\mathrm{d}h}&=\left(2\pi \rho(h)- \frac{z(h)}{2r(h)^2}\right)\frac{\mathrm{d}r(h)^2}{\mathrm{d}h}, \tag{2a} \end{align} with the radius squared $$r(h)^2$$ and the enclosed mass over radius $$z(h)\equiv m(r)/r$$. At the core one can expand this system in $$h_c-h$$: \begin{align} r(h)^2&= \frac{3}{2\pi(3p_c+\rho_c)}(h_c-h)+O((h_c-h)^2), \tag{3a}\\ z(h)&= \frac{2\rho_c}{3p_c+\rho_c}(h_c-h)+O((h_c-h)^2), \tag{3b} \end{align} which we can use to make a small step first step from $$h_c$$ to $$h_c-\Delta h$$ to start our numerical ODE integration from $$h_c-\Delta h$$ to avoid apparent divergencies in the ODE system (2) at $$h=h_c$$, where $$r(h_c)^2=0$$ and $$z(h_c)=0$$.

Lets see this in action for some polytropes \begin{align} k&=\kappa \bar{m} n_0^{1-\gamma},\tag{4a}\\ p(n)&=k n^\gamma, \tag{4b}\\ \rho(n)&= \bar{m}n+\frac{k}{\gamma-1}n^\gamma,\tag{4c}\\ n(h)&=\left(\frac{\bar{m}}{k}\frac{\gamma-1}{\gamma}(\mathrm{e}^h-1)\right)^{\frac{1}{\gamma-1}},\tag{4a} \end{align} with $$\bar{m}=1.66\times 10^{-27}\,\mathrm{kg}$$ and $$n_0=0.1\,\mathrm{fm}^{-3}$$. In the ODE system (2) we use $$p(n(h))$$ and $$\rho(n(h))$$ to close the system. Probing in $$h$$ between $$10^{-5}$$ and $$10$$ we obtain the following Mass-Radius curves

where dots mark the configurations with maximum mass per EoS, unstable configurations lie on the dotted branch, stable ones on the solid branch. Note that for $$\gamma=4/3$$ we find only unstable configurations -- i.e. the star with the larges mass on the M-R-curve is always the one with the biggest radius/lowest central pressure. This, the shape of the curves as well as Newtonian limits for those polytropes are discussed formally in

J. M. Heinzle, N. Rohr, and C. Uggla, Spherically symmetric relativistic stellar structures, Class. Quant. Grav. 20 (2003) 4567–4586, arXiv: gr-qc/0304012 [gr-qc].

The Mathematica code I used to compute this plot can be found here.

EDIT: To explain the mass-radius plot a bit better: here the corresponding plot for Mass over central density

For stable configuration $$M$$ rises with $$\rho_c$$. As per Eqs. (4b) and (4c): $$p_c<\rho_c$$ for all EoS here. The inspiraling of the M-R curves occurs at extreme densities and pressures -- $$10^4 \mathrm{MeV}\,\mathrm{fm}^{-3}$$ is around 100 times larger than nuclear density and exceed by far what we would expect even in the cores of the heaviest neutron stars -- nevertheless mathematically those are valid results, which are discussed at length by Heinzle et al. in Sec. 6. Mass-Radius Theorems and shown in Fig. 5 of the aforementioned reference.

• I updated the answer to make the MR plot and what is shown a bit clearer.
– N0va
Commented Aug 12, 2023 at 15:03
• Understandable: the implementation in $P(r)$ is just not very well suited for the more extreme cases/central pressures shown here. I can really recommend the one put forward here in $h$ and $\{r^2,z\}$. Its much easier to implement and one has to only adapt ones EoS to work in $h$.
– N0va
Commented Aug 12, 2023 at 15:32
• Is there other good choices of parametrization, instead of (1)-(3)?
– Cham
Commented Aug 12, 2023 at 18:22
• I am not aware of other popular options but one could use the EoS to get $\mathrm{d}h/\mathrm{d}q$, where $q$ is some thermodynamic variable and use it to rewrite Eqs. (2) as ODEs in $q$ and integrate in $q$. What makes $h$ so interesting especially in theoretical work ("inversion of TOV eq" -- the mapping between EoS and MR curve) is the numerical stability of the system and that Eqs. (3a) and (3b) are proper power series in $h_c-h$.
– N0va
Commented Aug 12, 2023 at 19:01
• Please, can you show a few plots for $1 < \gamma < \frac{4}{3}$? I'm getting some weird Mass-Radius relations, I don't know why yet.
– Cham
Commented Aug 24, 2023 at 0:17