2
$\begingroup$

I'm trying to numerically integrate the TOV (Tolman-Openheimer-Volkof) equations, using Mathematica. The code works, but I'm having issues with the polytropic equation of state (EoS). The equations to be solved are these (I'm using units such that $c \equiv 1$): \begin{align} \frac{dp}{dr} &= - G \frac{(M_r + 4 \pi p r^3)(\rho + p)}{r (r - 2 G M_r)}, \tag{1} \\ \frac{dM_r}{dr} &= \rho 4 \pi r^2, \tag{2} \\ \rho &= f(p), \tag{3} \end{align} where $p$ is the local pressure, $\rho$ is the local total energy density (not the mass density!), and $M_r$ is the total mass-energy enclosed inside a sphere of internal radius $r$ ($M_r = \int_0^r \rho 4 \pi r^2 dr$). In principle, the pressure $p$ can be as high as we wish, but some energy conditions asks that $p \le \rho$. The non-general relativistic state is defined by the condition $p \ll \rho$ and equations (1) reduces to the Newtonian version (with $\rho \approx \rho_m$, the mass density): $$\tag{4} \frac{dp}{dr} = - \frac{G M_r \rho_m}{r^2}. $$

The usual -- non-relativistic -- EoS of a polytrope is \begin{equation}\tag{5} p(\rho_m) = \kappa \rho_m^{\gamma} \qquad \Rightarrow \qquad \rho_m = \Bigl( \frac{p}{\kappa} \Bigr)^{\frac{1}{\gamma}}, \end{equation} where $\rho_m$ is the mass density (not the total energy density). Since $\rho = \rho_m + \rho_{\text{int}}$, and $p = (\gamma - 1)\rho_{\text{int}}$, we could write the non-relativistic EoS as this: \begin{equation}\tag{6} \rho(p) = \Bigl( \frac{p}{\kappa} \Bigr)^{\frac{1}{\gamma}} + \frac{p}{\gamma - 1}. \end{equation} This is consistent with equation (8) of this answer: Is total energy density in energy-momentum tensor equal to total energy density with relativistic mass density?. So I now could solve the TOV equations. Then, I have an issue for the ultra-relativistic fluid of adiabatic index $\gamma = \frac{4}{3}$. Equation (6) doesn't reduce to $\rho = 3 p$ in this case, and I cannot trust (6) for high pressure states. My numerical resolution agrees with the newtonian case, when pressure is low: $p \ll \rho$, but I don't trust it for high pressure.

So what should be the proper EoS $\rho = f(p)$ for a general relativistic polytrope fluid at any pressure $p \le \rho$, that reduces to $p = \frac{1}{3} \rho$ for $\gamma = \frac{4}{3}$? I guess there is none, and that we need to try various extrapolations. I feel that the adiabatic index $\gamma$ is well defined only for low pressure fluids: $p \ll \rho$, am I right?

So, what is a "polytrope" in a general relativistic regime? A fluid of EoS $p = \kappa \rho^{\Gamma}$ ($\rho_m$ replaced by $\rho$ and $\gamma$ replaced by another number)?

$\endgroup$
2
  • $\begingroup$ I do not fully understand the question. A relativistic adiabatic index is defined as $\gamma\equiv \frac{\rho+p}{p}\frac{d P}{d\rho}$ which is consistent with equation (6). If one considers the limit of large pressures in equation (6) for $\gamma=4/3$ we recover $\rho\rightarrow 3p$ since $(p/\kappa)^{3/4}$ is subleading in the limit $p\rightarrow\infty$ against the term $3p$. What am I missing? $\endgroup$
    – N0va
    Aug 10, 2023 at 14:42
  • $\begingroup$ @N0va, from where does that definition of 𝛾 comes out? We could assume the EoS $P(\rho) = \kappa \rho^{\Gamma}$ with \begin{equation}\Gamma \approx \frac{\gamma}{1 + P_c / \rho_c}\end{equation} from a "relativistic definition" \begin{equation}\Gamma = \frac{\rho}{P} \frac{dP}{d\rho} = \mathrm{cste}.\end{equation}(from this EoS and using your numbers from your answer to the other question, I get the radius $R = 19.3 \, \mathrm{km}$). $\endgroup$
    – Cham
    Aug 10, 2023 at 15:08

1 Answer 1

2
$\begingroup$

Given the brief exchange in the comments bellow the question I think I see the issue/confusion. Lets get some definitions out of the way first. We follow J. Heinzle et. al., 2003, Spherically symmetric relativistic stellar structures and define the Newtonian adiabatic index as $$\gamma_\mathrm{N}\equiv\frac{\rho}{p}\frac{\mathrm{d}p}{\mathrm{d}\rho} \tag{1}$$ and the relativistic adiabatic index $$\gamma\equiv\frac{1}{c^2}\frac{\varepsilon+c^2p}{p}\frac{\mathrm{d}p}{\mathrm{d}\varepsilon}, \tag{2}$$ following Eqs. (4) and (15) of the aforementioned reference. Note however that I introduced the rest mass density $\rho$ and the energy density $\varepsilon$. Note that in the relativistic case

$$\varepsilon=\rho c^2+\epsilon,\tag{3}$$ with the internal energy density $\epsilon$ and thus $\varepsilon\neq\rho$. Further more we note that $\varepsilon$ is of order $O(c^2)$. The TOV equation reads

$$ \frac{\mathrm{d}p(r)}{\mathrm{d}r} =-\frac{G M(r)}{r^2}\frac{\varepsilon(r)}{c^2}\left(1+\frac{p(r)}{\varepsilon(r)}\right)\left(1+\frac{4 \pi r^3 p(r)}{c^2M(r)}\right) \left(1-\frac{2\,G M(r)}{c^2 r}\right)^{-1}, \tag{4}$$ where we separated GR corrections and explicitly maintained the constants $G$ and $c$, see, e.g., R. Silbar and S. Reddy, 2005, Neutron Stars for Undergraduates.

The Newtonian limit ($c\rightarrow\infty$) of the TOV equation follows as $$ \frac{\mathrm{d}p}{\mathrm{d}r} =-\frac{G M(r)}{r^2}\rho(r),\tag{5} $$ note the change from energy to (rest) mass density. Similarly we can confirm that $\lim_{c\rightarrow\infty}\gamma=\gamma_\mathrm{N}$.

Now regarding polytropic EoS: a Newtonian polytrope may be defined as $$ P(\rho)=K\left(\frac{\rho}{\rho_0}\right)^{\gamma_\mathrm{N}}\tag{6} $$ and in fact a relativistic polytrope, see, e.g., S. Bonazzola et al., 1993, Axisymmetric rotating relativistic bodies: A new numerical approach for ’exact’ solutions, may be defined as $$ P(\rho)=K\left(\frac{\rho}{\rho_0}\right)^{\gamma},\tag{7} $$ note that this gives the pressure as a function of rest mass density, which assuming a mean particle rest mass of $\bar{m}$ is related to the particle number density $n$ by $\rho=\bar{m}n$ and thus $$ p(n)=K\left(\frac{n \bar{m}}{n_0 \bar{m}}\right)^{\gamma}=K\left(\frac{n}{n_0}\right)^{\gamma}.\tag{8} $$ For the TOV equation and EoS in the form of $p(\rho)$ or conversely $p(n)$ is not useful -- to close the system one requires $\varepsilon(p)$. The first law of thermodynamics (at zero temperature) i.e. the Gibbs–Duhem equation $$ \varepsilon(n) = -p(n)+\mu(n) n=-p(n) +\frac{\mathrm{d}\varepsilon(n)}{\mathrm{d}n}n \tag{9} $$ provides the necessary bridge, since Eq. (9) can be integrated using $p(n)$ from Eq. (8) fixing the integration constant on dimensional grounds and in contact to Eq. (3) as $\bar{m}c^2$ yielding $$ \varepsilon(n) = \frac{K}{\gamma-1}\left(\frac{n}{n_0}\right)^{\gamma} +\bar{m}c^2 n.\tag{10} $$ Inverting Eq. (8) and inserting $n(p)$ into Eq. (10) results in $$ \varepsilon(p) = \bar{m}c^2n_0 \left(\frac{p}{K}\right)^{1/\gamma}+ \frac{p}{\gamma-1} \tag{11}. $$

So far so good, no? As for 'trusting' a polytrope (however one may define it): polytropic EoS are thermodynamic consistent, usually causal EoS which are very convenient numerically. What they are not (in case of a single polytrope) are realistic EoS for physical neutron stars. The density range inside such stars ranges from the one of Iron (in the crust) to densities exceeding nuclear density (in their cores). An EoS for such a fast density range can not be modeled with 3 parameters as a smooth curve. See, e.g, my answer to the question Role of numeric values of polytropic index. Getting realistic radii (for typical masses) is especially difficult since the polytropes lack information about the highly complicated EoS in the neutron star crust.

$\endgroup$
2
  • $\begingroup$ This is all fine, but I still dislike (11) for a numerical integration, since it gives some issues (depending on the parameters). Mathematica give me some troubles for $\gamma \le \frac{4}{3}$ and high central pressure, and also for $\gamma \approx 1$. Maybe this is just some physical aspects that I don't properly understand. Some authors prefer to use the other relativistic generalization of the EoS (see in particular articles.adsabs.harvard.edu/pdf/1973ApJ...183..637B): $P(\varepsilon) = \kappa \varepsilon^{\Gamma}$ (with $\Gamma$ as in my comment above). $\endgroup$
    – Cham
    Aug 10, 2023 at 20:33
  • 1
    $\begingroup$ If you are after numerical stability of integration over vast pressure ranges I can not recommend enough the reformulation of the TOV equations of Lindbolm (arxiv.org/abs/gr-qc/9802072 in the appendix): the idea is to integrate in the logarithm of the chemical potential instead of radius: this solves a lot of issues like knowing the integration range beforehand, a better controlled initial condition and improved numerical stability. $\endgroup$
    – N0va
    Aug 10, 2023 at 21:26

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.