Constructing EOS and solving TOV equations in scalar tensor theory#
⚠️ This solver is still in development. Therefore, expect breaking changes in this notebook.
This example notebook shows how to construct the equation of state with the metamodel and speed-of-sound extension scheme parametrization used in the paper, as well as solve the TOV equations.
[1]:
import matplotlib.pyplot as plt
params = {
"text.usetex": True,
"font.family": "serif",
"font.serif": ["Computer Modern Serif"],
"xtick.labelsize": 16,
"ytick.labelsize": 16,
"axes.labelsize": 16,
"legend.fontsize": 16,
"legend.title_fontsize": 16,
}
plt.rcParams.update(params)
import jax.numpy as jnp
from jesterTOV.eos.metamodel.metamodel_CSE import MetaModel_with_CSE_EOS_model
from jesterTOV.eos.families import construct_family_ST, construct_family_ST_sol
import jesterTOV.utils as utils
---------------------------------------------------------------------------
ModuleNotFoundError Traceback (most recent call last)
Cell In[1], line 17
15 import jax.numpy as jnp
16 from jesterTOV.eos.metamodel.metamodel_CSE import MetaModel_with_CSE_EOS_model
---> 17 from jesterTOV.eos.families import construct_family_ST, construct_family_ST_sol
18 import jesterTOV.utils as utils
ModuleNotFoundError: No module named 'jesterTOV.eos.families'
Equation of state#
[ ]:
nsat = 0.16 # nuclear saturation density in fm^-3
# Define the EOS object, here we focus on Metamodel with CSE
eos = MetaModel_with_CSE_EOS_model(nmax_nsat=6.0)
# Define the nuclear empirical parameters (NEPs) -- all in MeV
NEP_dict = {
"E_sat": -16.0, # saturation parameters
"K_sat": 200.0,
"Q_sat": 0.0,
"Z_sat": 0.0,
"E_sym": 32.0, # symmetry parameters
"L_sym": 40.0,
"K_sym": -100.0,
"Q_sym": 0.0,
"Z_sym": 0.0,
}
# Define the breakdown density -- this is usually between 1-2 nsat
nbreak = 1.5 * nsat
NEP_dict["nbreak"] = nbreak
# Then we extend with some CSE grid points
ngrids = jnp.array([2.0, 3.0, 4.0, 5.0]) * nsat
cs2grids = jnp.array([0.5, 0.4, 0.3, 0.2]) # speed of sound squared at the grid points
# Now create the EOS -- returns a tuple with most useful EOS quantities
ns, ps, hs, es, dloge_dlogps, mu, cs2 = eos.construct_eos(NEP_dict, ngrids, cs2grids)
# Make a plot
plt.subplots(nrows=2, ncols=2, figsize=(12, 10))
# For the plot, let's make some conversions to more common units
ns_plots = ns / utils.fm_inv3_to_geometric / 0.16
es_plots = es / utils.MeV_fm_inv3_to_geometric
ps_plots = ps / utils.MeV_fm_inv3_to_geometric
# p(n)
plt.subplot(221)
plt.plot(ns_plots, ps_plots)
plt.xlabel(r"$n$ [$n_{\rm{sat}}$]")
plt.ylabel(r"$p$ [MeV/fm$^3$]")
# e(n)
plt.subplot(222)
plt.plot(ns_plots, es_plots)
plt.xlabel(r"$n$ [$n_{\rm{sat}}$]")
plt.ylabel(r"$e$ [MeV/fm$^3$]")
# cs2(n)
plt.subplot(223)
plt.plot(ns_plots, cs2)
plt.xlabel(r"$n$ [$n_{\rm{sat}}$]")
plt.ylabel(r"$c_s^2$")
plt.axvline(0.5, color="red", label="Crust-core transition")
plt.axvline(nbreak / nsat, color="black", label=r"$n_{\rm{break}}$")
plt.legend()
# p(e)
plt.subplot(224)
plt.plot(es_plots, ps_plots)
plt.xlabel(r"$e$ [MeV/fm$^3$]")
plt.ylabel(r"$p$ [MeV/fm$^3$]")
plt.tight_layout()
# plt.show() # uncomment to see the EOS
plt.close()
Neutron stars#
Solving the TOV equations to construct an M(R) curve
NOTE: Tidal deformabilities work in progress
[ ]:
# Solve TOV equations:
beta_ST = -4.5
# initial guesses
phi_c = -0.1
nu_c = -1.2
# Enter through EoS
eos_tuple = (ns, ps, hs, es, dloge_dlogps, beta_ST, phi_c, nu_c)
logpc, masses, radii, Lambdas = construct_family_ST(eos_tuple, ndat=50, min_nsat=1.0)
# Make a plot
plt.subplots(nrows=1, ncols=2, figsize=(12, 6))
# Limit masses to be above certain mass to make plot prettier
m_min = 0.5
mask = masses > m_min
masses = masses[mask]
radii = radii[mask]
Lambdas = Lambdas[mask]
# M(R) plot
plt.subplot(121)
plt.plot(radii, masses)
plt.xlabel(r"$R$ [km]")
plt.ylabel(r"$M$ [$M_\odot$]")
# Lambda(R) plot
plt.subplot(122)
plt.plot(masses, Lambdas)
plt.xlabel(r"$M$ [$M_\odot$]")
plt.ylabel(r"$\Lambda$")
plt.yscale("log")
plt.tight_layout()
plt.show()
plt.close()
/var/folders/xj/6sk96vdn15385fpjn9nd1rcw0000gp/T/ipykernel_57185/3867325403.py:30: UserWarning: Data has no positive values, and therefore cannot be log-scaled.
plt.yscale("log")
Solve the TOV equations for a single point in debug mode, to recover full details on the constructed solution
[ ]:
# Solve TOV equations:
beta_ST = -4.5
phi_c = -0.1
nu_c = -1.2
eos_tuple = (ns, ps, hs, es, dloge_dlogps, beta_ST, phi_c, nu_c)
logpc, masses, radii, Lambdas, sol_iter, sol_ext = construct_family_ST_sol(
eos_tuple, ndat=1, min_nsat=3.0
)
# Limit masses to be above certain mass to make plot prettier
m_min = 0.5
mask = masses > m_min
masses = masses[mask]
radii = radii[mask]
Lambdas = Lambdas[mask]
# Check out first batch
r_vals_float = sol_ext.ts[0]
M_vals_float = sol_ext.ys[0][0]
nu_vals_float = sol_ext.ys[1][0]
phi_vals_float = sol_ext.ys[2][0]
psi_vals_float = sol_ext.ys[3][0]
# r_vals_float = jnp.array(r_vals/1e3).tolist()
# M_vals_float = jnp.array(M_vals/utils.solar_mass_in_meter).tolist()
# nu_vals_float = jnp.array(nu_vals).tolist()
# phi_vals_float = jnp.array(phi_vals).tolist()
# psi_vals_float = jnp.array(psi_vals).tolist()
Iteration 0: ν∞=0.03536283122108169, φ∞=-0.00010855321412178872,νc=-1.2282902649768652, φc=-0.09991315742870258, M=2.0413830711185046
Iteration 1: ν∞=0.007147314621028367, φ∞=-0.00010416438302998107,νc=-1.2340081166736878, φc=-0.0998298259222786, M=2.041380728636001
Iteration 2: ν∞=0.0015011395234137959, φ∞=-9.996345995740054e-05,νc=-1.2352090282924189, φc=-0.09974985515431267, M=2.0413784871415914
Iteration 3: ν∞=0.0003689667411880949, φ∞=-9.594145331237605e-05,νc=-1.2355042016853692, φc=-0.09967310199166277, M=2.0413763461960315
Iteration 4: ν∞=0.00013988832495967648, φ∞=-9.2076033475678e-05,νc=-1.235616112345337, φc=-0.09959944116488223, M=2.041374653124587
Iteration 5: ν∞=9.157991077841086e-05, φ∞=-8.835654908023917e-05,νc=-1.2356893762739598, φc=-0.09952875592561804, M=2.0413734832384187
Iteration 6: ν∞=7.932350852884584e-05, φ∞=-8.479352629940824e-05,νc=-1.235752835080783, φc=-0.0994609211045785, M=2.0413723922143316
Iteration 7: ν∞=7.438966251948026e-05, φ∞=-8.13798634370388e-05,νc=-1.2358123468107984, φc=-0.09939581721382887, M=2.0413713770367776
Iteration 8: ν∞=7.102415319759875e-05, φ∞=-7.810898810100422e-05,νc=-1.2358691661333565, φc=-0.09933333002334807, M=2.041370427434754
Iteration 9: ν∞=6.807289457133039e-05, φ∞=-7.49745676749436e-05,νc=-1.2359236244490135, φc=-0.09927335036920812, M=2.0413695373257332
Iteration 10: ν∞=6.523523883028332e-05, φ∞=-7.197612011862787e-05,νc=-1.2359758126400777, φc=-0.09921576947311321, M=2.0413685638616945
Iteration 11: ν∞=6.240817401103205e-05, φ∞=-6.91168442296457e-05,νc=-1.2360257391792866, φc=-0.0991604759977295, M=2.041367275988285
Iteration 12: ν∞=5.985951488451935e-05, φ∞=-6.637556502278767e-05,νc=-1.2360736267911943, φc=-0.09910737554571127, M=2.041366045236096
Iteration 13: ν∞=5.745010180609623e-05, φ∞=-6.374707443054618e-05,νc=-1.2361195868726391, φc=-0.09905637788616684, M=2.041364868859203
Iteration 14: ν∞=5.514821233632721e-05, φ∞=-6.122642424969722e-05,νc=-1.236163705442508, φc=-0.09900739674676708, M=2.041363744087284
Iteration 15: ν∞=5.294390619990609e-05, φ∞=-5.8808906596874785e-05,νc=-1.236206060567468, φc=-0.09896034962148957, M=2.0413626683381287
Iteration 16: ν∞=5.083177294205616e-05, φ∞=-5.649004024195568e-05,νc=-1.2362467259858216, φc=-0.098915157589296, M=2.0413616391878286
Iteration 17: ν∞=4.8807479673988134e-05, φ∞=-5.426555816861293e-05,νc=-1.2362857719695608, φc=-0.09887174514276112, M=2.041360654358447
Iteration 18: ν∞=4.6867070835718934e-05, φ∞=-5.213139592007471e-05,νc=-1.2363232656262293, φc=-0.09883004002602506, M=2.0413597117069884
Iteration 19: ν∞=4.5006818787294255e-05, φ∞=-5.0083680691150934e-05,νc=-1.2363592710812592, φc=-0.09878997308147214, M=2.0413588092150166
Iteration 20: ν∞=4.32216828594558e-05, φ∞=-4.811882988960944e-05,νc=-1.2363938484275467, φc=-0.09875147801756046, M=2.0413579421489185
Iteration 21: ν∞=4.150699390469701e-05, φ∞=-4.6233506444196356e-05,νc=-1.2364270540226705, φc=-0.0987144912124051, M=2.0413571038474805
Iteration 22: ν∞=3.986630242508015e-05, φ∞=-4.442403024035638e-05,νc=-1.2364589470646106, φc=-0.0986789519882128, M=2.0413563007218856
Iteration 23: ν∞=3.829337552601112e-05, φ∞=-4.2687200695885874e-05,νc=-1.2364895817650314, φc=-0.09864480222765609, M=2.041355531152628
Iteration 24: ν∞=3.678467653474792e-05, φ∞=-4.101996342541244e-05,νc=-1.2365190095062593, φc=-0.09861198625691577, M=2.0413547936074417
Iteration 25: ν∞=3.533731640917672e-05, φ∞=-3.9419402757695775e-05,νc=-1.2365472793593866, φc=-0.09858045073470961, M=2.041354086635297
Iteration 26: ν∞=3.394863831510998e-05, φ∞=-3.7882734682726806e-05,νc=-1.2365744382700388, φc=-0.09855014454696343, M=2.0413534088611365
Iteration 27: ν∞=3.260925863869755e-05, φ∞=-3.6407769825948856e-05,νc=-1.2366005256769497, φc=-0.09852101833110266, M=2.0413527468107167
Iteration 28: ν∞=3.1326550643713926e-05, φ∞=-3.499167944350022e-05,νc=-1.2366255869174647, φc=-0.09849302498754786, M=2.041352106234289
Iteration 29: ν∞=3.009934846283759e-05, φ∞=-3.3631803048930797e-05,νc=-1.236649666396235, φc=-0.09846611954510871, M=2.0413514918137325
Iteration 30: ν∞=2.8922216687228075e-05, φ∞=-3.2325823774527636e-05,νc=-1.2366728041695847, φc=-0.0984402588860891, M=2.041350902406037
Iteration 31: ν∞=2.77924300155258e-05, φ∞=-3.107152661559323e-05,νc=-1.236695038113597, φc=-0.09841540166479662, M=2.0413503369245674
Iteration 32: ν∞=2.6707879409803605e-05, φ∞=-2.9866793492550747e-05,νc=-1.2367164044171248, φc=-0.09839150823000258, M=2.0413497943355643
Iteration 33: ν∞=2.566664971879627e-05, φ∞=-2.870959856970354e-05,νc=-1.2367369377369, φc=-0.09836854055114681, M=2.041349273655239
Iteration 34: ν∞=2.4666931560984396e-05, φ∞=-2.7598003840984757e-05,νc=-1.2367566712821487, φc=-0.09834646214807402, M=2.041348773946855
Iteration 35: ν∞=2.370700031442208e-05, φ∞=-2.6530154957077313e-05,νc=-1.2367756368824003, φc=-0.09832523802410836, M=2.041348294317955
Iteration 36: ν∞=2.2786450642334763e-05, φ∞=-2.5504189935429992e-05,νc=-1.2367938660429143, φc=-0.09830483467216002, M=2.0413478361666852
Iteration 37: ν∞=2.1904012478816048e-05, φ∞=-2.451831810714738e-05,νc=-1.2368113892528974, φc=-0.0982852200176743, M=2.0413474011405635
Iteration 38: ν∞=2.105407380985141e-05, φ∞=-2.3571113464122113e-05,νc=-1.2368282325119453, φc=-0.098266363126903, M=2.041346983491466
Iteration 39: ν∞=2.0237265567360193e-05, φ∞=-2.2661017521494637e-05,νc=-1.2368444223243993, φc=-0.0982482343128858, M=2.0413465824921033
Iteration 40: ν∞=1.9452624569598178e-05, φ∞=-2.1786537740317292e-05,νc=-1.2368599844240549, φc=-0.09823080508269355, M=2.041346197448373
Iteration 41: ν∞=1.8698916003766046e-05, φ∞=-2.094624451166058e-05,νc=-1.2368749435568578, φc=-0.09821404808708423, M=2.041345827697484
Iteration 42: ν∞=1.7974895879416903e-05, φ∞=-2.0138768277966584e-05,νc=-1.2368893234735614, φc=-0.09819793707246185, M=2.041345472606595
Iteration 43: ν∞=1.72793610753771e-05, φ∞=-1.9362796829271845e-05,νc=-1.2369031469624217, φc=-0.09818244683499844, M=2.0413451315709232
Iteration 44: ν∞=1.661115774541865e-05, φ∞=-1.8617072699579878e-05,νc=-1.236916435888618, φc=-0.09816755317683877, M=2.041344804012772
Iteration 45: ν∞=1.5969180641695188e-05, φ∞=-1.790039073275938e-05,νc=-1.2369292112331314, φc=-0.09815323286425257, M=2.041344489379729
Iteration 46: ν∞=1.5352371531555465e-05, φ∞=-1.7211595739755346e-05,νc=-1.2369414931303566, φc=-0.09813946358766076, M=2.041344187143673
Iteration 47: ν∞=1.4759717018651523e-05, φ∞=-1.6549580283407487e-05,νc=-1.2369533009039715, φc=-0.09812622392343404, M=2.0413438967995465
Iteration 48: ν∞=1.4190246453212301e-05, φ∞=-1.5913282575292695e-05,νc=-1.236964653101134, φc=-0.0981134932973738, M=2.04134361786415
Iteration 49: ν∞=1.3643030056482292e-05, φ∞=-1.530168446650071e-05,νc=-1.2369755675251792, φc=-0.0981012519498006, M=2.041343349875154
Iteration 50: ν∞=1.3117177053293457e-05, φ∞=-1.4713809540230986e-05,νc=-1.2369860612668218, φc=-0.09808948090216842, M=2.0413430923901266
Iteration 51: ν∞=1.2611833833742767e-05, φ∞=-1.4148721300413726e-05,νc=-1.2369961507338887, φc=-0.09807816192512808, M=2.04134284498545
Iteration 52: ν∞=1.212618243430475e-05, φ∞=-1.3605521434173566e-05,νc=-1.2370058516798361, φc=-0.09806727750798075, M=2.0413426072556295
Iteration 53: ν∞=1.1659438766967005e-05, φ∞=-1.3083348173654184e-05,νc=-1.2370151792308497, φc=-0.09805681082944183, M=2.041342378812239
Iteration 54: ν∞=1.121085129686416e-05, φ∞=-1.258137471648376e-05,νc=-1.237024147911887, φc=-0.09804674572966864, M=2.041342159283341
Iteration 55: ν∞=1.0779699453824761e-05, φ∞=-1.2098807736902692e-05,νc=-1.2370327716714502, φc=-0.09803706668347911, M=2.04134194831257
Iteration 56: ν∞=1.0365292314282615e-05, φ∞=-1.1634885958566849e-05,νc=-1.2370410639053016, φc=-0.09802775877471226, M=2.0413417455584217
Iteration 57: ν∞=9.966967403246844e-06, φ∞=-1.1188878786306096e-05,νc=-1.2370490374792242, φc=-0.09801880767168321, M=2.0413415506937564
Iteration 58: ν∞=9.584089293514187e-06, φ∞=-1.076008501625124e-05,νc=-1.237056704750659, φc=-0.0980101996036702, M=2.0413413634049724
Iteration 59: ν∞=9.216048490528475e-06, φ∞=-1.0347831596891759e-05,νc=-1.2370640775894515, φc=-0.0980019213383927, M=2.041341183391394
Iteration 60: ν∞=8.862260362846623e-06, φ∞=-9.95147244329606e-06,νc=-1.2370711673977417, φc=-0.09799396016043806, M=2.0413410103647744
Iteration 61: ν∞=8.522164083348861e-06, φ∞=-9.570387304690548e-06,νc=-1.2370779851290084, φc=-0.0979863038505943, M=2.041340844048836
Iteration 62: ν∞=8.195221418871856e-06, φ∞=-9.203980700036235e-06,νc=-1.2370845413061435, φc=-0.09797894066603427, M=2.041340684178427
Iteration 63: ν∞=7.880916068096636e-06, φ∞=-8.851680870516367e-06,νc=-1.237090846038998, φc=-0.09797185932133785, M=2.041340530499428
Iteration 64: ν∞=7.578752456714613e-06, φ∞=-8.512938809469037e-06,νc=-1.2370969090409634, φc=-0.09796504897029028, M=2.0413403827679697
Iteration 65: ν∞=7.288255026920414e-06, φ∞=-8.18722731398363e-06,νc=-1.237102739644985, φc=-0.09795849918843909, M=2.041340240750164
Iteration 66: ν∞=7.008967320710482e-06, φ∞=-7.874040089215015e-06,νc=-1.2371083468188415, φc=-0.09795219995636771, M=2.041340104221605
Iteration 67: ν∞=6.740451144202162e-06, φ∞=-7.57289089198097e-06,νc=-1.2371137391797569, φc=-0.09794614164365413, M=2.0413399729668886
Iteration 68: ν∞=6.482285946682494e-06, φ∞=-7.28331270002106e-06,νc=-1.2371189250085142, φc=-0.09794031499349411, M=2.0413398467794353
Iteration 69: ν∞=6.2340678744427554e-06, φ∞=-7.00485693670854e-06,νc=-1.2371239122628137, φc=-0.09793471110794474, M=2.041339725460855
Iteration 70: ν∞=5.995409267635799e-06, φ∞=-6.737092710700052e-06,νc=-1.2371287085902278, φc=-0.09792932143377618, M=2.041339608820816
Iteration 71: ν∞=5.7659378766177546e-06, φ∞=-6.479606098637388e-06,νc=-1.2371333213405291, φc=-0.09792413774889727, M=2.0413394966766054
Iteration 72: ν∞=5.545296234532593e-06, φ∞=-6.231999459643535e-06,νc=-1.2371377575775167, φc=-0.09791915214932956, M=2.041339388852786
Iteration 73: ν∞=5.333141095466417e-06, φ∞=-5.993890771495311e-06,νc=-1.237142024090393, φc=-0.09791435703671236, M=2.041339285180965
Iteration 74: ν∞=5.1291428288301275e-06, φ∞=-5.764912999236524e-06,νc=-1.237146127404656, φc=-0.09790974510631298, M=2.0413391854995084
Iteration 75: ν∞=4.932984742153949e-06, φ∞=-5.544713501453002e-06,νc=-1.2371500737924497, φc=-0.09790530933551182, M=2.041339089653033
Iteration 76: ν∞=4.7443628554387236e-06, φ∞=-5.3329534311937865e-06,νc=-1.237153869282734, φc=-0.09790104297276686, M=2.0413389974925855
Iteration 77: ν∞=4.562984975055567e-06, φ∞=-5.129307203908699e-06,νc=-1.237157519670714, φc=-0.09789693952700373, M=2.041338908874885
Iteration 78: ν∞=4.388570614328492e-06, φ∞=-4.933461946101232e-06,νc=-1.2371610305272056, φc=-0.09789299275744685, M=2.041338823662545
Iteration 79: ν∞=4.220850131644209e-06, φ∞=-4.745117006795538e-06,νc=-1.237164407207311, φc=-0.09788919666384141, M=2.0413387417233655
Iteration 80: ν∞=4.0595646862971214e-06, φ∞=-4.5639834487656915e-06,νc=-1.23716765485906, φc=-0.09788554547708239, M=2.0413386629306083
Iteration 81: ν∞=3.904465435509031e-06, φ∞=-4.389783599730176e-06,νc=-1.2371707784314083, φc=-0.09788203365020261, M=2.0413385871623446
Iteration 82: ν∞=3.7553134106793878e-06, φ∞=-4.222250590542499e-06,νc=-1.237173782682137, φc=-0.09787865584973017, M=2.041338514301568
Iteration 83: ν∞=3.611878946718078e-06, φ∞=-4.06112793413437e-06,νc=-1.2371766721852944, φc=-0.09787540694738286, M=2.0413384442358296
Iteration 84: ν∞=3.473941375586355e-06, φ∞=-3.906169111092913e-06,νc=-1.237179451338395, φc=-0.09787228201209398, M=2.0413383768570847
Iteration 85: ν∞=3.3412886401684282e-06, φ∞=-3.7571371767121126e-06,νc=-1.237182124369307, φc=-0.09786927630235262, M=2.041338312061488
Iteration 86: ν∞=3.21371703928988e-06, φ∞=-3.613804375252733e-06,νc=-1.2371846953429384, φc=-0.09786638525885241, M=2.0413382497493724
Iteration 87: ν∞=3.091030754215785e-06, φ∞=-3.475951784773489e-06,νc=-1.2371871681675417, φc=-0.0978636044974246, M=2.0413381898249248
Iteration 88: ν∞=2.9730415707072445e-06, φ∞=-3.3433689692715838e-06,νc=-1.2371895466007983, φc=-0.09786092980224918, M=2.041338132195989
Iteration 89: ν∞=2.8595686859976197e-06, φ∞=-3.2158536376346075e-06,νc=-1.237191834255747, φc=-0.09785835711933907, M=2.041338076774106
Iteration 90: ν∞=2.7504383165676514e-06, φ∞=-3.093211325402048e-06,νc=-1.2371940346064003, φc=-0.09785588255027876, M=2.0413380234742875
Iteration 91: ν∞=2.6454833865440213e-06, φ∞=-2.975255091778895e-06,νc=-1.2371961509931095, φc=-0.09785350234620534, M=2.041337972214794
Iteration 92: ν∞=2.54454334533351e-06, φ∞=-2.8618052198128364e-06,νc=-1.2371981866277857, φc=-0.09785121290202949, M=2.0413379229170903
Iteration 93: ν∞=2.447463854826767e-06, φ∞=-2.7526889367867865e-06,νc=-1.2372001445988696, φc=-0.09784901075088005, M=2.0413378755056564
Iteration 94: ν∞=2.3540966034216442e-06, φ∞=-2.6477401375963293e-06,νc=-1.2372020278761524, φc=-0.09784689255876998, M=2.041337829907959
Iteration 95: ν∞=2.264299045619733e-06, φ∞=-2.5467991230067528e-06,νc=-1.2372038393153888, φc=-0.09784485511947158, M=2.0413377860543602
Iteration 96: ν∞=2.1779340313313536e-06, φ∞=-2.449712361289873e-06,νc=-1.237205581662614, φc=-0.09784289534958254, M=2.0413377438777274
Iteration 97: ν∞=2.094869850417184e-06, φ∞=-2.3563322335668456e-06,νc=-1.2372072575584943, φc=-0.09784101028379569, M=2.0413377033136526
Iteration 98: ν∞=2.014979871735236e-06, φ∞=-2.266516806996268e-06,νc=-1.2372088695423917, φc=-0.09783919707035009, M=2.0413376643002605
Iteration 99: ν∞=1.9381423199190274e-06, φ∞=-2.180129616213063e-06,νc=-1.2372104200562477, φc=-0.09783745296665712, M=2.0413376267780294
Iteration 100: ν∞=1.8642400918841666e-06, φ∞=-2.097039452836699e-06,νc=-1.2372119114483213, φc=-0.09783577533509485, M=2.0413375906896505
Iteration 101: ν∞=1.7931607380595148e-06, φ∞=-2.0171201491717714e-06,νc=-1.2372133459769117, φc=-0.09783416163897551, M=2.041337555980217
Iteration 102: ν∞=1.724796001146764e-06, φ∞=-1.940250396919777e-06,νc=-1.2372147258137125, φc=-0.09783260943865797, M=2.0413375225967827
Iteration 103: ν∞=1.6590419153470534e-06, φ∞=-1.8663135465595915e-06,νc=-1.2372160530472447, φc=-0.09783111638782073, M=2.04133749048857
Iteration 104: ν∞=1.595798414236486e-06, φ∞=-1.7951974373826244e-06,νc=-1.237217329685976, φc=-0.09782968022987082, M=2.0413374596066216
Iteration 105: ν∞=1.5349693884964249e-06, φ∞=-1.7267942133727108e-06,νc=-1.237218557661487, φc=-0.09782829879450013, M=2.0413374299039666
Iteration 106: ν∞=1.476462376852122e-06, φ∞=-1.6610001623236558e-06,νc=-1.2372197388313884, φc=-0.09782696999437027, M=2.0413374013354093
Iteration 107: ν∞=1.4201885454301604e-06, φ∞=-1.5977155487389592e-06,νc=-1.2372208749822247, φc=-0.09782569182193128, M=2.0413373738576057
Iteration 108: ν∞=1.3660623513730204e-06, φ∞=-1.5368444727005418e-06,νc=-1.2372219678321057, φc=-0.09782446234635311, M=2.041337347428709
Iteration 109: ν∞=1.314001659781126e-06, φ∞=-1.4782947082061777e-06,νc=-1.2372230190334335, φc=-0.09782327971058655, M=2.0413373220086
Iteration 110: ν∞=1.263927431061464e-06, φ∞=-1.4219775694802655e-06,νc=-1.2372240301753783, φc=-0.09782214212853096, M=2.041337297558642
Iteration 111: ν∞=1.2157636925095618e-06, φ∞=-1.3678077700497444e-06,νc=-1.2372250027863323, φc=-0.09782104788231492, M=2.0413372740417066
Iteration 112: ν∞=1.1694373972443751e-06, φ∞=-1.315703292285896e-06,νc=-1.23722593833625, φc=-0.0978199953196811, M=2.0413372514221146
Iteration 113: ν∞=1.1248782909366404e-06, φ∞=-1.2655852613235852e-06,νc=-1.237226838238883, φc=-0.09781898285147203, M=2.0413372296655545
Iteration 114: ν∞=1.0820188013340934e-06, φ∞=-1.2173778257776075e-06,νc=-1.2372277038539239, φc=-0.0978180089492114, M=2.0413372087390087
Iteration 115: ν∞=1.0407939823064462e-06, φ∞=-1.1710080380359007e-06,νc=-1.2372285364891098, φc=-0.09781707214278097, M=2.041337188610756
Iteration 116: ν∞=1.0011413599199034e-06, φ∞=-1.1264057446609623e-06,νc=-1.2372293374021976, φc=-0.09781617101818524, M=2.041337169250283
Iteration 117: ν∞=9.63000815368915e-07, φ∞=-1.0835034817459897e-06,νc=-1.23723010780285, φc=-0.09781530421539984, M=2.0413371506281597
Iteration 118: ν∞=9.263146639236584e-07, φ∞=-1.0422363598638722e-06,νc=-1.2372308488545811, φc=-0.09781447042631194, M=2.0413371327162677
Iteration 119: ν∞=8.910273100795878e-07, φ∞=-1.0025419774313603e-06,νc=-1.2372315616764291, φc=-0.09781366839273, M=2.041337115487454
Iteration 120: ν∞=8.570853201037877e-07, φ∞=-9.643603207263804e-07,νc=-1.2372322473446853, φc=-0.09781289690447341, M=2.0413370989156108
From the output of debug mode, we can plot the exterior solution:
[ ]:
# Plot
fig, axs = plt.subplots(2, 2, figsize=(10, 8))
fig.suptitle("Solution from sol_ext", fontsize=14)
# M vs r
axs[0, 0].plot(r_vals_float, M_vals_float)
axs[0, 0].set_xlabel("Radius")
axs[0, 0].set_ylabel("Mass")
axs[0, 0].grid(True)
# nu vs r
axs[0, 1].plot(r_vals_float, nu_vals_float)
axs[0, 1].set_xlabel("Radius")
axs[0, 1].set_ylabel("nu")
axs[0, 1].grid(True)
# phi vs r
axs[1, 0].plot(r_vals_float, phi_vals_float)
axs[1, 0].set_xlabel("Radius")
axs[1, 0].set_ylabel("phi")
axs[1, 0].grid(True)
# psi vs r
axs[1, 1].plot(r_vals_float, psi_vals_float)
axs[1, 1].set_xlabel("Radius")
axs[1, 1].set_ylabel("psi")
axs[1, 1].grid(True)
plt.tight_layout()
plt.show()
And combine them to see the full profile
[ ]:
# Extract interior and exterior profiles
# y0 = (r0, m0, nu0_final, psi0, phi0_final)
r_int_m = sol_iter.ys[0][0]
M_int = sol_iter.ys[1][0]
nu_int = sol_iter.ys[2][0]
psi_int = sol_iter.ys[3][0]
phi_int = sol_iter.ys[4][0]
# y_surf = (M_s, nu_s, phi_s, psi_s)
r_ext_m = sol_ext.ts[0]
M_ext = sol_ext.ys[0][0]
nu_ext = sol_ext.ys[1][0]
phi_ext = sol_ext.ys[2][0]
psi_ext = sol_ext.ys[3][0]
# Convert to kilometers and solar masses
r_int = r_int_m / 1000
r_ext = r_ext_m / 1000
M_int_solar = M_int / utils.solar_mass_in_meter
M_ext_solar = M_ext / utils.solar_mass_in_meter
# Create composite plots with logarithmic x-axis
fig, axs = plt.subplots(2, 2, figsize=(14, 10))
fig.suptitle("Interior + Exterior Solution Profiles", fontsize=16)
# Mass Profile
axs[0, 0].semilogx(r_int, M_int_solar, "b-", linewidth=2, label="Interior")
axs[0, 0].semilogx(r_ext, M_ext_solar, "r--", linewidth=2, label="Exterior")
axs[0, 0].axvline(x=r_int[-1], color="k", linestyle=":", label="Surface")
axs[0, 0].set_xlabel("Radius [km]")
axs[0, 0].set_ylabel(r"Mass [$M_\odot$]")
axs[0, 0].legend()
axs[0, 0].grid(True)
# Metric Function ν
axs[0, 1].semilogx(r_int, nu_int, "b-", linewidth=2, label="Interior")
axs[0, 1].semilogx(r_ext, nu_ext, "r--", linewidth=2, label="Exterior")
axs[0, 1].axvline(x=r_int[-1], color="k", linestyle=":", label="Surface")
axs[0, 1].set_xlabel("Radius [km]")
axs[0, 1].set_ylabel(r"Metric function $\nu$")
axs[0, 1].legend()
axs[0, 1].grid(True)
# Scalar Derivative ψ
axs[1, 0].semilogx(r_int, psi_int, "b-", linewidth=2, label="Interior")
axs[1, 0].semilogx(r_ext, psi_ext, "r--", linewidth=2, label="Exterior")
axs[1, 0].axvline(x=r_int[-1], color="k", linestyle=":", label="Surface")
axs[1, 0].set_xlabel("Radius [km]")
axs[1, 0].set_ylabel(r"Scalar derivative $\psi$")
axs[1, 0].legend()
axs[1, 0].grid(True)
# Scalar Field φ
axs[1, 1].semilogx(r_int, phi_int, "b-", linewidth=2, label="Interior")
axs[1, 1].semilogx(r_ext, phi_ext, "r--", linewidth=2, label="Exterior")
axs[1, 1].axvline(x=r_int[-1], color="k", linestyle=":", label="Surface")
axs[1, 1].set_xlabel("Radius [km]")
axs[1, 1].set_ylabel(r"Scalar field $\phi$")
axs[1, 1].legend()
axs[1, 1].grid(True)
plt.tight_layout(rect=[0, 0, 1, 0.96])
plt.show()
Exploration of scalar-tensor effects#
Below, we show the solutions for a sweep of the beta parameter space.
[ ]:
import numpy as np
# Solve TOV equations for multiple beta values:
beta_values = [0.0, -4.0, -4.5, -5, -5.5, -6]
# initial guess
phi_c = -0.1
nu_c = -1.2
plt.subplots(nrows=1, ncols=2, figsize=(12, 6))
for i, beta_ST in enumerate(beta_values):
# Create EOS with current beta value
eos_tuple = (ns, ps, hs, es, dloge_dlogps, beta_ST, phi_c, nu_c)
logpc, masses, radii, Lambdas = construct_family_ST(
eos_tuple, ndat=50, min_nsat=1.0
)
# Limit masses to be above certain mass
m_min = 0.5
mask = masses > m_min
masses = masses[mask]
radii = radii[mask]
Lambdas = Lambdas[mask]
if i == 0:
c = "black"
zorder = 10
else:
c = None
zorder = 1
# Plot M-R curve for this beta
plt.subplot(121)
plt.plot(radii, masses, label=rf"$\beta={beta_ST:.2f}$", color=c, zorder=zorder)
plt.subplot(122)
# Plot Lambda(M) curve
plt.plot(masses, Lambdas, color=c, zorder=zorder)
plt.subplot(121)
plt.xlabel(r"$R$ [km]")
plt.ylabel(r"$M$ [$M_\odot$]")
plt.legend()
plt.subplot(122)
plt.xlabel(r"$M$ [$M_\odot$]")
plt.ylabel(r"$\Lambda$")
plt.yscale("log")
plt.show()
plt.close()
/var/folders/xj/6sk96vdn15385fpjn9nd1rcw0000gp/T/ipykernel_57185/1370823024.py:46: UserWarning: Data has no positive values, and therefore cannot be log-scaled.
plt.yscale("log")
[ ]: