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:
\begin{equation}\tag{1}
\frac{dp}{dr} = – G \frac{(M_r + 4 \pi p r^3)(\rho + p)}{r (r – 2 G M_r)},
\end{equation}
where $M_r(r)$ is the mass enclosed inside a sphere of radius $r$:
\begin{equation}\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,
\end{equation}
$\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:
\begin{equation}\tag{3}
p(\rho_{\text{mass}}) = \kappa \rho_{\text{mass}}^{\gamma},
\end{equation}
where $\gamma > 1$. Since $\rho = \rho_{\text{mass}} + p/(\gamma – 1)$, we may write the third equation as this:
\begin{equation}\tag{4}
\rho = \Bigl( \frac{p}{\kappa} \Bigr)^{1/\gamma} + \frac{p}{\gamma – 1}.
\end{equation}
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:
EDIT 2: Here's an example of the Mass-Radius relation curve:
Best Answer
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
and discussed in for example
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
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.