Constructing M-R diagram from TOV equations with Python

im trying to solve the modified TOV equations dP'(r)/dr = -1.474 ε'(r)M'(r)/r^2 (1+P'(r)/ε'(r))(1+11.210^-6 r^3 P'(r)/M'(r))[(1-2.948M'(r)/r)]^-1 dM'(r)/dr = bo * ε'(r) where M(r)=M'(r)Mo, ε(r)=ε'(r)εο , εο= 1MeVfm^-3, GMo/c^2=1.474 Km and 4pi/(Moc^2)=0.710^-40 s^2/(kgKm^2),bo=8.910^-7 km^-3 with the equation of state ε(P)=15.55P^0.666+76.71P^0.247. These modified equations are correct for sure cause they are taken from 2 different papers. The problem i have is that the radius doesnt change with different values of Pc whereas M does change but a little bit. The range of Pc i use is the same as in the paper: Pc=1-1200 Mev*fm^-3. Can someone help me? This is the code im using:

import numpy as np
from scipy.integrate import solve_ivp

# Constants for the modified TOV equations
epsilon_0 = 1  # MeV/fm^3
b0 = 8.9e-7  # km^-3

# Equation of state: ε(P) = 15.55 * P^0.666 + 76.71 * P^0.247
def energy_density(P):
    return 15.55 * P**0.666 + 76.71 * P**0.247

# Modified TOV equations: dP/dr and dM/dr
def tov_equations(r, y):
    P, M = y
    epsilon_r = energy_density(P) * epsilon_0  # Energy density from the equation of state
    if P <= 0:
        return [0, 0]  # Stop when pressure becomes zero

    # Modified TOV equations
    dP_dr = -1.474 * (epsilon_r * M * (1 + P / epsilon_r)) / r**2 * \
            (1 + 11.2e-6 * r**3 * P / M) * (1 - 2.948 * M / r)**-1
    dM_dr = b0 * epsilon_r

    return [dP_dr, dM_dr]

# Boundary conditions at r=0: central pressure and mass
def tov_initial_conditions(Pc):
    M_c = 1e-2  # Small initial mass at the center
    return [Pc, M_c]

# Solve TOV equations
def solve_tov(Pc):
    # Initial conditions: central pressure and mass
    y0 = tov_initial_conditions(Pc)

    # Define the radial range for integration (we'll stop when pressure drops to ~0)
    r_max = 50  # Maximum radius to integrate out to [km]
    r_eval = np.linspace(1e-4, r_max, 10000)  # Evaluation points in radius

    # Solve TOV equations using solve_ivp (Runge-Kutta method)
    solution = solve_ivp(tov_equations, [1e-4, r_max], y0, method='RK45', t_eval=r_eval, rtol=1e-6, atol=1e-6)

    # Extract radius, pressure, and mass from the solution
    radius = solution.t
    pressure = solution.y[0]
    mass = solution.y[1]

    # Find the radius where pressure becomes negligible (surface of the star)
    surface_indices = np.where(pressure <= 1e-10)[0]
    if len(surface_indices) == 0:
        # If no valid surface is found, use the last point
        surface_index = -1
    else:
        surface_index = surface_indices[0]

    R_star = radius[surface_index]  # Radius at the surface
    M_star = mass[surface_index]  # Mass at the surface

    # Handle cases where R_star or M_star may still be unrealistic (very small or zero)
    if R_star <= 0 or M_star <= 0:
        return None, None

    return R_star, M_star

# Generate central pressure values from 1 to 1200 MeV/fm^3
Pc_values = np.linspace(1, 1200, 400)

# Arrays to store the results
R_values = []
M_values = []

# Loop over each central pressure value, solve the TOV equations, and store the results
for Pc in Pc_values:
    R_star, M_star = solve_tov(Pc)
    if R_star is not None and M_star is not None:
        R_values.append(R_star)
        M_values.append(M_star)

# Print a sample of the results
print("Sample results (Radius in km and Mass in dimensionless units):")
for i in range(0, len(R_values), 50):  # Print every 50th result for sampling
    print(f"Central pressure: {Pc_values[i]:.2f} MeV/fm^3 -> Radius: {R_values[i]:.2f} km, Mass: {M_values[i]:.2f}")
Sample results (Radius in km and Mass in dimensionless units):
Central pressure: 1.00 MeV/fm^3 -> Radius: 0.08 km, Mass: 0.31
Central pressure: 151.25 MeV/fm^3 -> Radius: 0.08 km, Mass: 0.26
Central pressure: 301.50 MeV/fm^3 -> Radius: 0.08 km, Mass: 0.37
Central pressure: 451.75 MeV/fm^3 -> Radius: 0.08 km, Mass: 0.30
Central pressure: 602.00 MeV/fm^3 -> Radius: 0.08 km, Mass: 0.27
Central pressure: 752.25 MeV/fm^3 -> Radius: 0.08 km, Mass: 0.25
Central pressure: 902.50 MeV/fm^3 -> Radius: 0.08 km, Mass: 0.23
Central pressure: 1052.75 MeV/fm^3 -> Radius: 0.08 km, Mass: 0.22

You need to sign in to view this answers