Solving TOV equations that describes neutron stars in modified f(R, T) gravity

Question

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.

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