update SIDM parametric model

pyHalo/Halos/HaloModels/NFW_core_trunc.py

@@ -72,38 +72,29 @@ def __init__(self, mass, x, y, r3d, z,
         """
         See documentation in base class (Halos/halo_base.py)
         """
-        # This class must have a cross section model supplied through the kwargs_special argument
-        self._sidm_cross_section = args['SIDM_CROSS_SECTION']
-        _check_valid_cross_section(self._sidm_cross_section)
         super(TNFWCFieldHaloSIDM, self).__init__(mass, x, y, r3d, z,
                  sub_flag, lens_cosmo_instance, args,
                  truncation_class, concentration_class, unique_tag)
 
+    @property
+    def halo_effective_age(self):
+        return self.halo_age
+
     @property
     def sidm_timescale(self):
         """
         Computes the timescale given by Equation 2.2 in https://arxiv.org/pdf/2305.16176.pdf
         :return:
         """
-        if not hasattr(self, '_sidm_timescale'):
-            rho_s, rs, _ = self.nfw_params
-            C = 0.75
-            G = 4.3e-6
-            conversion_fac = 2.135788e-10
-            vmax = 1.65 * np.sqrt(G * rho_s * rs ** 2) # fixed for an NFW profile
-            v_scale = np.sqrt(4 * np.pi * G * rho_s * rs ** 2)
-            sigma_eff = self._sidm_cross_section.effective_cross_section(0.64 * vmax)
-            denom = sigma_eff * rho_s * v_scale * conversion_fac
-            self._sidm_timescale = self._scale_evolution_timescale * 150 / C / denom
-        return self._sidm_timescale
+        return self._args['sidm_timescale']
 
     @property
     def t_over_tc(self):
         """
         Computes the dimensionless timescale for the halo evolution
         :return:
         """
-        return self.halo_age / self.sidm_timescale
+        return self.halo_effective_age / self.sidm_timescale
 
     @property
     def lenstronomy_params(self):
@@ -130,22 +121,61 @@ def lenstronomy_params(self):
 
     def profile_args(self):
         """
-        TODO: the parameterization for the density profile presented by Yang et al. doesn't conserve mass.
         I temporarily fix this by rescaling the normalization to conserve mass at all times
 
         Computes the time-evolving profile parameters
         :return: the density normalization, scale radius, and core radius for thee SIDM halo
         """
         if not hasattr(self, '_profile_args'):
-            t = min(self._regularize_t_over_tc, self.t_over_tc)
-            rhos_0, rs_0, r200_0 = self.nfw_params
-            rho_s = rhos_0 * rho_s_evolution(t)
-            rs_kpc = rs_0 * rs_evolution(t)
-            rc_kpc = rs_0 * rc_evolution(t)
+            t = self.t_over_tc
             rt_kpc = self._truncation_class.truncation_radius_halo(self)
-            self._profile_args = (1.0 * rho_s, rs_kpc, rc_kpc, rt_kpc, r200_0)
+            rho_s, rs_kpc, rc_kpc = self.get_params(t)
+            r200_0 = self.nfw_params[-1]
+            self._profile_args = (rho_s, rs_kpc, rc_kpc, rt_kpc, r200_0)
         return self._profile_args
 
+    def get_params(self, t_over_tc):
+        """
+
+        :param t:
+        :return:
+        """
+        rhos_0, rs_0, r200_0 = self.nfw_params
+        t_cut = 1.001
+        t_max = 1.7
+        t_over_tc = min(t_max, t_over_tc)
+        t_min = 1e-2
+        if t_over_tc < t_min:
+            rho_s = rhos_0
+            rs_kpc = rs_0
+            rc_kpc = 0.0
+        elif t_over_tc < t_cut:
+            rho_s = rhos_0 * rho_s_evolution(t_over_tc)
+            rs_kpc = rs_0 * rs_evolution(t_over_tc)
+            rc_kpc = rs_0 * rc_evolution(t_over_tc)
+        else:
+            rhos_0, rs_0, r200_0 = self.nfw_params
+            rhos_last = rhos_0 * rho_s_evolution(t_cut)
+            rs_kpc_last = rs_0 * rs_evolution(t_cut)
+            rc_kpc_last = rs_0 * rc_evolution(t_cut)
+            profile_args = (rhos_last, rs_kpc_last, rc_kpc_last, 100 * r200_0, r200_0)
+            total_mass = self.mass_3d(r200_0, profile_args)
+            rc_kpc = rs_0 * rc_extrapolation(t_over_tc)
+            rs_kpc = rs_0 * rs_extrapolation(t_over_tc)
+            rho_s = self.rhos_extrapolation(total_mass, rs_kpc, rc_kpc)
+
+        return rho_s, rs_kpc, rc_kpc
+
+    def rhos_extrapolation(self, total_mass, rs_kpc, rc_kpc):
+        """
+
+        :return:
+        """
+        r200_0 = self.nfw_params[-1]
+        profile_args_integral = (1.0, rs_kpc, rc_kpc, 100 * r200_0, r200_0)
+        m_integral = self.mass_3d(r200_0, profile_args_integral)
+        return total_mass / m_integral
+
     @property
     def params_physical(self):
         """
@@ -164,8 +194,6 @@ class TNFWCSubhaloSIDM(TNFWCFieldHaloSIDM):
     """
     Implements a temporal evolution of the halo profile based on Yang et al. (2023)
     https://arxiv.org/pdf/2305.16176.pdf
-
-    TODO: add something here to account for tidal effects
     """
     def __init__(self, mass, x, y, r3d, z,
                  sub_flag, lens_cosmo_instance, args,
@@ -180,6 +208,12 @@ def __init__(self, mass, x, y, r3d, z,
                  sub_flag, lens_cosmo_instance, args,
                  truncation_class, concentration_class, unique_tag)
 
+    @property
+    def halo_effective_age(self):
+        if not hasattr(self, '_halo_effective_age'):
+            self._halo_effective_age = self.lens_cosmo.sidm_halo_effective_age(self.z, self.z_infall, self._args['lambda_t'])
+        return self._halo_effective_age
+
 def _check_valid_cross_section(cross_section):
     """
 
@@ -222,3 +256,24 @@ def rc_evolution(tr):
     """
 
     return 1.06419 * np.sqrt(tr) - 1.33207 * tr + 0.431133 * tr ** 2 - 0.148087 * tr ** 3 + 0.024455 * tr ** 4
+
+def rs_extrapolation(tr, t0=0.5, c1=0.47, c2=2.4):
+    """
+
+    :param tr:
+    :return:
+    """
+    return rs_evolution(t0) * np.exp(-c1 * (tr - t0) - c2 * (tr - t0) ** 3)
+
+def rc_extrapolation(t_over_tc, t0=0.5, c1=0.28):
+    """
+
+    :param t_over_tc:
+    :param t0:
+    :param c1:
+    :return:
+    """
+    rc = rc_evolution(t0) - c1 * (t_over_tc - t0)
+    if rc < 0.00001:
+        rc = 0.00001
+    return rc

pyHalo/Halos/HaloModels/powerlaw.py

@@ -83,8 +83,10 @@ def profile_args(self):
         See documentation in base class (Halos/halo_base.py)
         """
         if not hasattr(self, '_profile_args'):
-
-            concentration = self._concentration_class.nfw_concentration(self.mass, self.z_eval)
+            if 'c' not in list(self._args.keys()):
+                concentration = self._concentration_class.nfw_concentration(self.mass, self.z_eval)
+            else:
+                concentration = self._args['c']
             gamma = self._args['log_slope_halo']
             x_core_halo = self._args['x_core_halo']
             self._profile_args = (concentration, gamma, x_core_halo)

pyHalo/Halos/concentration.py

@@ -22,20 +22,24 @@ def __init__(self, cosmo, scatter=True, scatter_dex=0.2, *args, **kwargs):
         self._scatter = scatter
         self._scatter_dex = scatter_dex
 
-    def nfw_concentration(self, m, z):
+    def nfw_concentration(self, m, z, force_no_scatter=False):
         """
         Evaluates the concentration of a halo of mass 'm' at redshift z
         :param M: halo mass [M_sun]
         :param z: halo redshift
+        :param force_no_scatter: bool; if True, will return the median concentation
         :return:
         """
         if isinstance(m, float) or isinstance(m, int):
             c = float(self._evaluate_concentration(m, z))
         else:
             c = numpy.array([self._evaluate_concentration(mi, z) for mi in m])
         if self._scatter:
-            log_c = numpy.log(c)
-            c = numpy.random.lognormal(log_c, self._scatter_dex)
+            if force_no_scatter:
+                pass
+            else:
+                log_c = numpy.log(c)
+                c = numpy.random.lognormal(log_c, self._scatter_dex)
         if isinstance(c, float):
             c = max(c, self._universal_minimum)
         else:

pyHalo/Halos/halo_base.py

@@ -215,3 +215,15 @@ def mass_3d(self, rmax, profile_args=None):
             rmax = self.nfw_params[2]
         _integrand = lambda r: 4 * np.pi * r ** 2 * self.density_profile_3d(r, profile_args)
         return quad(_integrand, 0, rmax)[0]
+
+    def logarithmic_profile_slope(self, r, profile_args=None):
+        """
+
+        :param r:
+        :param profile_args:
+        :return:
+        """
+        density = self.density_profile_3d(r, profile_args)
+        log_density = np.log(density)
+        logr = np.log(r)
+        return np.gradient(log_density, logr)

pyHalo/Halos/lens_cosmo.py

@@ -328,3 +328,19 @@ def halo_dynamical_time(self, m_host, z, c_host):
         rho_average = m_host / volume
         g = 4.3e-6
         return 0.5427 / numpy.sqrt(g*rho_average)
+
+    def sidm_halo_effective_age(self, z, z_infall, lambda_t, zform=10.0):
+        """
+        Calculates a time since z = zform t_1 + t_2 where t_1 is the time from formation to infall, and t_2
+        is the time from infall to redshift z times lambda_t
+        :param z: halo redshift at the time of lensing
+        :param z_infall: infall redshift
+        :param lambda_t: rescales the passage of time since the halo becomes a subhalo
+        :param zform: formation redshift
+        :return: "age" in Gyr
+        """
+        if z_infall > 10:
+            z_infall = 10
+        time_formation_to_infall = self.cosmo.halo_age(z_infall, zform=zform)
+        time_infall_to_z = self.cosmo.halo_age(z, zform=z_infall)
+        return time_formation_to_infall + lambda_t * time_infall_to_z

pyHalo/realization_extensions.py

@@ -3,6 +3,7 @@
 from pyHalo.Halos.HaloModels.generalized_nfw import GeneralNFWSubhalo, GeneralNFWFieldHalo
 from pyHalo.single_realization import Realization
 from pyHalo.Halos.HaloModels.gaussianhalo import GaussianHalo
+from pyHalo.concentration_models import preset_concentration_models
 from pyHalo.Rendering.correlated_structure import CorrelatedStructure
 from pyHalo.Rendering.MassFunctions.delta_function import DeltaFunction
 from pyHalo.Rendering.MassFunctions.gaussian import Gaussian
@@ -40,7 +41,7 @@ def add_globular_clusters(self, log10_mgc_mean, log10_mgc_sigma, rendering_radiu
         :param log10_mgc_mean: the median log10(mass) of the GC's
         :param log10_mgc_sigma: the standard deviation of the Gaussian mass function for log10(m)
         :param rendering_radius_arcsec [arcsec]: sets the area around (center_x, center_y) where the GC's appear; GC's are
-        distributed uniformly in a circule centered at (center_x, center_y) with this radius
+        distributed uniformly in a circle centered at (center_x, center_y) with this radius
         :param gc_profile_args: the keyword arguments for the GC mass profile, must include "gamma", "gc_size_lightyear",
         "r_core_fraction", or the logarithmic power law slope, the size of the gc in light years, and the size of the core
         relative to the size of the GC
@@ -51,7 +52,6 @@ def add_globular_clusters(self, log10_mgc_mean, log10_mgc_sigma, rendering_radiu
         :param coordinate_center_y: center of rendering area in arcsec
         :return: an instance of Realization that includes the GC's
         """
-
         n = self.number_globular_clusters(log10_mgc_mean, log10_mgc_sigma, rendering_radius_arcsec, galaxy_Re,
                                           host_halo_mass, f)
         mfunc = Gaussian(n, log10_mgc_mean, log10_mgc_sigma)
@@ -79,7 +79,7 @@ def add_globular_clusters(self, log10_mgc_mean, log10_mgc_sigma, rendering_radiu
     def number_globular_clusters(self, log10_mgc_mean, log10_mgc_sigma, rendering_radius_arcsec,
                               galaxy_Re=10.0, host_halo_mass=10**13.3, f=3.4e-5):
         """
-
+        Compute the number of globular clusters from their mass function (see documentation in add_globular_clusters)
         :return:
         """
         m_gc = f * host_halo_mass
@@ -93,34 +93,62 @@ def number_globular_clusters(self, log10_mgc_mean, log10_mgc_sigma, rendering_ra
         n = int(mass_in_gc / integral)
         return n
 
-    def toSIDM(self, cross_section):
+    def toSIDM_single_halo(self, halo, t_c, subhalo_evolution_scaling):
+        """
+        Transform a single NFW profile to an SIDM profile based on the halo age, its status as a subhalo or a field halo,
+        and the core collapse timescale
+        :param halo: an instance of a Halo model
+        :param t_c: SIDM collapse timescale in Gyr
+        :param median_mc_relation: the concentration-mass relation used to create the CDM realization with scatter=False
+        :param mass_bin_center:
+        :return : an SIDM halo derived from the original profile
+        """
+        from pyHalo.Halos.HaloModels.NFW_core_trunc import TNFWCFieldHaloSIDM, TNFWCSubhaloSIDM
+
+        if halo.is_subhalo:
+
+            subhalo_flag = True
+            kwargs_profile = {'lambda_t': subhalo_evolution_scaling, 'sidm_timescale': t_c}
+            new_halo = TNFWCSubhaloSIDM(halo.mass, halo.x, halo.y, halo.r3d, halo.z, subhalo_flag,
+                                        halo.lens_cosmo, kwargs_profile,
+                                        halo._truncation_class, halo._concentration_class,
+                                        halo.unique_tag)
+            new_halo._z_infall = halo.z_infall
+        else:
+
+            subhalo_flag = False
+            kwargs_profile = {'lambda_t': 1.0, 'sidm_timescale': t_c}
+            new_halo = TNFWCFieldHaloSIDM(halo.mass, halo.x, halo.y, halo.r3d, halo.z, subhalo_flag,
+                                          halo.lens_cosmo, kwargs_profile,
+                                          halo._truncation_class, halo._concentration_class,
+                                          halo.unique_tag)
+        return new_halo
+
+    def toSIDM(self, mass_bin_list, core_collapse_timescale_list, subhalo_evolution_scaling):
         """
+        This takes a CDM relization and transforms it into an SIDM realization. The density profile follows
+        https://arxiv.org/pdf/2305.16176.pdf if t / t_c <= 1. For t / t_c > 1 we extrapolate
 
-        :param cross_section: a class with a method effective_cross_section(v) that takes as input a velocity scale
-        equal to 0.64 * v_max and returns an effective cross section in units cm^2 / gram
-        :return: a realization where all halos aree transformed into SIDM halos using the parameteric model by
-        Yang et al. (2023) (https://arxiv.org/pdf/2305.16176.pdf)
+        :param mass_bin_list: a list with len(core_collapse_timescale_list) of lists that specify log10mass ranges
+        :param core_collapse_timescale_list: the core collapse timescale in each mass bin in log10(Gyr)
+        :return:
         """
 
-        from pyHalo.Halos.HaloModels.NFW_core_trunc import TNFWCFieldHaloSIDM, TNFWCSubhaloSIDM
         sidm_halos = []
-        kwargs_profile = {'SIDM_CROSS_SECTION': cross_section}
         for halo in self._realization.halos:
-            if halo.is_subhalo:
-                new_halo = TNFWCSubhaloSIDM(halo.mass, halo.x, halo.y, halo.r3d, halo.z, False,
-                                              halo.lens_cosmo, kwargs_profile,
-                                              halo._truncation_class, halo._concentration_class,
-                                              halo.unique_tag)
-                new_halo._rperi_units_r200 = halo.rperi_units_r200
-                new_halo._time_since_infall = halo.time_since_infall
+            for bin_index, mass_bin in enumerate(mass_bin_list):
+                if np.log10(halo.mass) >= mass_bin[0] and np.log10(halo.mass) < mass_bin[1]:
+                    t_c = 10 ** core_collapse_timescale_list[bin_index]
+                    break
             else:
-                new_halo = TNFWCFieldHaloSIDM(halo.mass, halo.x, halo.y, halo.r3d, halo.z, False,
-                                              halo.lens_cosmo, kwargs_profile,
-                                              halo._truncation_class, halo._concentration_class,
-                                              halo.unique_tag)
-            new_halo._c = halo.c
-            new_halo._zeval = halo.z_eval
-            new_halo._halo_age = halo.halo_age
+                raise Exception('halo mass '+str(np.log10(halo.mass))+' not inside the minimum/maximum mass ranges')
+            c_halo = halo.c
+            c_median = halo._concentration_class.nfw_concentration(halo.mass,
+                                                                   halo.z_eval,
+                                                                   force_no_scatter=True)
+            delta_c = c_halo / c_median
+            concentration_factor = delta_c ** (7/2)
+            new_halo = self.toSIDM_single_halo(halo, concentration_factor * t_c, subhalo_evolution_scaling)
             sidm_halos.append(new_halo)
 
         new_realization = Realization.from_halos(sidm_halos, self._realization.lens_cosmo, self._realization.kwargs_halo_model,

pyHalo/single_realization.py

@@ -104,7 +104,6 @@ def __init__(self, masses, x, y, r3d, mdefs, z, subhalo_flag, lens_cosmo,
         if halos is None:
             for mi, xi, yi, r3di, mdefi, zi, sub_flag in zip(masses, x, y, r3d,
                            mdefs, z, subhalo_flag):
-
                 unique_tag = np.random.rand()
                 model = self._load_halo_model(mi, xi, yi, r3di, mdefi, zi, sub_flag, self.lens_cosmo,
                                               kwargs_halo_model, unique_tag)