astropy · caitlinbegbie · May 22, 2024 · Jun 24, 2024 · Jun 26, 2024 · Jun 26, 2024
diff --git a/changes/BDClusters_AtmIssues.ipynb b/changes/BDClusters_AtmIssues.ipynb
diff --git a/changes/BD_Clusters.ipynb b/changes/BD_Clusters.ipynb
diff --git a/changes/BD_Evolution.ipynb b/changes/BD_Evolution.ipynb
diff --git a/changes/BD_IMF.ipynb b/changes/BD_IMF.ipynb
diff --git a/changes/Exoplanet_IMF.ipynb b/changes/Exoplanet_IMF.ipynb
diff --git a/changes/Test_BD_Cluster.ipynb b/changes/Test_BD_Cluster.ipynb
diff --git a/changes/evolution_gaussian_fit.ipynb b/changes/evolution_gaussian_fit.ipynb
diff --git a/changes/mass_to_age.py b/changes/mass_to_age.py
diff --git a/spisea/atmospheres.py b/spisea/atmospheres.py
diff --git a/spisea/evolution.py b/spisea/evolution.py
diff --git a/spisea/ifmr.py b/spisea/ifmr.py
@@ -36,24 +36,25 @@ def get_Z(self, Fe_H):
         """
         return 10**(Fe_H - 1.85387)
 
-    def Kalirai_mass(self, MZAMS):
+    def Kalirai_mass(self, MZAMS):   # for white dwarfs
         """                                                                                                      
         From Kalirai+08 https://ui.adsabs.harvard.edu/abs/2008ApJ...676..594K/abstract
         1.16 < MZAMS < 6.5
         But we use this function for anything between 0.5 and 9 depending on the IFMR.
         FIXME: need to extend these ranges... explain extension somewhere? Paper maybe?
         """
-
+        
         result = 0.109*MZAMS + 0.394
 
         final = np.zeros(len(MZAMS))
 
-        bad_idx = np.where(MZAMS < 0.5)
+        bad_idx = np.where((MZAMS < 0.5) | (MZAMS >= 9))
         final[bad_idx] = -99
 
-        good_idx = np.where(MZAMS >= 0.5)
+        good_idx = np.where((MZAMS >= 0.5) & (MZAMS < 9))
         final[good_idx] = result[good_idx]
 
+
         return final
 
 
@@ -552,11 +553,12 @@ def generate_death_mass(self, mass_array, metallicity_array):
         * WD: typecode = 101
         * NS: typecode = 102
         * BH: typecode = 103
+        * BD: typecode = 90
         A typecode of value -1 means you're outside the range of 
         validity for applying the ifmr formula. 
         A remnant mass of -99 means you're outside the range of 
         validity for applying the ifmr formula.
-        Range of validity: MZAMS > 0.5 
+        Range of validity: MZAMS > 0.5 and 0.01 < MZAMS < 0.08
 
         Returns
         -------
@@ -569,7 +571,7 @@ def generate_death_mass(self, mass_array, metallicity_array):
         #output_array[1] holds the remnant type
         output_array = np.zeros((2, len(mass_array)))
 
-        codes = {'WD': 101, 'NS': 102, 'BH': 103}
+        codes = {'WD': 101, 'NS': 102, 'BH': 103, 'BD': 90}
 
         #create array to store the remnant masses generated by Spera equations
         rem_mass_array = np.zeros(len(mass_array))
@@ -647,6 +649,12 @@ def generate_death_mass(self, mass_array, metallicity_array):
         #remnant masses of stars with Z > 4.0e-3 and MZAMS > 10.0
         high_metal_high_mass_idx = np.where((Z_array > 4.0e-3) & (core_mass > 10.0))
         rem_mass_array[high_metal_high_mass_idx] = self.M_rem_high_metal_high_mass(Z_array[high_metal_high_mass_idx], core_mass[high_metal_high_mass_idx])
+
+        #classify brown dwarfs before checking for bad indices (in case they are looped into them)
+        BD_idx = np.where((mass_array >= 0.01) & (mass_array < 0.08))
+        rem_mass_array[BD_idx] = mass_array[BD_idx]
+        output_array[0][BD_idx] = rem_mass_array[BD_idx]
+        output_array[1][BD_idx] = codes['BD']
 
         #assign object types based on remnant mass
         bad_idx = np.where(rem_mass_array < 0) #outside the range of validity for the ifmr
@@ -721,7 +729,7 @@ def NS_mass(self, MZAMS):
         """                                                                                                      
         Drawing the NS mass from a Gaussian distrobuton based on observational data.
 
-        Gaussian fit by Emily Ramey and Sergiy Vasylyev of University of Caifornia, Berkeley using a
+        Gaussian fit by Emily Ramey and Sergiy Vasylyev of University of California, Berkeley using a
         sample of NSs from Ozel & Freire (2016) — J1811+2405 Ng et al. (2020), 
         J2302+4442 Kirichenko et al. (2018), J2215+5135 Linares et al. (2018), 
         J1913+1102 Ferdman & Collaboration (2017), J1411+2551 Martinez et al. (2017), 
@@ -751,14 +759,15 @@ def generate_death_mass(self, mass_array):
         * WD: typecode = 101
         * NS: typecode = 102
         * BH: typecode = 103
+        * BD: typecode = 90
 
         A typecode of value -1 means you're outside the range of 
         validity for applying the ifmr formula. 
 
         A remnant mass of -99 means you're outside the range of 
         validity for applying the ifmr formula.
 
-        Range of validity: 0.5 < MZAMS < 120
+        Range of validity: 0.01 < MZAMS < 0.08 and 0.5 < MZAMS < 120
 
         Returns
         -------
@@ -774,23 +783,33 @@ def generate_death_mass(self, mass_array):
         #Random array to get probabilities for what type of object will form
         random_array = np.random.randint(1, 1001, size = len(mass_array))
 
-        codes = {'WD': 101, 'NS': 102, 'BH': 103}
+        codes = {'WD': 101, 'NS': 102, 'BH': 103, 'BD': 90}
 
         """
         The id_arrays are to separate all the different formation regimes
         """
-        id_array0 = np.where((mass_array < 0.5) | (mass_array >= 120))
-        output_array[0][id_array0] = -99 * np.ones(len(id_array0))
-        output_array[1][id_array0]  = -1 * np.ones(len(id_array0))
 
+        #classifying brown dwarfs
+        id_array_BD = np.where((mass_array >= 0.01) & (mass_array < 0.08))
+        output_array[0][id_array_BD] = mass_array[id_array_BD]
+        output_array[1][id_array_BD] = codes['BD']
+
+        #classifying invalid mass ranges
+        id_array0 = np.where((mass_array < 0.01) | ((mass_array >= 0.08) & (mass_array < 0.5)) | (mass_array >= 120))
+        output_array[0][id_array0] = -99
+        output_array[1][id_array0] = -1
+
+        #classifying white dwarfs
         id_array1 = np.where((mass_array >= 0.5) & (mass_array < 9))
         output_array[0][id_array1] = self.Kalirai_mass(mass_array[id_array1])
         output_array[1][id_array1]= codes['WD']
 
+        #classifying neutron stars
         id_array2 = np.where((mass_array >= 9) & (mass_array < 15))
         output_array[0][id_array2] = self.NS_mass(mass_array[id_array2])
         output_array[1][id_array2] = codes['NS']
 
+        #classifying black holes
         id_array3_BH = np.where((mass_array >= 15) & (mass_array < 17.8) & (random_array > 679))
         output_array[0][id_array3_BH] = self.BH_mass_low(mass_array[id_array3_BH], 0.9)
         output_array[1][id_array3_BH] = codes['BH']
@@ -843,6 +862,7 @@ def generate_death_mass(self, mass_array):
         output_array[0][id_array10_NS] = self.NS_mass(mass_array[id_array10_NS])
         output_array[1][id_array10_NS] = codes['NS']
 
+
         return(output_array)
 
 class IFMR_N20_Sukhbold(IFMR):

diff --git a/spisea/imf/imf.py b/spisea/imf/imf.py
@@ -47,7 +47,7 @@ class IMF(object):
         If None, no multiplicity is assumed. Otherwise, use 
         multiplicity object to create multiple star systems.
     """
-    def __init__(self, massLimits=np.array([0.1,150]), multiplicity=None):
+    def __init__(self, massLimits=np.array([0.01,150]), multiplicity=None):
         self._multi_props = multiplicity
         self._mass_limits = massLimits
 
@@ -199,7 +199,10 @@ def generate_cluster(self, totalMass, seed=None):
     def calc_multi(self, newMasses, compMasses, newSystemMasses, newIsMultiple, CSF, MF):
         """
         Helper function to calculate multiples more efficiently.
-        We will use array operations as much as possible
+        We will use array operations as much as possible.
+        Uses Fontanive+18 parameters for brown dwarf masses 
+        (M <= 0.08 M_sun) while keeping default parameters for
+        all other stellar primaries.
         """
         # Identify multiple systems, calculate number of companions for
         # each 
@@ -210,18 +213,37 @@ def calc_multi(self, newMasses, compMasses, newSystemMasses, newIsMultiple, CSF,
             n_comp_arr[too_many] = self._multi_props.CSF_max
         primary = newMasses[idx]
 
+        # limiting BD companions to 1 (mass-based)
+        bd_mask = primary <= 0.08
+        n_comp_arr[bd_mask & (n_comp_arr > 1)] = 1
+
         # We will deal with each number of multiple system independently. This is
         # so we can put in uniform arrays in _multi_props.random_q.
         num = np.unique(n_comp_arr)
         for ii in num:
             tmp = np.where(n_comp_arr == ii)[0]
+            prim_subset = primary[tmp]
+
+            # define masks based on stellar or substellar range
+            bd_mask = prim_subset <= 0.08
+            star_mask = ~bd_mask
 
             if ii == 1:
                 # Single companion case
-                q_values = self._multi_props.random_q(np.random.rand(len(tmp)))
+                q_values = np.empty(len(tmp))
+
+                if np.any(star_mask):
+                    rand_vals = np.random.rand(star_mask.sum())
+                    q_values[star_mask] = self._multi_props.random_q(rand_vals)
+
+                if np.any(bd_mask):
+                    rand_vals = np.random.rand(bd_mask.sum())
+                    b = 1.0 + 6.1   #gamma from Fontanive+18
+                    q_values[bd_mask] = (rand_vals * (1.0 - self._multi_props.q_min ** b) + 
+                                         self._multi_props.q_min ** b) ** (1.0 / b)
 
                 # Calculate mass of companion
-                m_comp = q_values * primary[tmp]
+                m_comp = q_values * prim_subset
 
                 # Only keep companions that are more than the minimum mass. Update
                 # compMasses, newSystemMasses, and newIsMultiple appropriately 
@@ -233,22 +255,30 @@ def calc_multi(self, newMasses, compMasses, newSystemMasses, newIsMultiple, CSF,
                 bad = np.where(m_comp < self._mass_limits[0])[0]
                 newIsMultiple[idx[tmp[bad]]] = False                
             else:
-                # Multple companion case
-                q_values = self._multi_props.random_q(np.random.rand(len(tmp), ii))
-
-                # Calculate masses of companions
-                m_comp = np.multiply(q_values, np.transpose([primary[tmp]]))
-
-                # Update compMasses, newSystemMasses, and newIsMultiple appropriately
+                # Multi-companion case
                 for jj in range(len(tmp)):
-                    m_comp_tmp = m_comp[jj]
+                    prim = prim_subset[jj]
+
+                    # Finding q values of stellar and substellar primaries
+                    if prim <= 0.08:
+                        # BD case (Fontanive+18)
+                        b = 1.0 + 6.1
+                        rand_vals = np.random.rand(ii)
+                        q_values = (rand_vals * (1.0 - self._multi_props.q_min ** b) +
+                                    self._multi_props.q_min ** b) ** (1.0 / b)
+                    else:
+                        # Stellar case (Duchene & Kraus)
+                        q_values = self._multi_props.random_q(np.random.rand(ii))
+
+                    # Calculate masses of companions & update compMasses, newSystemMasses, and newIsMultiple appropriately
+                    m_comp_tmp = q_values * prim
                     compMasses[idx[tmp[jj]]] = m_comp_tmp[m_comp_tmp >= self._mass_limits[0]]
                     newSystemMasses[idx[tmp[jj]]] += compMasses[idx[tmp[jj]]].sum()
-
-                    # Double check for the case when we drop all companions.
-                    # This happens a lot near the minimum allowed mass.
+
+                    # Drop system if no valid companions remain
                     if len(compMasses[idx[tmp[jj]]]) == 0:
                         newIsMultiple[idx[tmp[jj]]] = False
+
 
         return compMasses, newSystemMasses, newIsMultiple
 
@@ -684,12 +714,29 @@ class Weidner_Kroupa_2004(IMF_broken_powerlaw):
     Mass range is 0.01 M_sun - inf M_sun.
     """
     def __init__(self, multiplicity=None):
-        massLimits = np.array([0.01, 0.08, 0.5, 1, np.inf])
+        massLimits = np.array([0.01, 0.08, 0.5, 1, 120])
         powers = np.array([-0.3, -1.3, -2.3, -2.35])
 
         IMF_broken_powerlaw.__init__(self, massLimits, powers,
                                      multiplicity=multiplicity)
-
+
+class Salpeter_Kirkpatrick_2024(IMF_broken_powerlaw):
+    """
+    Define combined IMF from Kirkpatrick (2024) and Salpeter (1955) to allow
+    inclusion of the brown dwarf mass range.
+    Mass range: 
+        * 0.01 M_sun - 8 M_sun: Kirkpatrick 2024
+        <https://ui.adsabs.harvard.edu/abs/2024ApJS..271...55K/abstract>`_.
+        * 8 M_sun - 120 M_sun: Salpeter 1955
+        <https://ui.adsabs.harvard.edu/abs/1955ApJ...121..161S/abstract>`_.
+    """
+    def __init__(self, multiplicity=None):
+        massLimits = np.array([0.01, 0.05, 0.22, 0.55, 8, 120])
+        powers = np.array([-0.6, -0.25, -1.3, -2.3, -2.35])
+
+        IMF_broken_powerlaw.__init__(self, massLimits, powers,
+                                     multiplicity=multiplicity)
+
 ##################################################
 # 
 # Generic functions -- see if we can move these up.

diff --git a/spisea/imf/multiplicity.py b/spisea/imf/multiplicity.py
@@ -41,6 +41,11 @@ class MultiplicityUnresolved(object):
 
                 MF(mass) = MF_amp * (mass ** MF_power)
 
+    However, in the brown dwarf mass regime, it is currently recognized 
+    that only binaries are possible, and the MF decreases dissimilarly 
+    to higher masses (> 0.08 solar masses). The values for this range 
+    are given by Aberasturi et al. (2014) and Fontanive et al. (2023).
+
     **Companion Star Fraction** -- the expected number of companions in
     a multiple system. The companion star fraction (CSF) also 
     changes with mass and this dependency can be described as
@@ -52,6 +57,9 @@ class MultiplicityUnresolved(object):
     value, CSF_max. The actual number of companions is drawn 
     from a Poisson distribution with an expectation value of CSF.
 
+    In the brown dwarf regime we impose an assumption that only
+    binary systems are possible due to current literature trends.
+
     **Mass Ratio (Q)** -- The ratio between the companion star 
     mass and primary star mass, Q = (m_comp / m_prim ) has
     a probability density function described by a powerlaw::
@@ -116,6 +124,9 @@ def multiplicity_fraction(self, mass):
         Given a star's mass, determine the probability that the star is in a
         multiple system (multiplicity fraction = MF).
 
+        Modified to allow binary fraction to decrease in brown dwarf regime.
+        Supported by Aberasturi et al. (2014) and Fontanive et al. (2018).
+
         Parameters
         ----------
         mass : float or numpy array
@@ -133,6 +144,13 @@ def multiplicity_fraction(self, mass):
         if np.isscalar(mf):
             if mf > 1:
                 mf = 1
+            # physically override mf for brown dwarfs
+            if (mass <= 0.08) & (mass > 0.06):
+                mf = 0.16
+            if (mass <= 0.06) & (mass > 0.02):
+                mf = 0.08
+            if (mass < 0.02):
+                mf = 0
         else:
             mf[mf > 1] = 1
 
@@ -141,7 +159,8 @@ def multiplicity_fraction(self, mass):
     def companion_star_fraction(self, mass):
         """
         Given a star's mass, determine the average number of
-        companion stars (companion star fraction = CSF).
+        companion stars (companion star fraction = CSF). For 
+        brown dwarfs we impose a hard limit of one companion.
 
         Parameters
         ----------
@@ -160,8 +179,12 @@ def companion_star_fraction(self, mass):
         if np.isscalar(csf):
             if csf > self.CSF_max:
                 csf = self.CSF_max
+            if (mass <= 0.08):
+                csf = self.multiplicity_fraction(mass)
         else:
             csf[csf > self.CSF_max] = self.CSF_max
+            bd = mass <= 0.08
+            csf[bd] = self.multiplicity_fraction(mass[bd])
 
         return csf
 
@@ -198,6 +221,10 @@ def random_companion_count(self, x, CSF, MF):
         """
         Helper function: calculate number of companions.
         """
+        # bd stipulation since mf=0
+        if MF <= 0:
+            return 0
+
         n_comp = 1 + np.random.poisson((CSF / MF) - 1)
 
         if self.companion_max == True:
@@ -210,6 +237,8 @@ class MultiplicityResolvedDK(MultiplicityUnresolved):
     """
     Sub-class of MultiplicityUnresolved that adds semimajor axis and eccentricity information 
     for multiple objects from distributions described in Duchene and Kraus 2013
+
+    For brown dwarf regime, mean separation and std are given by Fontanive et al. (2018).
 
     Parameters
     --------------
@@ -246,6 +275,8 @@ def log_semimajoraxis(self, mass):
         Generate the semimajor axis for a given mass. The mean and standard deviation of a given mass are determined 
         by fitting the data from fitting the semimajor axis data as a function of mass in table 1 of Duchene and Kraus 2013.
         Then a random semimajor axis is drawn from a log normal distribution with that mean and standard deviation.
+
+        The brown dwarf range is described by mean seperation and std values given in Fontanive et al. (2018).
 
         Parameters
         ----------
@@ -265,6 +296,9 @@ def log_semimajoraxis(self, mass):
             log_a_std = log_a_std_func(np.log10(2.9)) #sigma_log(a)
         if log_a_std < 0.1:
             log_a_std = 0.1
+        if mass <= 0.08:
+            log_a_mean = np.log10(2.9)
+            log_a_std = 0.21
 
         log_semimajoraxis = np.random.normal(log_a_mean, log_a_std)
         while 10**log_semimajoraxis > 2000 or log_semimajoraxis < -2: #AU