Source code for neost.eos.polytropes

import numpy as np
from math import pow
from scipy.interpolate import UnivariateSpline,interp1d
from scipy.integrate import odeint

from . base import BaseEoS

from .. import global_imports

c = global_imports._c
G = global_imports._G
Msun = global_imports._M_s
pi = global_imports._pi
rho_ns = global_imports._rhons
dyncm2_to_MeVfm3 = global_imports._dyncm2_to_MeVfm3
gcm3_to_MeVfm3 = global_imports._gcm3_to_MeVfm3
oneoverfm_MeV = global_imports._oneoverfm_MeV


[docs] class PolytropicEoS(BaseEoS): """ Class representing a polytropic equation of state object. Parameters ---------- rho_t: float The transition density between the crust EOS and the high density parameterization in cgs. ceft: bool If True a low-density cEFT parameterization is used. ceft_method: str The name of the cEFT calculations used at low density. Can be one of 'Hebeler', 'Drischler', 'Lynn', 'Keller-N2LO', 'Keller-N3L0', or 'Tews'. adm_type: str The name of the ADM particle type. Can be 'None', 'Bosonic', or 'Fermionic' dm_halo: bool If True, ADM halos will be allowed. two_fluid_tidal: bool If True, the two-fluid tidal deformability solver will be used. See TidalDef.py file Methods ------- update(eos_params, max_edsc=True) Update the EoS object with a given set of parameters get_eos_crust() Construct the crust of the equation of state, with or without cEFT. get_eos() Construct the high-density parameterization of the equation of state. eos_core_pp() Function to compute the polytropic equation of state parameterization. add_adm_eos() Function to compute the ADM EoS, whether it be bosonic or fermionic in nature. plot() Plot the equation of state. plot_massradius() Plot the mass-radius curve of the equation of state. """ def __init__(self, crust, rho_t,adm_type = 'None',dm_halo = False,two_fluid_tidal = False): super(PolytropicEoS, self).__init__(crust, rho_t) self.eos_name = 'polytropes' self.param_names = ['gamma1', 'gamma2', 'gamma3', 'rho_t1', 'rho_t2'] self.adm_type = adm_type self.dm_halo = dm_halo self.two_fluid_tidal = two_fluid_tidal if self.ceft is True: self.param_names.append('ceft') if self.adm_type not in ['Bosonic', 'Fermionic', 'None']: raise TypeError('ADM model not recognized. Choose either "Bosonic" or "Fermionic". If ADM is not needed, select "None".') if self.adm_type in ['Bosonic','Fermionic']: self.param_names +=['mchi','gchi_over_mphi', 'adm_fraction']
[docs] def get_eos(self): self.gammas = list(map(self.eos_params.get, ['gamma1', 'gamma2', 'gamma3'])) self.rho_ts = list(map(self.eos_params.get, ['rho_t1', 'rho_t2'])) self._rho_core = np.zeros(297) self._rho_core[0:99] = np.linspace(self.rho_t / rho_ns, self.rho_ts[0], 100)[1::] self._rho_core[99:198] = np.linspace(self.rho_ts[0], self.rho_ts[1], 100)[1::] self._rho_core[198::] = np.logspace(np.log10(self.rho_ts[1]), 2.2, 100)[1::] self._pres_core = np.zeros(len(self._rho_core)) for i, e in enumerate(self._rho_core): self._pres_core[i] = self.eos_core_pp(e, self.P_t * dyncm2_to_MeVfm3) totalrho = np.hstack([self._rho_crust, self._rho_core * rho_ns]) totalpres = np.hstack([self._pres_crust, self._pres_core / dyncm2_to_MeVfm3]) eps0 = self._eds_crust[-1] prho = 0 totalrho,indicies = np.unique(totalrho,return_index = True) totalpres = totalpres[indicies] try: prho = UnivariateSpline(totalrho, totalpres, k=2, s=0) except ValueError: prho = interp1d(totalrho,totalpres,kind = 'linear',fill_value = 'extrapolate') result = odeint(self.edens, eps0, totalrho[totalrho >= self._rho_crust[-1]], args=tuple([prho])) self._eds_core = result.flatten()[1::] totaleps = np.hstack([self._eds_crust, self._eds_core]) self.pressures = totalpres #oringally in cgs units self.energydensities = totaleps self.massdensities = totalrho self.eos = UnivariateSpline(self.energydensities, self.pressures, k=1, s=0) if self.adm_type in ['Bosonic', 'Fermionic']: self.add_adm_eos()
[docs] def eos_core_pp(self, rho, P_t): P_ts, k = (np.zeros(len(self.gammas)) for i in range(2)) P_ts[0] = P_t k[0] = P_t / ((self.rho_t / rho_ns)**self.gammas[0]) P_ts[1] = k[0] * self.rho_ts[0]**self.gammas[0] k[1] = P_ts[1] / (self.rho_ts[0]**self.gammas[1]) P_ts[2] = k[1] * self.rho_ts[1]**self.gammas[1] k[2] = P_ts[2] / (self.rho_ts[1]**self.gammas[2]) self.pres_ts = P_ts[1::] / dyncm2_to_MeVfm3 if rho <= self.rho_ts[0]: pres = k[0] * rho**self.gammas[0] if self.rho_ts[0] < rho <= self.rho_ts[1]: pres = k[1] * rho**self.gammas[1] if self.rho_ts[1] < rho: pres = k[2] * rho**self.gammas[2] return pres
[docs] def polytropic_func(self, rho, K, index): return K * rho**index
[docs] def add_adm_eos(self): def joinarray(ar1,ar2): return np.concatenate((ar1,ar2)) _c = 2.99792458e8 # m/s _hbar = 1.0545718e-34 # J*s MeV_to_Joules = 1.60218e-13 # J self.mchi = self.eos_params['mchi'] * MeV_to_Joules / _c**2. self.gchi_over_mphi = (self.eos_params['gchi_over_mphi'] / MeV_to_Joules * _c**2.) self.adm_fraction = self.eos_params['adm_fraction'] # Bosonic Dark matter component # geometric if self.adm_type == 'Bosonic': #print('Using bosonic EoS') pos_nchi = np.logspace(-5, 55, 5000) self.energydensities_dm = (self.mchi * _c**2 * pos_nchi + 1. / 2. * pow(self.gchi_over_mphi,2) * _hbar**3 * _c**(-1) * pos_nchi**2) *10 / c**2 #Factor of 10 is here to convert from SI base units to cgs base units # /c**2 is here since this is base units of g/(cm s^2) not g/cm^3 self.pressures_dm = (1. / 2. * pow(self.gchi_over_mphi,2) * _hbar**3 * _c**(-1.) * pos_nchi**2) *10 #Factor of 10 is here to convert from SI base units to cgs base units self.massdensities_dm = pos_nchi * self.mchi * pow(10,-3) #Factor of 10^-3 is here to convert from SI base units to cgs base units if self.adm_type == 'Fermionic': #Fermionic Dark matter component # geometric units #print('Using fermionic EoS') pts = 5000 #default point value lower_bound = -2.5 #default lower bound upper_bound = 2.5 #default upper bound if self.eos_params['mchi'] < 10: upper_bound = 6 if self.eos_params['mchi']>=pow(10,3): lower_bound = -14 upper_bound = 1 if self.eos_params['mchi']>= pow(10,5): lower_bound = -10 upper_bound = 1 x= np.logspace(lower_bound,upper_bound,pts) self.energydensities_dm = np.zeros(pts) self.pressures_dm = np.zeros(pts) # APPROXIMATION EoS FOR VALUES OF x < 10**2 z = x[x<pow(10,-2)] interact_term = 1./2.*pow(self.gchi_over_mphi,2)*self.mchi**2*z**6/(3*np.pi**2)**2 self.energydensities_dm[0:len(z)] = (_c**5*self.mchi**4/_hbar**3*(z**3/(3*np.pi**2) + interact_term))*10 / c**2 #Factor of 10 is here to convert from SI base units to cgs base unit # /c**2 is here since this is base units of g/(cm s^2) not g/cm^3 self.pressures_dm[0:len(z)] = (_c**5*self.mchi**4/_hbar**3*(z**5/(15*np.pi**2) + interact_term))*10 #Factor of 10 is here to convert from SI base units to cgs base unit # EXACT EoS FOR VALUES OF x >= 10**-2 y = x[x>=pow(10,-2)] const = (self.mchi**4*_c**5/(8*np.pi**2*_hbar**3)) A= const*(2*y**3+y)*np.sqrt(1+y**2) B = const*np.log(y+np.sqrt(1+y**2)) C = (1/(18*np.pi**4))*pow(self.gchi_over_mphi,2)*(_c**5/_hbar**3)*(self.mchi*y)**6 self.energydensities_dm[len(z)::] = (A-B+C)*10 / c**2 #Factor of 10 is here to convert from SI base units to cgs base unit # /c**2 is here since this is base units of g/(cm s^2) not g/cm^3 self.pressures_dm[len(z)::] = (const*np.sqrt(1+y**2)*((2./ 3.*y**3)-y)+B+C)*10 #Factor of 10 is here to convert from SI base units to cgs base unit #10 ADDED ON 11/18. This is added to convert kg/ms^2 to g/cms^2 nchi = 1./(3*np.pi**2)*(self.mchi*_c*x/_hbar)**3 self.massdensities_dm = nchi * self.mchi * pow(10,-3) #Factor of 10^-3 is here to convert from SI base units to cgs base units
[docs] def check_constraints(self): check = True #required because there are checks in speedofsound.py file, but no constraints are needed for this file if self.reach_fraction == True: check = True else: check = False return check