Skip to content

ENH: gravitational waves artale rates #474

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
wants to merge 23 commits into
base: module/gravitational_waves
Choose a base branch
from
Open

ENH: gravitational waves artale rates #474

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
23 commits
Select commit Hold shift + click to select a range
c5e1007
initial commit of artale models
itrharrison Jul 19, 2021
075265e
tests pass
itrharrison Jul 19, 2021
1bd53ed
fixed doc example
itrharrison Jul 19, 2021
e8437f9
fixed doc example, MSun units
itrharrison Jul 19, 2021
c2ef42e
added sfr and metallicity model docs and tests
itrharrison Jul 19, 2021
8d23946
put back doctest
itrharrison Jul 19, 2021
7877207
fixed sfr docstring
itrharrison Jul 21, 2021
cc259b9
added metallicity function to docs autosummary
itrharrison Jul 21, 2021
5cfcf4e
Update skypy/gravitational_waves/merger_rate.py
itrharrison Jul 22, 2021
de63b1c
Update skypy/gravitational_waves/merger_rate.py
itrharrison Jul 22, 2021
94d16d3
Update skypy/gravitational_waves/tests/test_merger_rate.py
itrharrison Jul 22, 2021
31e69f8
Update skypy/gravitational_waves/merger_rate.py
itrharrison Jul 22, 2021
e8afd1d
Update skypy/gravitational_waves/merger_rate.py
itrharrison Jul 22, 2021
c713d81
Update skypy/gravitational_waves/merger_rate.py
itrharrison Jul 22, 2021
bfac3de
Update skypy/gravitational_waves/merger_rate.py
itrharrison Jul 22, 2021
2464202
Update skypy/gravitational_waves/merger_rate.py
itrharrison Jul 22, 2021
97cab80
Update skypy/gravitational_waves/merger_rate.py
itrharrison Jul 22, 2021
c322ffa
Update skypy/gravitational_waves/merger_rate.py
itrharrison Jul 22, 2021
21b2f19
Update skypy/gravitational_waves/merger_rate.py
itrharrison Jul 22, 2021
15f951d
removed private functions, doc examples
itrharrison Jul 22, 2021
b7e15f9
couple fo small typos
itrharrison Jul 22, 2021
cd58d76
changed docs to specify astropy.Quantity
itrharrison Nov 19, 2021
42b922c
Merge branch 'module/gravitational_waves' into module/gravitational_w…
itrharrison Nov 19, 2021
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
3 changes: 3 additions & 0 deletions skypy/gravitational_waves/__init__.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -10,6 +10,9 @@
:toctree: ../api/

b_band_merger_rate
m_star_merger_rate
m_star_sfr_merger_rate
m_star_sfr_metallicity_merger_rate

"""

Expand Down
285 changes: 267 additions & 18 deletions skypy/gravitational_waves/merger_rate.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -5,11 +5,15 @@

"""

import numpy as np
from astropy import constants, units


__all__ = [
'b_band_merger_rate',
'm_star_merger_rate',
'm_star_sfr_merger_rate',
'm_star_sfr_metallicity_merger_rate'
]


Expand All @@ -29,19 +33,279 @@
'high': 20}
}

artale_tables = {
'NS-NS': {
'redshift': [0.1, 1.0, 2.0, 6.0],
'alpha1': [1.038, 1.109, 1.050, 1.027],
'alpha1_err': [0.001, 0.001, 0.001, 0.003],
'alpha2': [-6.09, -6.214, -5.533, -5.029],
'alpha2_err': [0.010, 0.006, 0.006, 0.021],
'beta1': [0.800, 0.964, 0.977, 1.113],
'beta1_err': [0.002, 0.001, 0.001, 0.003],
'beta2': [0.323, 0.155, 0.068, -0.070],
'beta2_err': [0.002, 0.001, 0.001, 0.002],
'beta3': [-3.555, -4.819, -4.874, -5.764],
'beta3_err': [0.018, 0.013, 0.011, 0.026],
'gamma1': [0.701, 0.896, 1.018, 1.137],
'gamma1_err': [0.002, 0.002, 0.002, 0.004],
'gamma2': [0.356, 0.184, 0.048, -0.082],
'gamma2_err': [0.002, 0.001, 0.001, 0.002],
'gamma3': [0.411, 0.222, -0.103, -0.053],
'gamma3_err': [0.005, 0.003, 0.002, 0.004],
'gamma4': [-1.968, -3.795, -5.451, -6.104],
'gamma4_err': [0.026, 0.019, 0.017, 0.037],
},
'NS-BH': {
'redshift': [0.1, 1.0, 2.0, 6.0],
'alpha1': [0.824, 0.873, 0.913, 0.965],
'alpha1_err': [0.001, 0.001, 0.001, 0.002],
'alpha2': [-4.731, -4.478, -4.401, -4.315],
'alpha2_err': [0.008, 0.008, 0.007, 0.018],
'beta1': [0.711, 0.813, 0.871, 0.985],
'beta1_err': [0.002, 0.002, 0.002, 0.003],
'beta2': [0.150, 0.064, 0.039, -0.017],
'beta2_err': [0.002, 0.002, 0.001, 0.001],
'beta3': [-3.536, -3.900, -4.019, -4.490],
'beta3_err': [0.016, 0.018, 0.014, 0.024],
'gamma1': [0.833, 1.074, 1.084, 0.978],
'gamma1_err': [0.002, 0.002, 0.002, 0.003],
'gamma2': [0.101, -0.058, -0.068, -0.013],
'gamma2_err': [0.002, 0.002, 0.001, 0.002],
'gamma3': [-0.461, -0.788, -0.535, 0.016],
'gamma3_err': [0.004, 0.004, 0.003, 0.004],
'gamma4': [-5.434, -7.733, -7.055, -4.386],
'gamma4_err': [0.022, 0.023, 0.018, 0.034],
},
'BH-BH': {
'redshift': [0.1, 1.0, 2.0, 6.0],
'alpha1': [0.807, 0.813, 0.831, 0.933],
'alpha1_err': [0.001, 0.001, 0.001, 0.004],
'alpha2': [-4.310, -3.845, -3.600, -4.190],
'alpha2_err': [0.006, 0.008, 0.008, 0.026],
'beta1': [0.812, 0.840, 0.858, 1.053],
'beta1_err': [0.001, 0.002, 0.002, 0.004],
'beta2': [-0.006, -0.029, -0.026, -0.098],
'beta2_err': [0.001, 0.002, 0.001, 0.002],
'beta3': [-4.358, -4.109, -3.850, -5.213],
'beta3_err': [0.013, 0.018, 0.015, 0.034],
'gamma1': [0.921, 1.134, 1.135, 1.131],
'gamma1_err': [0.001, 0.002, 0.002, 0.005],
'gamma2': [-0.051, -0.172, -0.167, -0.137],
'gamma2_err': [0.001, 0.002, 0.001, 0.002],
'gamma3': [-0.404, -0.839, -0.681, -0.171],
'gamma3_err': [0.003, 0.004, 0.003, 0.005],
'gamma4': [-6.049, -8.338, -7.758, -6.321],
'gamma4_err': [0.018, 0.024, 0.020, 0.047],
}
}


def m_star_merger_rate(redshift,
m_star,
population):
r"""M_star merger rate.
Model of Artale et al (2020), equation (1) with parameters
from Tables I, II and III.

Compact binary merger rates as a power law function of a galaxy's
stellar mass.

Parameters are redshift dependent, with linear interpolation
between the simulated points of z={0.1, 1, 2, 6}.

Parameters
----------
redshift : (ngal,) array_like
The redshifts of the galaxies to generate merger
rates for.
m_star : (ngal,) astropy.Quantity
The stellar mass of the galaxies to generate merger
rates for, in units of stellar mass.
population : {'NS-NS', 'NS-BH', 'BH-BH'}
Compact binary population to get rate for.
'NS-NS' is neutron star - neutron star
'NS-BH' is neutron star - black hole
'BH-BH' is black hole - black hole

Returns
-------
merger_rate : (ngal,) astropy.Quantity
Merger rates for the galaxies in units of Gigayear^-1

Notes
-----

References
----------
.. Artale et al. 2020, MNRAS,
Volume 491, Issue 3, p.3419-3434 (2020)
https://arxiv.org/abs/1910.04890

"""

alpha1 = np.interp(redshift,
artale_tables[population]['redshift'],
artale_tables[population]['alpha1'])
alpha2 = np.interp(redshift,
artale_tables[population]['redshift'],
artale_tables[population]['alpha2'])

m_star = m_star.to(units.Msun).value

n_gw = 10.**(alpha1 * np.log10(m_star) + alpha2)

return n_gw / units.Gyr


def m_star_sfr_merger_rate(redshift,
m_star,
sfr,
population):
r"""M_star SFR merger rate.
Model of Artale et al (2020), equation (2) with parameters
from Tables I, II and III.

Compact binary merger rates as a power law function of a galaxy's
stellar mass.

Parameters are redshift dependent, with linear interpolation
between the simulated points of z={0.1, 1, 2, 6}.

Parameters
----------
redshift : (ngal,) array_like
The redshifts of the galaxies to generate merger
rates for.
m_star : (ngal,) astropy.Quantity
The stellar mass of the galaxies to generate merger
rates for, in units of stellar mass.
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Same as above

sfr : (ngal,) astropy.Quantity
The star formation rate of the galaxies to generate
merger rates for, in units of stellar mass per year
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

You do a unit transformation for sfr as well. So do we need to specify here?

population : {'NS-NS', 'NS-BH', 'BH-BH'}
Compact binary population to get rate for.
'NS-NS' is neutron star - neutron star
'NS-BH' is neutron star - black hole
'BH-BH' is black hole - black hole

Returns
-------
merger_rate : astropy.Quantity
Merger rates for the galaxies in units of Gigayear^-1

Notes
-----

References
----------
.. Artale et al. 2020, MNRAS,
Volume 491, Issue 3, p.3419-3434 (2020)
https://arxiv.org/abs/1910.04890

"""

beta1 = np.interp(redshift,
artale_tables[population]['redshift'],
artale_tables[population]['beta1'])
beta2 = np.interp(redshift,
artale_tables[population]['redshift'],
artale_tables[population]['beta2'])
beta3 = np.interp(redshift,
artale_tables[population]['redshift'],
artale_tables[population]['beta3'])

m_star = m_star.to(units.Msun).value
sfr = sfr.to(units.Msun / units.year).value

n_gw = 10.**(beta1 * np.log10(m_star) + beta2 * np.log10(sfr) + beta3)

return n_gw / units.Gyr


def m_star_sfr_metallicity_merger_rate(redshift,
m_star,
sfr,
Z,
population):
r"""M_star SFR Metallicity merger rate.
Model of Artale et al (2020), equation (3) with parameters
from Tables I, II and III.

Compact binary merger rates as a power law function of a galaxy's
stellar mass.

Parameters are redshift dependent, with linear interpolation
between the simulated points of z={0.1, 1, 2, 6}.

Parameters
----------
redshift : (ngal,) array_like
The redshifts of the galaxies to generate merger
rates for.
m_star : (ngal,) astropy.Quantity
The stellar mass of the galaxies to generate merger
rates for, in units of stellar mass.
sfr : (ngal,) astropy.Quantity
The star formation rate of the galaxies to generate
merger rates for, in units of stellar mass per year
Z : (ngal,) array_like
The metallicity of the galaxies to generate merger
rates for.
population : {'NS-NS', 'NS-BH', 'BH-BH'}
Compact binary population to get rate for.
'NS-NS' is neutron star - neutron star
'NS-BH' is neutron star - black hole
'BH-BH' is black hole - black hole

Returns
-------
merger_rate : (ngal,) astropy.Quantity
Merger rates for the galaxies in units of Gigayear^-1

Notes
-----

References
----------
.. Artale et al. 2020, MNRAS,
Volume 491, Issue 3, p.3419-3434 (2020)
https://arxiv.org/abs/1910.04890

"""

gamma1 = np.interp(redshift,
artale_tables[population]['redshift'],
artale_tables[population]['gamma1'])
gamma2 = np.interp(redshift,
artale_tables[population]['redshift'],
artale_tables[population]['gamma2'])
gamma3 = np.interp(redshift,
artale_tables[population]['redshift'],
artale_tables[population]['gamma3'])
gamma4 = np.interp(redshift,
artale_tables[population]['redshift'],
artale_tables[population]['gamma4'])

m_star = m_star.to(units.Msun).value
sfr = sfr.to(units.Msun / units.year).value

n_gw = 10.**(gamma1 * np.log10(m_star) + gamma2 * np.log10(sfr) + gamma3 * np.log10(Z) + gamma4)

return n_gw / units.Gyr


def b_band_merger_rate(luminosity,
population='NS-NS',
optimism='low'):

r"""Model of Abadie et al (2010), Table III

Compact binary merger rates as a linear function of a galaxies
Compact binary merger rates as a linear function of a galaxy's
B-band luminosity.

Parameters
----------
luminosity : (ngal,) array-like
luminosity : (ngal,) astropy.Quantity
The B-band luminosity of the galaxies to generate merger
rates for, in units of solar luminosity.
population : {'NS-NS', 'NS-BH', 'BH-BH'}
Expand All @@ -55,7 +319,7 @@ def b_band_merger_rate(luminosity,

Returns
-------
merger_rate : array_like
merger_rate : astropy.Quantity
Merger rates for the galaxies in units of year^-1

Notes
Expand All @@ -67,21 +331,6 @@ def b_band_merger_rate(luminosity,
Volume 27, Issue 17, article id. 173001 (2010)
https://arxiv.org/abs/1003.2480

Examples
--------
>>> import numpy as np
>>> from skypy.gravitational_waves import b_band_merger_rate

Sample 100 luminosity values near absolute magnitude -20.5.

>>> luminosities = 10.**(-0.4*(-20.5 + np.random.randn(100)))

Generate merger rates for these luminosities.

>>> rates = b_band_merger_rate(luminosities,
... population='NS-NS',
... optimism='low')

"""

# Convert luminosity to units of L_10 defined in Abadie et. al. 2010
Expand Down
Loading