Added new kick prescriptions

docs/settings_examples/plot_kickflag.py

@@ -0,0 +1,129 @@
+"""
+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.
+"""
+
+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()

src/cosmic/data/cosmic-settings.json

@@ -1019,6 +1019,14 @@
                         "name": 6,
                         "description": "Follows <a class='reference external' href='https://scixplorer.org/abs/2020MNRAS.499.3214M/abstract'>Mandel &amp; Mueller 2020</a>, 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 <a class='reference external' href='https://ui.adsabs.harvard.edu/abs/2017ApJ...837...84J/abstract'>Janka (2017)</a>, with the exact form as in Eq. 2 of <a class='reference external' href='https://ui.adsabs.harvard.edu/abs/2026ApJ...997...52C/abstract'>Chattaraj et al. (2026)</a>, which is observationally 'calibrated' to the Galactic (field) double neutron star population."
+                    },
+                    {
+                        "name": 8,
+                        "description": "Follows <a class='reference external' href='https://ui.adsabs.harvard.edu/abs/2023MNRAS.522.3972R/abstract'>Richards et al. (2023)</a> 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",
                         "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"

src/cosmic/src/kick.f

@@ -122,6 +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 epsilon_5, f_kin, beta_nu, alpha_ej
       integer i
       logical ECSN_or_USSN
 * Output
@@ -155,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
@@ -280,8 +278,31 @@ 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
+* Asymmetric ejecta / neutrino-driven kick prescription
+* Gravitational remnant mass = m1n
+* Star mass pre-collapse = m1
+* 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)
+             vk2 = vk*vk
+          elseif(abskickflag.eq.8)then
+* From Richards+2023, Eq 1 (improvement upon Bray & Eldridge 2016)
+* Calibrated against single pulsars, double neutron stars, low kick velocities of ultra stripped supernova
+             alphakick = 115.0d0
+             betakick = 15.0d0
+             vk = alphakick * ((m1-m1n)/m1n) + betakick * (1.4d0/m1n)
+             vk2 = vk*vk
           endif
 
 * If a massless remnant is produced then artificially set the kick to