EOSGenerators module guide

In this notebook, we will showcase all the EOS we currently implemented in our package.

RMF EOS

First import all of the package that will be used here.

import EOSgenerators.crust_EOS as crust
import EOSgenerators.RMF_EOS as RMF
import EOSgenerators.Polytrope_EOS as Polytrope
import EOSgenerators.Strangeon_EOS as Strangeon
import TOVsolver.main as main
import matplotlib.pyplot as plt
import numpy as np
import math
from TOVsolver.unit import g_cm_3, dyn_cm_2, km, Msun

(a) Load crust EOS and prepare interface EOS

Define constants, and load the crust EOS into this notebook, since when we solve the equation of state to TOV equation, we need full scope of the equation of state. Here, crust eos file is ‘Tolos_crust_out.txt’, while we need to times g_cm_3 for energy density and dyn_cm_2 for pressure, that will be easier for treatment, no need worry about G and c.

Then just connect the crust part with the interface part(the part between core equation of state and crust) by a PolyInterpolate function

Tolos_crust_out = np.loadtxt("Test_Case/Tolos_crust_out.txt")
eps_crust_T_out = Tolos_crust_out[:, 3] * g_cm_3
pres_crust_T_out = Tolos_crust_out[:, 4] * dyn_cm_2
eps_com, pres_com = crust.PolyInterpolate(eps_crust_T_out, pres_crust_T_out)

(b) Defining the paramters

Defining the free prameters that we need to compute the equation of state, you could check our intro part of documentation or just check the original paper of us about details. Huang et al,2023

After defined the equation of state parameters, we could call the ‘compute_EOS’ function from RMF, to compute full core equation of state.

m_sig = 495.0 / 197.33
m_w = 3.96544
m_rho = 3.86662

g_sigma = math.sqrt(107.5751)
g_omega = math.sqrt(182.3949)
g_rho = math.sqrt(206.4260)

kappa = 3.09114168 / 197.33
lambda_0 = -0.00168015405
zeta = 0.024
Lambda_w = 0.045
theta = np.array(
    [m_sig, m_w, m_rho, g_sigma, g_omega, g_rho, kappa, lambda_0, zeta, Lambda_w]
)
eps, pres = RMF.compute_EOS(eps_com, pres_com, theta)

(c) Connecting Core and crust EOS

These two line can easily connect the Equation of state we computed for core part by setting the parameters and the crust equation of state, to a full equation of state that prepared for next step output

eps_total = np.array([*eps_com, *eps])
pres_total = np.array([*pres_com, *pres])

Here below is one of the possiblity, use our defined function, ‘OutputMR’ to compute out what the mass radius curve corresponding to this equation of state.

MR = main.OutputMR("", eps_total, pres_total).T

fig, ax = plt.subplots(1, 1, figsize=(9, 6))
ax.plot(MR[0] / km, MR[1] / Msun, lw=2)
ax.set_ylabel(r"M [$M_{\odot}$]", fontsize=16)
ax.set_xlabel("R [km]", fontsize=16)
ax.set_xlim(8.0, 20.0)
ax.set_ylim(0, 3)
ax.tick_params(top=1, right=1, which="both", direction="in", labelsize=14)
ax.tick_params(top=1, right=1, which="both", direction="in", labelsize=14)
fig.tight_layout()
plt.show()

Illustrating the power of Numba accelerated EoS Generation

We shall now run the same code but with the Numba accelerated version, comparing the times for EoS Generation.

Keep in mind all previous definitions are still valid, all we need to change is a single line of code.

import EOSgenerators.fastRMF_EoS as FastRMF
import time

Tolos_crust_out = np.loadtxt("Tolos_crust_out.txt")
eps_crust_T_out = Tolos_crust_out[:, 3] * g_cm_3
pres_crust_T_out = Tolos_crust_out[:, 4] * dyn_cm_2
eps_com, pres_com = crust.PolyInterpolate(eps_crust_T_out, pres_crust_T_out)

m_sig = 495.0 / 197.33
m_w = 3.96544
m_rho = 3.86662

g_sigma = math.sqrt(107.5751)
g_omega = math.sqrt(182.3949)
g_rho = math.sqrt(206.4260)

kappa = 3.09114168 / 197.33
lambda_0 = -0.00168015405
zeta = 0.024
Lambda_w = 0.045
theta = np.array(
    [m_sig, m_w, m_rho, g_sigma, g_omega, g_rho, kappa, lambda_0, zeta, Lambda_w]
)

(a) Comparing the time it takes to generate the EoS

The reason why we test the time it takes to generate the EoS here is due to the fact that Numba is compiling the code, and as such the first run will take considerably longer to complete, but all other runs will be greatly faster.

Warning: Do not rerun this cell, as the compilation time will be gone.

### Normal RMF ####

first_run_time_normal = time.time()
eps_fast, pres_fast = RMF.compute_EOS(eps_com, pres_com, theta)
first_run_time_normal = time.time() - first_run_time_normal

time_for_20_runs_normal = time.time()
for _ in range(20):
    eps, pres = RMF.compute_EOS(eps_com, pres_com, theta)
time_for_20_runs_normal = time.time() - time_for_20_runs_normal


##### Fast RMF #####

first_run_time = time.time()
eps_fast, pres_fast = FastRMF.compute_EOS(eps_com, pres_com, theta)
first_run_time = time.time() - first_run_time

time_for_20_runs = time.time()
for _ in range(20):
    eps_fast, pres_fast = FastRMF.compute_EOS(eps_com, pres_com, theta)
time_for_20_runs = time.time() - time_for_20_runs


print(f"--------------------------------------------------------------------------------------")
print(f"|" + " "*84 + "|")
print(f"|" + " "*32 + "Time Comparion Table" + " "*32 + "|")
print(f"|" + " "*84 + "|")
print(f"|" + " "*17 + f"Time for first run  with RMF-Normal  = {first_run_time_normal*1000:6.0f} ms" + " "*19 + "|")
print(f"|" + " "*17 + f"Time for first run  with RMF-Fast    = {first_run_time*1000       :6.0f} ms" + " "*19 + "|")
print(f"|" + " "*84 + "|")
print(f"|" + " "*17 + f"Time for   20  runs with RMF-Normal  = {time_for_20_runs_normal*1000:6.0f} ms" + " "*19 + "|")
print(f"|" + " "*17 + f"Time for   20  runs with RMF-Fast    = {time_for_20_runs*1000       :6.0f} ms" + " "*19 + "|")
print(f"|" + " "*84 + "|")
print(f"|" + " "*17 + f"Time for  all  runs with RMF-Normal  = {(first_run_time_normal + time_for_20_runs_normal)*1000:6.0f} ms" + " "*19 + "|")
print(f"|" + " "*17 + f"Time for  all  runs with RMF-Fast    = {(first_run_time        + time_for_20_runs       )*1000:6.0f} ms (Has Compilation) " + "|")
print(f"|" + " "*84 + "|")
print(f"|" + " "*17 + f"Time for   1   run  with RMF-Normal  = {time_for_20_runs_normal*1000/20:6.0f} ms  " + " "*17 + "|")
print(f"|" + " "*17 + f"Time for   1   run  with RMF-Fast    = {time_for_20_runs       *1000/20:6.0f} ms" + " "*19 + "|")
print(f"|" + " "*84 + "|")
print(f"|" + " "*84 + "|")
print(f"|" + " "*30 + f"Total Speed Up = {time_for_20_runs_normal/time_for_20_runs:.2f} x" + " "*30 + "|" )
print(f"--------------------------------------------------------------------------------------")

--------------------------------------------------------------------------------------
|                                                                                    |
|                                Time Comparion Table                                |
|                                                                                    |
|                 Time for first run  with RMF-Normal  =    700 ms                   |
|                 Time for first run  with RMF-Fast    =   2525 ms                   |
|                                                                                    |
|                 Time for   20  runs with RMF-Normal  =  14176 ms                   |
|                 Time for   20  runs with RMF-Fast    =    194 ms                   |
|                                                                                    |
|                 Time for  all  runs with RMF-Normal  =  14875 ms                   |
|                 Time for  all  runs with RMF-Fast    =   2719 ms (Has Compilation) |
|                                                                                    |
|                 Time for   1   run  with RMF-Normal  =    709 ms                   |
|                 Time for   1   run  with RMF-Fast    =     10 ms                   |
|                                                                                    |
|                                                                                    |
|                              Total Speed Up = 73.13 x                              |
--------------------------------------------------------------------------------------

eps_total_fast = np.array([*eps_com, *eps_fast])
pres_total_fast = np.array([*pres_com, *pres_fast])

compare_fig1, ax = plt.subplots(1, 1, figsize=(9,6))

ax.plot(eps_total_fast, pres_total_fast, lw=2, label="Fast Version")
ax.plot(eps_total, pres_total, lw=6, label="Normal Version", alpha=0.3)
ax.set_ylabel(r"P $[fm^{-4}]$", fontsize=16)
ax.set_xlabel(r"$\epsilon$ $[fm^{-4}]$", fontsize=16)
ax.set_xlim(0, 6)
ax.set_ylim(0, 2)
ax.tick_params(top=1, right=1, which="both", direction="in", labelsize=14)
ax.tick_params(top=1, right=1, which="both", direction="in", labelsize=14)
plt.legend()
compare_fig1.tight_layout()
plt.show()

MR_fast = main.OutputMR("", eps_total_fast, pres_total_fast).T

compare_fig2, ax = plt.subplots(1, 1, figsize=(9, 6))

ax.plot(MR_fast[0] / km, MR_fast[1] / Msun, lw=2, label="Fast Version")
ax.plot(MR[0] / km, MR[1] / Msun, lw=6, label="Normal Version", alpha=0.3)
ax.set_ylabel(r"M [$M_{\odot}$]", fontsize=16)
ax.set_xlabel("R [km]", fontsize=16)
ax.set_xlim(8.0, 20.0)
ax.set_ylim(0, 3)
ax.tick_params(top=1, right=1, which="both", direction="in", labelsize=14)
ax.tick_params(top=1, right=1, which="both", direction="in", labelsize=14)
plt.legend()
compare_fig2.tight_layout()
plt.show()

Polytrope EOS

First import all the package that will be used.

import numpy as np
import math
import TOVsolver.main as main
import matplotlib.pyplot as plt
from scipy.integrate import ode
import EOSgenerators.Polytrope_EOS as Polytrope
import EOSgenerators.Strangeon_EOS as Strangeon
import EOSgenerators.crust_EOS as crust
from TOVsolver.unit import  km, Msun, MeV,fm,g_cm_3,dyn_cm_2, G,c

The polytropic equation of state (Polytrope) is the most commonly used EOS in astrophysics, its advantage is assumed to be capable to simulate all EOS shape with given parameters. Here, we provide one of the computation for polytrope.

(a) Connecting crust with inner crust

For polytrope, we first connecting it with crust as we did in RMF

Tolos_crust_out = np.loadtxt("Test_Case/Tolos_crust_out.txt")
eps_crust_T_out = Tolos_crust_out[:, 3] * g_cm_3
pres_crust_T_out = Tolos_crust_out[:, 4] * dyn_cm_2
eps_com, pres_com = crust.PolyInterpolate(eps_crust_T_out, pres_crust_T_out)

Define a constant “Saturation density”, so that we know the EOS span from which lower limit to upper limit.

rho_ns = 267994004080000.03# saturation density in cgs

We define a set of parameters, first three in theta are the parameters control the slope of each polytrope, and the last two are transition point from one segment to another. Details see the documentation about polytrope

theta = [1.0, 2.5, 3.0, 1*rho_ns, 4.5*rho_ns]# define the parameters
density = np.logspace(14.3,15.6,1000) # the core density from 10^(14.3) g/cm^3 to 10^(15.6) g/cm^3
pressures = Polytrope.compute_EOS(density, theta)
eps_total = np.array([*eps_com, *density*g_cm_3])
pres_total = np.array([*pres_com, *pressures])

(b) Solve TOV with this EOS

Here below we use the strangeon matter EOS to compute the mass radius curve.

The following code calculates the mass and radius of the strange stars for the given EOS of the polytrope using our main.OutputMR

MR = main.OutputMR('',eps_total,pres_total)

fig, ax = plt.subplots(1, 1, figsize=(9, 6))
ax.plot(MR[:,0]/km , MR[:,1]/Msun, lw=2)

ax.set_ylabel(r"M [$M_{\odot}$]", fontsize=16)
ax.set_xlabel("R [km]", fontsize=16)
ax.set_xlim(8,16)
ax.tick_params(top=1, right=1, which="both", direction="in", labelsize=14)
ax.tick_params(top=1, right=1, which="both", direction="in", labelsize=14)
fig.tight_layout()
plt.show()

Strangeon matter EOS

First import all the package that will be used.

import numpy as np
import math
import TOVsolver.main as main
import matplotlib.pyplot as plt
from scipy.integrate import ode
from EOSgenerators import EOS

from TOVsolver.unit import  km, Msun, MeV,fm

The strangeon matter EOS describes the strongly interacting matter in the solid state, and the compact star can be a self-bound quark star composed of strangeon matter. Therefore, in the following we present the bare quark star EOS without a crust.

Note that the strangeon matter parameters we input are in the natural unit system. The quentity returned by the Strangeon_compute_EOS function is in the Geometric Unit System. Thus, here, the parameter epsilon and ns are in the units of \(MeV/fm^{-3}\). So follow the unite conversion rule, if we are using this unit, we just time this unit to do compute

(a) Defining the EOS parameters

Define the strangeon matter EOS that we will use to calculate the energy densities and pressures of the system. The following is an example of how to define the theta and Nq. For example, theta=[50, 0.24] means the example values for epsilon and ns: epsilon in units of MeV, ns in units of \(fm^-3\). Nq is an integer, e.g. Nq=9, 12, 15, 18, 21, 24, 27

For the definition of parameter n, the minimum value is 3ns/Nq, the maximum value is about 0.168*3/Nq. Then we can generate,for example, 1000 points for the input of the srangeon matter EOS. Special note: If you are using this EOS, need a very fine grid for the EOS. From our test, 1000 points for this EOS is the minimum requirement

Once we define the EOS parameters already, we can use the Strangeon_compute_EOS function to calcute the energy_density and pressure.

Nq=18
epsilon=50
ns=0.24
theta = np.array([Nq, epsilon, ns])
n_min = 3 * theta[2] / theta[0]
n_max = 0.16 * 8 * 3 / theta[0]
n_values = np.linspace(n_min, n_max, 2000)

energy_densities, pressures = Strangeon.compute_EOS(n_values, theta)
print("n values:", n_values)
print("Energy densities:", energy_densities)
print("Pressures:", pressures)

n values: [0.04       0.04008671 0.04017342 ... 0.21315991 0.21324662 0.21333333]
Energy densities: [  204.61935484   205.0631334    205.50734538 ... 46616.1606005
 46711.98281964 46807.96367981]
Pressures: [0.00000000e+00 1.98864796e-01 4.00749542e-01 ... 1.88749504e+05
 1.89139331e+05 1.89529799e+05]

(b) Solve TOV with this EOS

Here below we use the strangeon matter EOS to compute the mass radius curve.

The following code calculates the mass and radius of the strange stars for the given EOS of the strange matter EOS. Since our energy_density and pressure are all in same unit \(MeV/fm^{-3}\) so all we need is to time this unit

Input a given central pressure and central energy density, you will obtain the radius where the pressure is zero and the mass at that radius.

Solve the TOV equations using each central pressure and energy density

MR= main.OutputMR('',energy_densities* MeV/fm**3   , pressures* MeV/fm**3)

fig, ax = plt.subplots(1, 1, figsize=(9, 6))
ax.plot(MR[:,0]/km , MR[:,1]/Msun, lw=2)

ax.set_ylabel(r"M [$M_{\odot}$]", fontsize=16)
ax.set_xlabel("R [km]", fontsize=16)
# ax.set_xlim(8.0, 20.0)
# ax.set_ylim(0, 3)
ax.tick_params(top=1, right=1, which="both", direction="in", labelsize=14)
ax.tick_params(top=1, right=1, which="both", direction="in", labelsize=14)
fig.tight_layout()
plt.show()

Speed of sound EOS

First import all the package that will be used.

import numpy as np
import TOVsolver.main as main
import matplotlib.pyplot as plt
from EOSgenerators import SpeedofSound_EOS
import EOSgenerators.crust_EOS as crust
from TOVsolver.unit import g_cm_3, dyn_cm_2, km, Msun

The Speed of sound EOS describes the core EOS of a compact star, so it should be connected with the crust EOS to form a full EOS. See https://arxiv.org/abs/1812.08188 for details.

Tolos_crust_out = np.loadtxt("Test_Case/Tolos_crust_out.txt")
eps_crust_T_out = Tolos_crust_out[:, 3] * g_cm_3
pres_crust_T_out = Tolos_crust_out[:, 4] * dyn_cm_2
eps_com, pres_com = crust.PolyInterpolate(eps_crust_T_out, pres_crust_T_out)

(a) Defining the EOS parameters

To construct the Speed of sound EOS, we need to specify the outer crust EOS and the interface EOS, and then connect them with the core EOS.

x_last = eps_com[-1]
y_last = pres_com[-1]
dydx_last = (pres_com[-1] - pres_com[-2]) / (eps_com[-1] - eps_com[-2])
CS_EOS = SpeedofSound_EOS.compute_EOS(x_last, y_last, dydx_last)

In order to generate the parameters of the Speed of sound EOS, we need a list of 5 uniform random numbers between 0 and 1.

After generate the parameters, we use function check_a to check if the parameters are valid.

cs_a = CS_EOS.gen_a((0.2, 0.2, 0.3, 0.4, 0.5))
print(CS_EOS.check_a(cs_a))

True

(b) Solve TOV with this EOS

Here below we use the speed of sound EOS to compute the mass radius curve.

First we calculate the core EOS.

core_e_min = x_last
core_e_max = 2e16 * g_cm_3
n_core_e = 1000
core_e = np.geomspace(core_e_min, core_e_max, n_core_e)
core_p = CS_EOS.cal_core_p(core_e, cs_a)

Then we concat the core EOS with the crust EOS.

full_e = np.concatenate((eps_com, core_e[1:]))
full_p = np.concatenate((pres_com, core_p[1:]))
plt.figure(dpi=200)
plt.plot(full_e / g_cm_3, full_p / dyn_cm_2)
plt.xscale("log")
plt.yscale("log")
plt.xlabel(r"$\epsilon\quad [g/cm^3]$")
plt.ylabel(r"$P\quad [dyn/cm^2]$")
plt.show()

Finally we compute the mass radius curve.

MR = main.OutputMR("", full_e, full_p)
plt.figure(dpi=200)
plt.plot(MR[:, 0] / km, MR[:, 1] / Msun)
plt.xlim(6, 17)
plt.ylim(0, 2)
plt.xlabel("Radius [km]")
plt.ylabel("Mass [Msun]")
plt.show()