Skip to content

Commit a57d7b1

Browse files
itrharrisonntessore
itrharrison
authored and
ntessore
committed
initial commit of artale models
1 parent 7d1bca4 commit a57d7b1

File tree

3 files changed

+231
-3
lines changed

3 files changed

+231
-3
lines changed

‎skypy/gravitational_waves/__init__.py

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -10,6 +10,7 @@
1010
:toctree: ../api/
1111
1212
b_band_merger_rate
13+
m_star_merger_rate
1314
1415
"""
1516

‎skypy/gravitational_waves/merger_rate.py

Lines changed: 211 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -5,11 +5,13 @@
55
66
"""
77

8+
import numpy as np
89
from astropy import constants, units
910

1011

1112
__all__ = [
1213
'b_band_merger_rate',
14+
'm_star_merger_rate'
1315
]
1416

1517

@@ -29,14 +31,222 @@
2931
'high': 20}
3032
}
3133

34+
artale_tables = {
35+
'NS-NS' : {
36+
'redshift': [0.1, 1.0, 2.0, 6.0],
37+
'alpha1': [1.038, 1.109, 1.050, 1.027],
38+
'alpha1_err': [0.001, 0.001, 0.001, 0.003],
39+
'alpha2': [-6.09, -6.214, -5.533, -5.029],
40+
'alpha2_err': [0.010, 0.006, 0.006, 0.021],
41+
'beta1': [0.800, 0.964, 0.977, 1.113],
42+
'beta1_err': [0.002, 0.001, 0.001, 0.003],
43+
'beta2': [0.323, 0.155, 0.068, -0.070],
44+
'beta2_err': [0.002, 0.001, 0.001, 0.002],
45+
'beta3': [-3.555, -4.819, -4.874, -5.764],
46+
'beta3_err': [0.018, 0.013, 0.011, 0.026],
47+
'gamma1': [0.701, 0.896, 1.018, 1.137],
48+
'gamma1_err': [0.002, 0.002, 0.002, 0.004],
49+
'gamma2': [0.356, 0.184, 0.048, -0.082],
50+
'gamma2_err': [0.002, 0.001, 0.001, 0.002],
51+
'gamma3': [0.411, 0.222, -0.103, -0.053],
52+
'gamma3_err': [0.005, 0.003, 0.002, 0.004],
53+
'gamma4': [-1.968, -3.795, -5.451, -6.104],
54+
'gamma4_err': [0.026, 0.019, 0.017, 0.037],
55+
}
56+
'NS-BH' : {
57+
'redshift': [0.1, 1.0, 2.0, 6.0],
58+
'alpha1': [0.824, 0.873, 0.913, 0.965],
59+
'alpha1_err': [0.001, 0.001, 0.001, 0.002],
60+
'alpha2': [-4.731, -4.478, -4.401, -4.315],
61+
'alpha2_err': [0.008, 0.008, 0.007, 0.018],
62+
'beta1': [0.711, 0.813, 0.871, 0.985],
63+
'beta1_err': [0.002, 0.002, 0.002, 0.003],
64+
'beta2': [0.150, 0.064, 0.039, -0.017],
65+
'beta2_err': [0.002, 0.002, 0.001, 0.001],
66+
'beta3': [-3.536, -3.900, -4.019, -4.490],
67+
'beta3_err': [0.016, 0.018, 0.014, 0.024],
68+
'gamma1': [0.833, 1.074, 1.084, 0.978],
69+
'gamma1_err': [0.002, 0.002, 0.002, 0.003],
70+
'gamma2': [0.101, -0.058, -0.068, -0.013],
71+
'gamma2_err': [0.002, 0.002, 0.001, 0.002],
72+
'gamma3': [-0.461, -0.788, -0.535, 0.016],
73+
'gamma3_err': [0.004, 0.004, 0.003, 0.004],
74+
'gamma4': [-5.434, -7.733, -7.055, -4.386],
75+
'gamma4_err': [0.022, 0.023, 0.018, 0.034],
76+
}
77+
'BH-BH' : {
78+
'redshift': [0.1, 1.0, 2.0, 6.0],
79+
'alpha1': [0.807, 0.813, 0.831, 0.933],
80+
'alpha1_err': [0.001, 0.001, 0.001, 0.004],
81+
'alpha2': [-4.310, -3.845, -3.600, -4.190],
82+
'alpha2_err': [0.006, 0.008, 0.008, 0.026],
83+
'beta1': [0.812, 0.840, 0.858, 1.053],
84+
'beta1_err': [0.001, 0.002, 0.002, 0.004],
85+
'beta2': [-0.006, -0.029, -0.026, -0.098],
86+
'beta2_err': [0.001, 0.002, 0.001, 0.002],
87+
'beta3': [-4.358, -4.109, -3.850, -5.213],
88+
'beta3_err': [0.013, 0.018, 0.015, 0.034],
89+
'gamma1': [0.921, 1.134, 1.135, 1.131],
90+
'gamma1_err': [0.001, 0.002, 0.002, 0.005],
91+
'gamma2': [-0.051, -0.172, -0.167, -0.137],
92+
'gamma2_err': [0.001, 0.002, 0.001, 0.002],
93+
'gamma3': [-0.404, -0.839, -0.681, -0.171],
94+
'gamma3_err': [0.003, 0.004, 0.003, 0.005],
95+
'gamma4': [-6.049, -8.338, -7.758, -6.321],
96+
'gamma4_err': [0.018, 0.024, 0.020, 0.047],
97+
}
98+
}
99+
100+
101+
def m_star_merger_rate(redshift,
102+
m_star,
103+
population):
104+
r"""Model of Artale et al (2020), equation (1) with parameters
105+
from Tables I, II and III.
106+
107+
Compact binary merger rates as a log-log function of a galaxy's
108+
stellar mass.
109+
110+
Parameters are redshift dependent, with linear interpolation
111+
between the simulated points of z={0.1, 1, 2, 6}.
112+
113+
Parameters
114+
----------
115+
redshift : (ngal,) array-like
116+
The redshifts of the galaxies to generate merger
117+
rates for.
118+
m_star : (ngal,) array-like
119+
The stellar mass of the galaxies to generate merger
120+
rates for, in units of stellar mass.
121+
population : {'NS-NS', 'NS-BH', 'BH-BH'}
122+
Compact binary population to get rate for.
123+
'NS-NS' is neutron star - neutron star
124+
'NS-BH' is neutron star - black hole
125+
'BH-BH' is black hole - black hole
126+
127+
Returns
128+
-------
129+
merger_rate : array_like
130+
Merger rates for the galaxies in units of Gigayear^-1
131+
132+
Notes
133+
-----
134+
135+
References
136+
----------
137+
.. Artale et al. 2020, MNRAS,
138+
Volume 491, Issue 3, p.3419-3434 (2020)
139+
https://arxiv.org/abs/1910.04890
140+
141+
Examples
142+
--------
143+
>>> import numpy as np
144+
>>> from skypy.gravitational_waves import m_star_merger_rate
145+
146+
Sample 100 stellar masses values near 10^9 solar masses.
147+
148+
>>> stellar_masses = 10.**(9.0 + np.random.randn(100))
149+
150+
Generate merger rates for these luminosities.
151+
152+
>>> rates = m_star_merger_rate(redshifts,
153+
... stellar_masses,
154+
... population='NS-NS')
155+
156+
"""
157+
158+
alpha1 = np.interp(redshift,
159+
artale_tables[population]['redshift'],
160+
artale_tables[population]['alpha1'])
161+
alpha2 = np.interp(redshift,
162+
artale_tables[population]['redshift'],
163+
artale_tables[population]['alpha2'])
164+
165+
return _m_star_merger_rate(m_star, alpha1, alpha2)
166+
167+
def m_star_sfr_merger_rate(redshift,
168+
m_star,
169+
sfr,
170+
population):
171+
172+
beta1 = np.interp(redshift,
173+
artale_tables[population]['redshift'],
174+
artale_tables[population]['beta1'])
175+
beta2 = np.interp(redshift,
176+
artale_tables[population]['redshift'],
177+
artale_tables[population]['beta2'])
178+
beta3 = np.interp(redshift,
179+
artale_tables[population]['redshift'],
180+
artale_tables[population]['beta3'])
181+
182+
return _m_star_sfr_merger_rate(m_star, sfr, beta1, beta2, beta3)
183+
184+
def m_star_sfr_metallicity_merger_rate(redshift,
185+
m_star,
186+
sfr,
187+
Z,
188+
population):
189+
190+
gamma1 = np.interp(redshift,
191+
artale_table_I[population]['redshift'],
192+
artale_table_I[population]['gamma1'])
193+
gamma2 = np.interp(redshift,
194+
artale_table_I[population]['redshift'],
195+
artale_table_I[population]['gamma2'])
196+
gamma3 = np.interp(redshift,
197+
artale_table_I[population]['redshift'],
198+
artale_table_I[population]['gamma3'])
199+
gamma4 = np.interp(redshift,
200+
artale_table_I[population]['redshift'],
201+
artale_table_I[population]['gamma4'])
202+
203+
return _m_star_sfr_metallicity_merger_rate(m_star, sfr, Z, gamma1, gamma2, gamma3, gamma4)
204+
205+
def _m_star_merger_rate(m_star,
206+
alpha1,
207+
alpha2):
208+
209+
m_star = m_star.to(units.Msun).value
210+
211+
n_gw = 10.**(alpha1 * np.log10(m_star) + alpha2)
212+
213+
return n_gw / units.Gyr
214+
215+
def _m_star_sfr_merger_rate(m_star,
216+
sfr,
217+
beta1,
218+
beta2,
219+
beta3):
220+
221+
m_star = m_star.to(units.Msun).value
222+
sfr = sfr.to(units.Msun / units.year)
223+
224+
n_gw = 10.**(beta1 * np.log10(m_star) + beta2 * np.log10(sfr) + beta3)
225+
226+
return n_gw / units.Gyr
227+
228+
def _m_star_sfr_metallicity_merger_rate(m_star,
229+
sfr,
230+
Z,
231+
gamma1,
232+
gamma2,
233+
gamma3,
234+
gamma4):
235+
236+
m_star = m_star.to(units.Msun).value
237+
sfr = sfr.to(units.Msun / units.year)
238+
239+
n_gw = 10.**(gamma1 * np.log10(m_star) + gamma2 * np.log10(sfr) + gamma3 * np.log10(Z) + gamma4)
240+
241+
return n_gw / units.Gyr
32242

33243
def b_band_merger_rate(luminosity,
34244
population='NS-NS',
35245
optimism='low'):
36246

37247
r"""Model of Abadie et al (2010), Table III
38248
39-
Compact binary merger rates as a linear function of a galaxies
249+
Compact binary merger rates as a linear function of a galaxy's
40250
B-band luminosity.
41251
42252
Parameters

‎skypy/gravitational_waves/tests/test_merger_rate.py

Lines changed: 19 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -1,7 +1,7 @@
11
import numpy as np
22
from astropy import constants, units
3-
from skypy.gravitational_waves import b_band_merger_rate
4-
from skypy.gravitational_waves.merger_rate import abadie_table_III
3+
from skypy.gravitational_waves import b_band_merger_rate, m_star_merger_rate
4+
from skypy.gravitational_waves.merger_rate import abadie_table_III, artale_tables
55

66

77
def test_abadie_rates():
@@ -17,3 +17,20 @@ def test_abadie_rates():
1717
L_B_rate = b_band_merger_rate(L_10, population='NS-NS', optimism='low')
1818
table_value = abadie_table_III['NS-NS']['low']
1919
assert np.isclose(L_B_rate.to_value(1/units.year), table_value, rtol=1e-5)
20+
21+
def test_artale_rates():
22+
23+
# Check the number of merger rates returned
24+
redshifts = np.random.uniform(0., 3., 100)
25+
stellar_masses = 10.**(9.0 + np.random.randn(100))
26+
rates = m_star_merger_rate(redshifts, stellar_masses, population='NS-NS')
27+
assert len(rates) == len(stellar_masses)
28+
29+
# Test that a luminosity of M_sol returns a
30+
# rate that matches the value in Artale Table I
31+
m_sol = constants.Msun.to_value('Msun')
32+
m_sol_rate = m_star_merger_rate(0.0, m_sol, population='NS-NS')
33+
alpha1 = artale_tables['NS-NS']['alpha1'][0]
34+
alpha2 = artale_tables['NS-NS']['alpha2'][0]
35+
table_value = 10.**(alpha1 * np.log10(m_sol) + alpha2)
36+
assert np.isclose(m_sol_rate.to_value(1/units.Gyr), table_value, rtol=1e-5)

0 commit comments

Comments
 (0)