Sorry for the long post, tldr at bottom.

I'm trying to use standard RK4 code in C/C++ to solve a coupled system of 2 modified TOV equations in f(R,T) gravity and reproduce some of the results of this paper. The two equations are:

$\displaystyle\frac{dm}{dr}=4\pi r^2\rho-\frac{r^2}{4c^2}h(T)$

$\displaystyle \frac{dP}{dr}=-\frac{G}{r^2}\left(\rho+\frac{P}{c^2}\right)\left[m+\frac{4\pi r^3P}{c^2}+\frac{r^3}{2}\right(\frac{h(T)}{2c^2}+\left(\rho+\frac P{c^2}\right)h_T\Bigg)\Bigg]\left(1-\frac{2Gm}{c^2r}\right)^{-1}$

where $h(T) = 2\kappa\lambda T + \kappa^2\xi T^2$ and $h_T$ is the derivative of $h(T)$ with respect to $T$, $h_T = 2\kappa\lambda + 2\kappa^2\xi T$ for parameters $\xi, \lambda$ and $\kappa = \displaystyle\frac{8\pi G}{c^4}$. $T = -\rho c^2+3p$ is the trace of the energy-momentum tensor. For $h(T) = 0$, we get the original TOV equations in General Relativity.

The idea to integrate is: start integration from some radius $r$ close to zero and continue integrating till pressure $p$ becomes zero. The $r$ value at which pressure becomes zero gives the radius of the star.

Since there are two equations and three unknowns $(r, \rho, p)$, we require another equation to solve the system and we use the polytropic equation of state $p=K\rho^\gamma$. Boundary conditions: $r = r_0$ which is close to zero, $\rho = \rho_C$ for some central density at $r = r_0$ and the pressure can be found from equation of state.

The paper uses $\gamma=2$ and $K = 7.1\times 10^{-18} kg^{-1} km^{-5} s^{-2}$. But my code isn't working with this value and it stops at the first iteration itself.

I used chatgpt to find that a good $K$ value for central density $\rho_C = 1.5\times10^{18} kg/m^3$ is $K = 2E-2$. I modified it to $K = 0.711E-2$ to match the results of original TOV equations (i.e. $h(T)=0$) of the paper. On using this value, the code runs okay and gives some lowest pressure $P_0$ after which it becomes negative. In order to not make the post too long, I'm attaching my code as pastebin link https://pastebin.com/pKsBwQ71.

Compared to the paper, I've only change the value of the constant $K$. One of my doubts is that the pressure value is of the order of $10^{23}$ before it becomes negative of about the same order in the next iteration. Shouldn't it be smooth? And when I change the paramters $\lambda, \xi$, I'm getting no change unlike the results of the paper. I've tried to contact the authors but got no response.

TLDR: problem with code to solve TOV equations in f(R, T) gravity. What should the adiabatic constant value for polytropic equation of state $p = K\rho^2$ if I'm including all factors like $G, c$ in SI units in the equations?

Any help would be appreciated.

  • $\begingroup$ "my code isn't working with this value and it stops at the first iteration itself." You don't understand your code then. What specifically happens when it stops? What does "isn't working" mean concretely? What have you already done to understand why it stops? You can change parameters all day long but how much confidence do you have in the correctness of your results if for some parameters the code just stops? $\endgroup$ Commented Jul 10 at 19:32
  • $\begingroup$ @WolfgangBangerth By isn't working, I mean it's not giving the desired results. There is a square root of the pressure term involved in the computation and thus the code will give nan if the pressure becomes negative, which it does at some point. I don't fully understand the background of the code, like underflow or overflow errors which might be happening given the range of values. I found that by increasing K, the code runs for a while before giving the final positive Pressure value after which P becomes negative. So it is possibly related to stability or stiffness of the system. $\endgroup$
    – hidenori
    Commented Jul 11 at 7:50
  • $\begingroup$ I continue to think that you don't understand what your code is doing. If it is computing the square root of a negative number, you need to figure out why that is happening, not just increase a parameter in hopes that it doesn't happen. I will also continue to ask "how much confidence do you have in the correctness of your results if for some parameters the code just stops?" $\endgroup$ Commented 2 days ago
  • $\begingroup$ @WolfgangBangerth Yeah, I'm trying to figure out why this is happening. As for the confidence, I'm not so sure I understand what you mean. The way I check the correctness of the results is by comparing the output of my code with the standard results available. Since, there aren't many free parameters in the code, I feel changing the ones that are free, would change the numerical stability of the system and thus change the way the code runs and its results. I have checked that the root structure of my code is logically okay. I'll try to check it again. $\endgroup$
    – hidenori
    Commented yesterday

1 Answer 1


(Edit from my previous clearly wrong answer) Good afternoon, For me the units of $K$ (lets call them $[x]$ ) in the paper are not clear, since for polytropic gases

$$[Pa]=[x][kg\; m^{-3}]$$ $$[kg\;m^{-1}\;s^{-2}]=[x][kg\; m^{-3}]$$ therefore $$[x]=[m^2\;s^{-2}]$$

Another way around is to consider it with the (more familiar to me) ideal gases EOS

$$p=\rho R T$$

with $R$ being the gas constant and $T$ the temperature. By comparison of ideal gases EOS and polytropic gases EOS we can conclude that, since $[J]=[kg\;m^2\;s^{-2}]$

$$[x]=[R T]=[J\;kg^{-1}\;K^{-1}\;K]=[m^2\;s^{-2}]$$

On the other side, why aren't you advancing $T$ (trace of the energy-momentum tensor) inside the RK loop? For me, maybe you are computing the stages $s>1$ using values of pressure from different levels. Apparently, after computing each stage $k_s$, you could re-compute $T$ with the same pressure value you will use to compute $k_{s+1}$.

If this is true, I do not expect that has more effect than lowering the convergence of the method.

PD: My apologizes again for not being commenting instead of posting as answer, my rank in this forum is not enough

  • $\begingroup$ The exponent of $km$ is $-5$ so on converting it to SI, $K$ will decrease further. Might be a typo in the preprint. But what I don't understand is, why my code is not working. Assuming it was a typo, I'm doing everything according to the paper. $\endgroup$
    – hidenori
    Commented Jul 10 at 6:56
  • $\begingroup$ A correction: the equation of state is $p=K\rho^2$. So you'd get $[x] = kg^{-1} m^{5}s^{-2}$. In the paper, it is $km^{-5}$ in the unit which I think is a typo (negative exponent). I tried advancing T in the loop, but it has no effect. $\endgroup$
    – hidenori
    Commented Jul 10 at 20:44

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.