Added metallicity dependence to martinez26 and added martinez26_ecsn

src/cosmic/data/cosmic-settings.json

@@ -123,7 +123,11 @@
                     },
                     {
                         "name": "martinez26",
-                        "description": "Uses sana12 for massive binaries (\\(m_1 ≥ 8 {\\rm M_{\\odot}}\\)) with upper limit at 3000 days and raghavan10 for low-mass binaries. Used in <a href='https://ui.adsabs.harvard.edu/abs/2025arXiv251123285M/abstract'>Martinez+2026</a>."
+                        "description": "Uses <code class='docutils literal notranslate'><span class='pre'>sana12</span></code> for binaries with a primary massive enough to undergo an FeCCSN with an upper limit at 3000 days and <code class='docutils literal notranslate'><span class='pre'>raghavan10</span></code> for low mass binaries. The mass threshold is metallicity dependent, and assumes zsun = 0.02, and standard winds. Used in <a href='https://ui.adsabs.harvard.edu/abs/2025arXiv251123285M/abstract'>Martinez+2026</a>."
+                    },
+                    {
+                        "name": "martinez26_ecsn",
+                        "description": "Same as <code class='docutils literal notranslate'><span class='pre'>martinez26</span></code> but with the <code class='docutils literal notranslate'><span class='pre'>sana12</span></code> power law orbital period distribution for primaries massive enough to undergo an ECSN instead of an FeCCSN."
                     },
                     {
                         "name": "custom",

src/cosmic/sample/sampler/independent.py

@@ -75,8 +75,8 @@ def get_independent_sampler(
         Model to sample eccentricity; choices include: thermal, uniform, sana12
 
     porb_model : `str` or `dict`
-        Model to sample orbital period; choices include: log_uniform, sana12, renzo19, raghavan10, moe19, martinez26
-        or a custom power law distribution defined with a dictionary with keys "min", "max", and "slope"
+        Model to sample orbital period; choices include: log_uniform, sana12, renzo19, raghavan10, moe19, martinez26,
+        martinez26_ecsn or a custom power law distribution defined with a dictionary with keys "min", "max", and "slope"
         (e.g. {"min": 0.15, "max": 0.55, "slope": -0.55}) would reproduce the Sana+2012 distribution
 
     qmin : `float`
@@ -836,11 +836,14 @@ def sample_porb(self, mass1, mass2, rad1, rad2, porb_model, porb_max=None, size=
             `Moe+2019 <https://ui.adsabs.harvard.edu/abs/2019ApJ...875...61M/abstract>_`
             martinez26 : piecewise model with a power law orbital period following
             `Sana+2012 <https://ui.adsabs.harvard.edu/abs/2012Sci...337..444S/abstract>_`
-            between 0.15 < log(P/day) < log(3000) for binaries with m1 >= 8Msun and following
+            between 0.15 < log(P/day) < log(3000) for binaries with primaries large enough to undergo an FeCCSN and following
             `Raghavan+2010 <https://ui.adsabs.harvard.edu/abs/2010ApJS..190....1R/abstract>_`
             with a log normal orbital period in days with mean_logP = 4.9 and sigma_logP = 2.3 between
-            0 < log10(P/day) < 9 for binaries with m1 < 8Msun. Used in
-            `Martinez+2026 <https://ui.adsabs.harvard.edu/abs/2025arXiv251123285M/abstract>_`.
+            0 < log10(P/day) < 9 for binaries with low mass primaries. Used in
+            `Martinez+2026 <https://ui.adsabs.harvard.edu/abs/2025arXiv251123285M/abstract>_`. Assumes zsun = 0.02 for the 
+            metallicity dependence to predict outcomes correctly.
+            martinez26_ecsn : same as martinez26 but with the Sana+2012 power law orbital period distribution
+            for primaries massive enough to undergo an ECSN instead of an FeCCSN.
             Custom power law distribution defined with a dictionary with keys "min", "max", and "slope"
             (e.g. porb_model={"min": 0.15, "max": 0.55, "slope": -0.55}) would reproduce the
             Sana+2012 distribution.
@@ -1041,13 +1044,83 @@ def get_logP_dist(nsamp, norm_wide, norm_close, mu=4.4, sigma=2.1):
             porb = 10**logP_dist 
             aRL_over_a = a_min / utils.a_from_p(porb,mass1,mass2) 
 
-        elif porb_model == "martinez26":
-            # martinez+26 model: use sana12 for mass1 >= 8.0 and raghavan10 for mass1 < 8.0
-            import scipy
+        elif porb_model == "martinez26" or porb_model == "martinez26_ecsn":
+            # martinez26 model: use sana12 for high mass primaries and raghavan10 for low mass primaries. We fit a metallicity dependent minimum mass for an FeCCSN for 'martinez26,
+            # and we fit for the minmum ecsn mass for 'martinez26_ecsn' using data from the population grid from Martinez+2026.
+            try:
+                met = kwargs.pop('met')
+            except:
+                raise ValueError(
+                    "You have chosen martinez26 or martinez26_ecsn for the orbital period distribution which is a metallicity-dependent distribution. "
+                    "Please specify a metallicity for the population."
+                    )
 
+            import scipy
+
             # Create mask for high-mass and low-mass systems
-            (ind_massive,) = np.where(mass1 >= 8.0)
-            (ind_lowmass,) = np.where(mass1 < 8.0)
+            def zams_threshold_mass(metallicity, include_ecsn=False, eps=1e-4):
+                z_by_zsun = np.array([
+                    0.005, 0.0061, 0.0074, 0.009, 0.011, 0.01335, 0.01625, 0.0198,
+                    0.0241, 0.02935, 0.03575, 0.0435, 0.05295, 0.0645, 0.0785,
+                    0.09555, 0.1163, 0.1416, 0.1724, 0.20985, 0.25545, 0.311,
+                    0.3786, 0.4609, 0.56105, 0.683, 0.83145, 1.0, 1.2322, 1.5
+                ])
+                fe_ccsne_low = np.array([
+                    6.800704688282865, 6.800897411951605, 6.800897411951605,
+                    6.801021970762962, 6.8025286354545775, 6.802769553556676,
+                    6.803089427396513, 6.803538527040527, 6.804243674222698,
+                    6.80467206993016, 6.80526072924524, 6.8059692527155,
+                    6.806721013425345, 6.81111108784108, 6.917788936421563,
+                    7.011947874916608, 7.08861989201708, 7.149902100414557,
+                    7.204919186852334, 7.238086587027953, 7.329104218495623,
+                    7.421043416504588, 7.514067557855885, 7.667766797129948,
+                    7.714185711362864, 7.866353261910183, 7.938786132196538,
+                    8.03205271544496, 8.222270760441472, 8.266114315589816
+                ]) - eps
+                incl_ecsne_low = np.array([
+                    6.400693921416194, 6.372129640761526, 6.3238994412467235,
+                    6.309580871018385, 6.318731679747097, 6.351048092438493,
+                    6.414738750162605, 6.44081895134046, 6.4754364194354626,
+                    6.520842448126225, 6.579606889898161, 6.626748726227675,
+                    6.697069204683146, 6.81111108784108, 6.917788936421563,
+                    7.011947874916608, 7.08861989201708, 7.149902100414557,
+                    7.204919186852334, 7.198121695356181, 7.221524796053249,
+                    7.350139621902266, 7.4732604008166295, 7.667766797129948,
+                    7.714185711362864, 7.866353261910183, 7.938786132196538,
+                    8.03205271544496, 8.222270760441472, 8.266114315589816
+                ]) - eps
+
+                masses = incl_ecsne_low if include_ecsn else fe_ccsne_low
+
+                #normalize all metallicities with zsun=0.02
+                metallicity_normalized = metallicity / 0.02
+
+                # take boundary value outside range
+                if metallicity_normalized <= z_by_zsun.min():
+                    return masses[0]
+                if metallicity_normalized >= z_by_zsun.max():
+                    return masses[-1]
+
+                # interpolate in log(Z)
+                logZ = np.log10(z_by_zsun)
+                logZ_target = np.log10(metallicity_normalized)
+
+                # find bracketing indices
+                idx = np.searchsorted(logZ, logZ_target) - 1
+
+                z1, z2 = logZ[idx], logZ[idx + 1]
+                m1, m2 = masses[idx], masses[idx + 1]
+
+                # linear interpolation in logZ
+                frac = (logZ_target - z1) / (z2 - z1)
+                mass_interp = m1 + frac * (m2 - m1)
+
+                return mass_interp
+
+            include_ecsn = porb_model == "martinez26_ecsn"
+            threshold_mass = zams_threshold_mass(met, include_ecsn=include_ecsn)
+            (ind_massive,) = np.where(mass1 >= threshold_mass)
+            (ind_lowmass,) = np.where(mass1 < threshold_mass)
 
             # Initialize porb array
             porb = np.zeros(size)

src/cosmic/tests/test_sample.py

@@ -261,7 +261,7 @@ def test_msort(self):
     def test_sample_porb(self):
         # next do Sana12
         np.random.seed(4)
-        mass1, total_mass = SAMPLECLASS.sample_primary(primary_model='kroupa01', size=100000)
+        mass1, total_mass = SAMPLECLASS.sample_primary(primary_model='kroupa01', size=200000)
         mass2 = SAMPLECLASS.sample_secondary(primary_mass = mass1, qmin=0.1)
         rad1 = SAMPLECLASS.set_reff(mass=mass1, metallicity=0.02)
         rad2 = SAMPLECLASS.set_reff(mass=mass2, metallicity=0.02)
@@ -319,23 +319,47 @@ def test_sample_porb(self):
         self.assertTrue(np.round(log_porb_mean, 1) >= MEAN_RAGHAVAN-0.15)
         self.assertEqual(np.round(log_porb_sigma, 0), np.round(SIGMA_RAGHAVAN, 0))
 
-        # next check martinez26
+        # next check martinez26 for low mass
+        met = 2e-4 # 1/100 solar
+        feccsn_mass, ecsn_mass = 6.8, 6.4
         porb,aRL_over_a = SAMPLECLASS.sample_porb(
-            mass1, mass2, rad1, rad2, 'martinez26', size=mass1.size
+            mass1, mass2, rad1, rad2, 'martinez26', size=mass1.size, met=met
         )
-        # the part of the model with m1 < 8 M_sun should follow the Raghavan10 distribution
-        porb_low_mass = porb[mass1 < 8]
+        # the part of the model with m1 < 6.8 M_sun should follow the Raghavan10 distribution
+        porb_low_mass = porb[mass1 < feccsn_mass]
         log_porb_mean = np.mean(np.log10(porb_low_mass))
         log_porb_sigma = np.std(np.log10(porb_low_mass))
         self.assertTrue(np.round(log_porb_mean, 1) >= MEAN_RAGHAVAN-0.15)
         self.assertEqual(np.round(log_porb_sigma, 0), np.round(SIGMA_RAGHAVAN, 0))
-        # the part of the model with m1 >= 8 M_sun should follow the same power law as Sana12
-        m1_high = mass1+8
-        rad1_high = SAMPLECLASS.set_reff(mass=m1_high, metallicity=0.02)
+
+        # check martinez26_ecsn for low mass
+        porb,aRL_over_a = SAMPLECLASS.sample_porb(
+            mass1, mass2, rad1, rad2, 'martinez26_ecsn', size=mass1.size, met=met
+        )
+        # the part of the model with m1 < 5.3 M_sun should follow the Raghavan10 distribution
+        porb_low_mass = porb[mass1 < ecsn_mass]
+        log_porb_mean = np.mean(np.log10(porb_low_mass))
+        log_porb_sigma = np.std(np.log10(porb_low_mass))
+        self.assertTrue(np.round(log_porb_mean, 1) >= MEAN_RAGHAVAN-0.15)
+        self.assertEqual(np.round(log_porb_sigma, 0), np.round(SIGMA_RAGHAVAN, 0))
+
+        # check martinez_26 for high mass
+        m1_high = mass1+feccsn_mass
+        rad1_high = SAMPLECLASS.set_reff(mass=m1_high, metallicity=met)
+        porb,aRL_over_a = SAMPLECLASS.sample_porb(
+            m1_high, mass2, rad1_high, rad2, 'martinez26', size=m1_high.size, met=met
+        )
+        porb_high_mass = porb[(m1_high >= feccsn_mass) & (np.log10(porb) > 0.5)]
+        power_slope = power_law_fit(np.log10(porb_high_mass), n_bins=25)
+        self.assertEqual(np.round(power_slope, 2), SANA12_PORB_POWER_LAW)
+
+        # check martinez_26_ecsn for high mass
+        m1_high = mass1+ecsn_mass
+        rad1_high = SAMPLECLASS.set_reff(mass=m1_high, metallicity=met)
         porb,aRL_over_a = SAMPLECLASS.sample_porb(
-            m1_high, mass2, rad1_high, rad2, 'martinez26', size=m1_high.size
+            m1_high, mass2, rad1_high, rad2, 'martinez26_ecsn', size=m1_high.size, met=met
         )
-        porb_high_mass = porb[(m1_high >= 8) & (np.log10(porb) > 0.5)]
+        porb_high_mass = porb[(m1_high >= ecsn_mass) & (np.log10(porb) > 0.5)]
         power_slope = power_law_fit(np.log10(porb_high_mass), n_bins=25)
         self.assertEqual(np.round(power_slope, 2), SANA12_PORB_POWER_LAW)