import numpy as np
import astropy.units as u
import astropy.constants as c
from scipy.optimize import fsolve
from scipy.interpolate import UnivariateSpline, interp1d
import lal
import lalsimulation as lalsim
from .data import PACKAGE_FILENAMES
[docs]
def compute_isco(chi_bh):
'''
This function takes as input the aligned spin component of the BH and
returns the innermost stable circular orbit radius normalized by the mass
of the BH.
'''
if not np.all([np.abs(chi_bh) <= 1.0]):
raise ValueError('|chi1| must be less than 1.0')
Z1 = 1.0 + ((1.0 - chi_bh**2)**(1./3.))*((1 + chi_bh)**(1./3.) +
(1 - chi_bh)**(1./3.))
Z2 = np.sqrt(3*chi_bh**2 + Z1**2)
r_isco = 3 + Z2 - np.sign(chi_bh) * np.sqrt((3 - Z1)*(3 + Z1 + 2*Z2))
return r_isco
def _lalsim_neutron_star_radius(m, max_mass, fam):
'''Return neutron star radius in meters'''
if m > max_mass:
return m * 1500 # 1 solar mass = 1500m
else:
try:
return lalsim.SimNeutronStarRadius(m * lal.MSUN_SI, fam)
except RuntimeError:
# FIXME handle RuntimeError for edge cases
# Raised if m is close to max_mass
# GSL function failed: interpolation error (errnum=1)
# XLAL Error - <GSL function> (interp.c:150): Generic failure
return m * 1500
def _r_ns_from_lal_simulation(m_ns, eosname):
eos = lalsim.SimNeutronStarEOSByName(eosname)
fam = lalsim.CreateSimNeutronStarFamily(eos)
max_mass = lalsim.SimNeutronStarMaximumMass(fam)/c.M_sun.value
try:
iter(m_ns)
R_ns = np.array([_lalsim_neutron_star_radius(m, max_mass, fam)
for m in m_ns])
except TypeError:
R_ns = _lalsim_neutron_star_radius(m_ns, max_mass, fam)
return R_ns, max_mass
def _compactness_baryon_mass(m_ns, r_ns):
C_ns = c.G * m_ns * c.M_sun/(r_ns * u.m * c.c**2)
C_ns = C_ns.value
d1 = 0.619
d2 = 0.1359
BE = (d1*C_ns + d2*C_ns*C_ns) * m_ns # arXiv:1601.06083 (Eq: 21)
m2_b = m_ns + BE # Baryonic mass - Gravitational mass = Binding Energy
return [C_ns, m2_b]
[docs]
def max_mass_from_eosname(eosname):
if eosname == "2H":
max_mass = 2.834648092299807
else:
eos = lalsim.SimNeutronStarEOSByName(eosname)
fam = lalsim.CreateSimNeutronStarFamily(eos)
max_mass = lalsim.SimNeutronStarMaximumMass(fam)/c.M_sun.value
return max_mass
[docs]
def computeCompactness(M_ns, eosname='2H', max_mass=None):
'''
Return the neutron star compactness as a function of mass
and equation of state or radius
Parameters
----------
M_ns : array_like
Neutron star mass in solar masses
eosname : str or interp1d
Neutron star equation of state to be used
max_mass : float
Maximum mass of neutron star.
Returns
-------
[C_ns, m2_b, max_mass]
Compactness, baryon mass and maximum neutron star mass
in solar masses.
Notes
-----
The radius and maximum mass of the neutron star is
inferred based on the equation of state supplied.
Max mass only needs to be supplied for EoS marginalization.
Examples
--------
>>> computeCompactness(2.8)
[array(0.298), array(3.354), 2.834]
>>> computeDiskMass.computeCompactness(2.9, eosname='AP4')
[0.5, 0.0, 2.212]
>>> m_ns = np.array([1.1, 1.2, 1.3])
>>> computeDiskMass.computeCompactness(m_ns, eosname='AP4')
[array([0.141, 0.154, 0.167]), array([1.199, 1.318, 1.439]), 2.212]
'''
if isinstance(eosname, interp1d):
# find R as a function of M
R_ns = eosname(M_ns)
C_ns, m2_b = _compactness_baryon_mass(M_ns, R_ns)
elif eosname != '2H':
# infer radius and maximum mass based on lalsimulation EoS
R_ns, max_mass = _r_ns_from_lal_simulation(M_ns, eosname)
C_ns, m2_b = _compactness_baryon_mass(M_ns, R_ns)
else:
with open(PACKAGE_FILENAMES['equil_2H.dat'], 'rb') as f:
M_g, M_b, Compactness = np.loadtxt(f, unpack=True)
max_mass = max_mass_from_eosname("2H")
s = UnivariateSpline(M_g, Compactness, k=5)
s_b = UnivariateSpline(M_g, M_b, k=5)
C_ns = s(M_ns)
m2_b = s_b(M_ns)
try:
C_ns[M_ns > max_mass] = 0.5 # BH compactness set to 0.5
m2_b[M_ns > max_mass] = 0.0 # BH baryon mass set to 0.0
except TypeError: # if C_ns is not an array
if M_ns > max_mass:
C_ns = 0.5
m2_b = 0.0
return [C_ns, m2_b, max_mass]
'''
For Newtonian: Obtained from arXiv 1807.00011
alpha = 0.406
beta = 0.139
gamma = 0.255
delta = 1.761
For Kerr: Obtained from private message from F. Foucart
alpha = 0.24037857
beta = 0.15184471
gamma = 0.35607169
delta = 1.81491784
'''
[docs]
def xi_implicit(xi, kappa, chi1, eta):
xi_implicit = 3*(xi**2 - 2*kappa*xi + kappa**2*chi1**2)
xi_implicit /= (xi**2 - 3*kappa*xi + 2*chi1*np.sqrt(kappa**3*xi))
xi_implicit -= eta*xi**3
return xi_implicit
[docs]
def computeDiskMass(m1, m2, chi1, chi2, eosname='2H', kerr=False,
R_ns=None, max_mass=None):
"""
This function computes the remnant disk mass after the coalescence using
the equation (4) arXiv 1807.00011.
Parameters
----------
m1, m2 : array_like
primary and secondary mass(es)
chi1, chi2 : array_like
primary and secondary spin(s)
eosname : str
Name of the equation of state to be used. `AP4`
`None` when supplying `R_ns`.
kerr : bool
Supply to use the relativistic tidal parameter. See Fishbone (1971).
R_ns : float
Radius of the secondary in m, assuming it is a neutron star.
max_mass : float
Maximum mass of a neutron star. To be supplied if not supplying EoS.
Example
-------
>>> computeDiskMass(5.0, 2.0, 0.99, 0.)
0.6321412881595185
>>> m1 = np.array([5., 6., 7.])
>>> m2 = np.array([1.0, 1.2, 1.6])
>>> chi1 = np.zeros(3)
>>> chi2 = np.zeros(3)
>>> computeDiskMass(m1, m2, chi1, chi2)
array([0.12833991, 0.05054819, 0.])
>>> computeDiskMass(m1, m2, chi1, chi2, eosname='AP4')
array([0.00851525, 0., 0.])
>>> max_mass = 3.0 # in m_sun
>>> r_ns = 15000. # in meters
>>> computeDiskMass(m1, m2, chi1, chi2,
... R_ns=r_ns, max_mass=max_mass, eosname=None)
array([0.12265712, 0.04272054, 0.])
>>> # m1=2.0, m2=2.0 is a BNS event assuming 2H EOS
>>> # DiskMass == 1.0 is a ad hoc value assigned for BNS events
>>> masses_spins = np.array([2.0, 2.0, 0., 0.])
>>> computeDiskMass(masses_spins[0], masses_spins[1],
... masses_spins[2], masses_spins[3], eosname="2H")
1.0
>>> masses_spins = np.array([5.0, 2.0, 0.99, 0.])
>>> computeDiskMass(masses_spins[0], masses_spins[1],
... masses_spins[2], masses_spins[3], eosname="2H")
0.6321412881595185
Notes
-----
The primary mass, by convention, is larger one. If arrays are supplied,
`m1`, `m2`, `chi1`, `chi1` should be of the same size.
"""
try:
_ = iter(m1)
assert m1.size == m2.size == chi1.size == chi2.size, \
"masses and spins must have the same dimensions"
# Making sure that the supplied m1 >= m2 ##
largerMass1 = m1 >= m2
largerMass2 = m2 > m1
m_primary = np.zeros_like(m1)
m_secondary = np.zeros_like(m2)
chi_primary = np.zeros_like(chi1)
chi_secondary = np.zeros_like(chi2)
m_primary[largerMass1] = m1[largerMass1]
m_primary[~largerMass1] = m2[largerMass2]
m_secondary[~largerMass2] = m2[~largerMass2]
m_secondary[largerMass2] = m1[~largerMass1]
chi_primary[largerMass1] = chi1[largerMass1]
chi_primary[~largerMass1] = chi2[~largerMass1]
chi_secondary[~largerMass2] = chi2[~largerMass2]
chi_secondary[largerMass2] = chi1[~largerMass1]
except TypeError:
if m1 >= m2:
m_primary = m1
m_secondary = m2
chi_primary = chi1
chi_secondary = chi2
else:
m_primary = m2
m_secondary = m1
chi_primary = chi2
chi_secondary = chi1
[C_ns, m2_b, max_mass] = computeCompactness(m_secondary, eosname=eosname,
max_mass=max_mass)
eta = m_primary*m_secondary/(m_primary + m_secondary)**2
BBH = m_secondary > max_mass
BNS = m_primary < max_mass
if not isinstance(BNS, np.ndarray):
if BNS or BBH:
return float(not (BBH) or BNS)
if not kerr:
alpha = 0.406
beta = 0.139
gamma = 0.255
delta = 1.761
firstTerm = alpha*(1 - 2*C_ns)/(eta**(1/3))
secondTerm = beta*compute_isco(chi1)*C_ns/eta
thirdTerm = gamma
else:
alpha = 0.24037857
beta = 0.15184471
gamma = 0.35607169
delta = 1.81491784
kappa = (m_primary/m_secondary)*C_ns
try:
Num = len(kappa)
xi = fsolve(xi_implicit, 100*np.ones(Num), args=(kappa, chi1, eta))
except TypeError:
xi = fsolve(xi_implicit, 100, args=(kappa, chi1, eta))
firstTerm = alpha*xi*(1 - 2*C_ns)
secondTerm = beta*compute_isco(chi1)*C_ns/eta
thirdTerm = gamma
combinedTerms = firstTerm - secondTerm + thirdTerm
if hasattr(combinedTerms, 'ndim') and combinedTerms.ndim >= 1:
lessthanzero = combinedTerms < 0.0
combinedTerms[lessthanzero] = 0.0
else:
combinedTerms = 0 if combinedTerms < 0 else combinedTerms
M_rem = (combinedTerms**delta)*m2_b
if isinstance(BNS, np.ndarray):
M_rem[BBH] = 0.0 # BBH is always EM-Dark
M_rem[BNS] = 1.0 # Ad hoc value, assuming BNS is always EM-Bright
return M_rem