From 5f58872aefc443a4f6431999354dfc7f69e655b8 Mon Sep 17 00:00:00 2001 From: Abhishek Date: Fri, 29 May 2026 10:15:32 -0500 Subject: [PATCH 1/9] Added new kick prescriptions --- src/cosmic/src/kick.f | 19 +++++++++++++++++++ 1 file changed, 19 insertions(+) diff --git a/src/cosmic/src/kick.f b/src/cosmic/src/kick.f index 67299594a..d25abab20 100644 --- a/src/cosmic/src/kick.f +++ b/src/cosmic/src/kick.f @@ -281,6 +281,25 @@ SUBROUTINE kick_pfahl(kw,m1,m1c,m1n,m2,ecc,sep,jorb,vk,sn,r2, * Use kick scaling from Bray & Eldridge 2016, Eq. 1 vk = alphakick * ((m1-m1n)/m1n) + betakick vk2 = vk*vk + elseif(abskickflag.eq.7)then +* Asymmetric ejecta / neutrino-driven kick prescription +* Gravitational remnant mass = m1n +* Star mass pre-collapse = m1 +* Ejecta mass = m1 - m1n --> update this +* Baseline values assumed: epsilon_5 = 1.0, f_kin = 0.1, beta_nu = 0.1 + epsilon_5 = 1.0d0 + f_kin = 0.1d0 + beta_nu = 0.1d0 + alpha_ej = 0.01d0 + vk = 21.d0 * SQRT(epsilon_5 * f_kin * beta_nu) * (alpha_ej / 0.01d0) * ((m1 - m1n) / 0.1d0) * (1.5d0 / m1n) + vk2 = vk*vk + elseif(abskickflag.eq.8)then +* From Richardson et al. +* Calibrated against single pulsars, double neutron stars, low kick velocities of ultra stripped supernova + alpha = 115.0d0 + beta = 15.0d0 + vk = alpha * ((m1c - m1n) / m1n) + beta + vk2 = vk*vk endif * If a massless remnant is produced then artificially set the kick to From 998fc8626427c515e5e915c8a794431a4902f27d Mon Sep 17 00:00:00 2001 From: Abhishek Date: Mon, 15 Jun 2026 13:44:21 -0500 Subject: [PATCH 2/9] bug fix --- src/cosmic/src/kick.f | 8 ++++++-- 1 file changed, 6 insertions(+), 2 deletions(-) diff --git a/src/cosmic/src/kick.f b/src/cosmic/src/kick.f index d25abab20..fd6dcfd53 100644 --- a/src/cosmic/src/kick.f +++ b/src/cosmic/src/kick.f @@ -121,6 +121,8 @@ SUBROUTINE kick_pfahl(kw,m1,m1c,m1n,m2,ecc,sep,jorb,vk,sn,r2, real*8 LRL_prev_dot_h, LRL_dot_h_prev, unsigned_psi real*8 disberg_mean real*8 mu_mm + real*8 mean_mns, mean_mej, alphakick, betakick + real*8 epsilon_5, f_kin, beta_nu, alpha_ej, alpha, beta integer i logical ECSN_or_USSN * Output @@ -285,13 +287,15 @@ SUBROUTINE kick_pfahl(kw,m1,m1c,m1n,m2,ecc,sep,jorb,vk,sn,r2, * Asymmetric ejecta / neutrino-driven kick prescription * Gravitational remnant mass = m1n * Star mass pre-collapse = m1 -* Ejecta mass = m1 - m1n --> update this +* Ejecta mass = m1 - m1n --> update this maybe? * Baseline values assumed: epsilon_5 = 1.0, f_kin = 0.1, beta_nu = 0.1 epsilon_5 = 1.0d0 f_kin = 0.1d0 beta_nu = 0.1d0 alpha_ej = 0.01d0 - vk = 21.d0 * SQRT(epsilon_5 * f_kin * beta_nu) * (alpha_ej / 0.01d0) * ((m1 - m1n) / 0.1d0) * (1.5d0 / m1n) + vk = 21.d0 * SQRT(epsilon_5 * f_kin * beta_nu) * + & (alpha_ej / 0.01d0) * ((m1 - m1n) / 0.1d0) * + & (1.5d0 / m1n) vk2 = vk*vk elseif(abskickflag.eq.8)then * From Richardson et al. From 075c16d233fdb6f2533a1e972eea66730992c509 Mon Sep 17 00:00:00 2001 From: Abhishek Date: Wed, 17 Jun 2026 14:32:23 -0500 Subject: [PATCH 3/9] Updated JSON file with kickflags 7 and 8 --- src/cosmic/data/cosmic-settings.json | 8 ++++++++ 1 file changed, 8 insertions(+) diff --git a/src/cosmic/data/cosmic-settings.json b/src/cosmic/data/cosmic-settings.json index 6f7936db8..402653269 100644 --- a/src/cosmic/data/cosmic-settings.json +++ b/src/cosmic/data/cosmic-settings.json @@ -1019,6 +1019,14 @@ "name": 6, "description": "Follows Mandel & Mueller 2020, where the kick velocity is drawn based on the mass of the formed NS or BH, with scaling parameters that are set with 'mm_mu_ns' and 'mm_mu_bh'." }, + { + "name": 7, + "description": "Draws natal kick magnitudes following the gravitational tug-boat mechanism in asymmetric core-collapse explosions from Janka (2017) [https://ui.adsabs.harvard.edu/abs/2017ApJ...837...84J/abstract], with the exact form as in Eq. 2 of Chattaraj et al. (2026) [https://ui.adsabs.harvard.edu/abs/2026ApJ...997...52C/abstract], which is observationally 'calibrated' to the Galactic (field) double neutron star population." + }, + { + "name": 8, + "description": "Follows Richards et al. (2023) [https://ui.adsabs.harvard.edu/abs/2023MNRAS.522.3972R/abstract], which constrains the natal kicks on neutron stars from a combination of observations and the double neutron star merger rate." + }, { "name": "negative values", "description": "Same as above settings but using the old Kiel & Hurley 2009 prescription for changing the orbital configuration of the binary, available for reproducibility purposes but not recommended for new work" From cf5026abcc9adabb76028629e238dce27e6564ac Mon Sep 17 00:00:00 2001 From: Abhishek Date: Wed, 17 Jun 2026 14:55:28 -0500 Subject: [PATCH 4/9] updated hyperlink --- src/cosmic/data/cosmic-settings.json | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/cosmic/data/cosmic-settings.json b/src/cosmic/data/cosmic-settings.json index 402653269..536046ada 100644 --- a/src/cosmic/data/cosmic-settings.json +++ b/src/cosmic/data/cosmic-settings.json @@ -1021,11 +1021,11 @@ }, { "name": 7, - "description": "Draws natal kick magnitudes following the gravitational tug-boat mechanism in asymmetric core-collapse explosions from Janka (2017) [https://ui.adsabs.harvard.edu/abs/2017ApJ...837...84J/abstract], with the exact form as in Eq. 2 of Chattaraj et al. (2026) [https://ui.adsabs.harvard.edu/abs/2026ApJ...997...52C/abstract], which is observationally 'calibrated' to the Galactic (field) double neutron star population." + "description": "Draws natal kick magnitudes following the gravitational tug-boat mechanism in asymmetric core-collapse explosions from Janka (2017), with the exact form as in Eq. 2 of Chattaraj et al. (2026), which is observationally 'calibrated' to the Galactic (field) double neutron star population." }, { "name": 8, - "description": "Follows Richards et al. (2023) [https://ui.adsabs.harvard.edu/abs/2023MNRAS.522.3972R/abstract], which constrains the natal kicks on neutron stars from a combination of observations and the double neutron star merger rate." + "description": "Follows Richards et al. (2023), which constrains the natal kicks on neutron stars from a combination of observations and the double neutron star merger rate." }, { "name": "negative values", From fe80bcbefa6fdd28d07d82369bfb3e915b7ab51f Mon Sep 17 00:00:00 2001 From: Tom Wagg Date: Tue, 23 Jun 2026 17:32:41 -0400 Subject: [PATCH 5/9] allow compilation by avoiding double initialisation --- src/cosmic/src/kick.f | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/src/cosmic/src/kick.f b/src/cosmic/src/kick.f index 93fbef303..798faf63d 100644 --- a/src/cosmic/src/kick.f +++ b/src/cosmic/src/kick.f @@ -122,8 +122,7 @@ SUBROUTINE kick_pfahl(kw,m1,m1c,m1n,m2,ecc,sep,jorb,vk,sn,r2, real*8 LRL_prev_dot_h, LRL_dot_h_prev, unsigned_psi real*8 disberg_mean real*8 mu_mm - real*8 mean_mns, mean_mej, alphakick, betakick - real*8 epsilon_5, f_kin, beta_nu, alpha_ej, alpha, beta + real*8 epsilon_5, f_kin, beta_nu, alpha_ej integer i logical ECSN_or_USSN * Output From bd55e31448ade11de432d6ee02710813ba39b7f4 Mon Sep 17 00:00:00 2001 From: Tom Wagg Date: Tue, 23 Jun 2026 17:34:33 -0400 Subject: [PATCH 6/9] use same structure as Bray, also add 1.4/m1n factor for beta from Eq. 1 --- src/cosmic/src/kick.f | 13 ++++++------- 1 file changed, 6 insertions(+), 7 deletions(-) diff --git a/src/cosmic/src/kick.f b/src/cosmic/src/kick.f index 798faf63d..9099880ad 100644 --- a/src/cosmic/src/kick.f +++ b/src/cosmic/src/kick.f @@ -156,9 +156,6 @@ SUBROUTINE kick_pfahl(kw,m1,m1c,m1n,m2,ecc,sep,jorb,vk,sn,r2, * Set values for mean NS mass and mean ejecta as in Giacobbo & Mapelli 2020 mean_mns = 1.2d0 mean_mej = 9.0d0 -* Set values for alpha and beta as in Bray & Eldridge 2016 - alphakick = 70.0d0 - betakick = 120.0d0 if(using_cmc.eq.0)then * check if we have supplied a randomseed for this SN from kick_info @@ -281,6 +278,8 @@ SUBROUTINE kick_pfahl(kw,m1,m1c,m1n,m2,ecc,sep,jorb,vk,sn,r2, vk2 = vk*vk elseif(abskickflag.eq.4)then * Use kick scaling from Bray & Eldridge 2016, Eq. 1 + alphakick = 70.0d0 + betakick = 120.0d0 vk = alphakick * ((m1-m1n)/m1n) + betakick vk2 = vk*vk elseif(abskickflag.eq.7)then @@ -298,11 +297,11 @@ SUBROUTINE kick_pfahl(kw,m1,m1c,m1n,m2,ecc,sep,jorb,vk,sn,r2, & (1.5d0 / m1n) vk2 = vk*vk elseif(abskickflag.eq.8)then -* From Richardson et al. +* From Richards+2023, Eq 1 (improvement upon Bray & Eldridge 2016) * Calibrated against single pulsars, double neutron stars, low kick velocities of ultra stripped supernova - alpha = 115.0d0 - beta = 15.0d0 - vk = alpha * ((m1c - m1n) / m1n) + beta + alphakick = 115.0d0 + betakick = 15.0d0 + vk = alphakick * ((m1-m1n)/m1n) + betakick * (1.4d0/m1n) vk2 = vk*vk endif From 8e21445eed91dad10db2a0e5065019f6c2e8f42d Mon Sep 17 00:00:00 2001 From: Tom Wagg Date: Tue, 23 Jun 2026 17:36:45 -0400 Subject: [PATCH 7/9] slight clarification in the settings --- src/cosmic/data/cosmic-settings.json | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/cosmic/data/cosmic-settings.json b/src/cosmic/data/cosmic-settings.json index 536046ada..fca181857 100644 --- a/src/cosmic/data/cosmic-settings.json +++ b/src/cosmic/data/cosmic-settings.json @@ -1025,7 +1025,7 @@ }, { "name": 8, - "description": "Follows Richards et al. (2023), which constrains the natal kicks on neutron stars from a combination of observations and the double neutron star merger rate." + "description": "Follows Richards et al. (2023) Eq. 1, which constrains the natal kicks on neutron stars from a combination of observations and the double neutron star merger rate. Similar to kickflag==4, but with different parameters." }, { "name": "negative values", From 6d3f0482f21548157bd7f22989283b3c5428c484 Mon Sep 17 00:00:00 2001 From: Tom Wagg Date: Tue, 23 Jun 2026 22:05:40 -0400 Subject: [PATCH 8/9] first draft of kickflag plot for docs --- docs/settings_examples/plot_kickflag.py | 129 ++++++++++++++++++++++++ 1 file changed, 129 insertions(+) create mode 100644 docs/settings_examples/plot_kickflag.py diff --git a/docs/settings_examples/plot_kickflag.py b/docs/settings_examples/plot_kickflag.py new file mode 100644 index 000000000..07484743e --- /dev/null +++ b/docs/settings_examples/plot_kickflag.py @@ -0,0 +1,129 @@ +""" +Comparing Natal Kick Velocity Prescriptions +=========================================== + +This example demonstrates how changing the `kickflag` parameter can affect the natal kick velocity distributions of compact objects. +We separate this into neutron stars and black holes. +""" + +import numpy as np +import pandas as pd +import time +import matplotlib.pyplot as plt +from cosmic.sample.initialbinarytable import InitialBinaryTable +from cosmic.evolve import Evolve +from cosmic.output import COSMICOutput + + +#---------------------------------------------------------------------------------- +#---------------------------------------------------------------------------------- +## You'll want to edit this part locally to use your own BSEDict and style sheet! + +import sys +sys.path.append("..") +import generate_default_bsedict +BSEDict = generate_default_bsedict.get_default_BSE_settings(to_python=True) + +import matplotlib.pyplot as plt +plt.style.use("../_static/gallery.mplstyle") +#---------------------------------------------------------------------------------- +#---------------------------------------------------------------------------------- + +flag_labels = { + 1: r"kickflag=1: Hobbs+2005 (Maxwellian, $\sigma=265\ \rm km\ s^{-1}$)", + 2: r"kickflag=2: Giacobbo & Mapelli (2020) Eq. 1", + 3: r"kickflag=3: Giacobbo & Mapelli (2020) Eq. 2", + 4: r"kickflag=4: Bray & Eldridge (2016)", + 5: r"kickflag=5: Disberg & Mandel (2025)", + 6: r"kickflag=6: Mandel & Müller (2020)", + 7: r"kickflag=7: Janka (2017) (Eq. 2 in Chattaraj+2026)", + 8: r"kickflag=8: Richards+2023" +} + +BINS = np.linspace(0, 800, 40) + +# sample a population of neutron stars and black holes + +ibt = InitialBinaryTable.sampler( + 'independent', [13, 14], [13, 14], + binfrac_model=1.0, primary_model='kroupa01', + ecc_model='sana12', porb_model='sana12', + qmin=-1, SF_start=13700.0, SF_duration=0.0, + met=0.002, size=2000 +)[0] + +ns_kick_results = {} +bh_kick_results = {} + +# go through each kickflag +for flag in flag_labels: + start_time = time.perf_counter() + print(f"Evolving population with kickflag = {flag}...") + + # update BSEDict with new kickflag + BSEDict['kickflag'] = flag + + # evolve a population and store output in results + bpp, bcm, initC, kick_info = Evolve.evolve( + initialbinarytable=ibt, + BSEDict=BSEDict, + progress=True, + nproc=4 + ) + results = COSMICOutput(bpp=bpp, bcm=bcm, initC=initC, kick_info=kick_info) + + # assign each row in the bpp a unique row_num + results.bpp["row_num"] = np.arange(len(results.bpp)) + + # find the rows that produce the first and second SNe + first_pre_SN_rows = results.bpp[results.bpp["evol_type"] == 15] + second_pre_SN_rows = results.bpp[results.bpp["evol_type"] == 16] + + # convert them to the post SN rows by increasing row_num by one to get the remnant type + first_SN_type = results.bpp[results.bpp["row_num"].isin(first_pre_SN_rows["row_num"] + 1)]["kstar_1"].values + second_SN_type = results.bpp[results.bpp["row_num"].isin(second_pre_SN_rows["row_num"] + 1)]["kstar_2"].values + + # split kicks by the first and second + first_kick = results.kick_info[results.kick_info["star"] == 1]["natal_kick"].values + second_kick = results.kick_info[results.kick_info["star"] == 2]["natal_kick"].values + + # use information on remnant type to separate into NSs and BHs + ns_kick_results[flag] = np.concatenate([first_kick[first_SN_type == 13], second_kick[second_SN_type == 13]]) + bh_kick_results[flag] = np.concatenate([first_kick[first_SN_type == 14], second_kick[second_SN_type == 14]]) + + end_time = time.perf_counter() + elapsed_time = end_time - start_time + print(f" [Time taken: {elapsed_time:.4f} seconds]") + +# plot the results for NSs and BHs separately +for kick_results, title in zip([ns_kick_results, bh_kick_results], + ["Neutron star kicks", "Black hole kicks"]): + + # separate into four flags at a time to clean up the plot + for flags_to_plot, subtitle in zip([[1, 2, 3, 4], [5, 6, 7, 8]], + ["kickflags 1-4", "kickflags 5-8"]): + + fig, ax = plt.subplots() + + for flag in flags_to_plot: + ax.hist(kick_results[flag], bins=BINS, label=flag_labels[flag], density=True, + histtype="step", lw=2, color=f"C{flag - 1}") + + ax.set( + xlabel=r"Natal kick, $v_{\rm kick}$ [$\rm km\ s^{-1}$]", + ylabel=r"Probability density", + title=f"{title} ({subtitle})" + ) + + ax.grid(True, linestyle=":", alpha=0.5, color="gray") + + ax.legend( + loc="upper right", + fontsize=8.5, + facecolor='white', + edgecolor='none', + ) + + plt.tight_layout() + + plt.show() From 8319a7bd95dbbf086cbf13bad92d7f2841a06fa5 Mon Sep 17 00:00:00 2001 From: Tom Wagg Date: Tue, 23 Jun 2026 22:08:29 -0400 Subject: [PATCH 9/9] update title to match others --- docs/settings_examples/plot_kickflag.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/docs/settings_examples/plot_kickflag.py b/docs/settings_examples/plot_kickflag.py index 07484743e..56f5b72c5 100644 --- a/docs/settings_examples/plot_kickflag.py +++ b/docs/settings_examples/plot_kickflag.py @@ -1,6 +1,6 @@ """ -Comparing Natal Kick Velocity Prescriptions -=========================================== +kickflag +======== This example demonstrates how changing the `kickflag` parameter can affect the natal kick velocity distributions of compact objects. We separate this into neutron stars and black holes.