Merge pull request #65 from dangilman/sidm_evolution

Sidm evolution

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/PresetModels/sidm.py

@@ -83,3 +83,72 @@ class for details
         probabilities_subhalos, probabilities_field_halos, kwargs_sub_function, kwargs_field_function)
     sidm = extension.add_core_collapsed_halos(index_collapsed, collapsed_halo_profile, **kwargs_collapsed_profile)
     return sidm
+
+def SIDM_parametric(z_lens, z_source, log10_mass_ranges, collapse_timescales, subhalo_time_scaling=1.0,
+        sigma_sub=0.025, log10_sigma_sub=None, log_mlow=6., log_mhigh=10.,
+        concentration_model_subhalos='DIEMERJOYCE19', kwargs_concentration_model_subhalos={},
+        concentration_model_fieldhalos='DIEMERJOYCE19', kwargs_concentration_model_fieldhalos={},
+        truncation_model_subhalos='TRUNCATION_GALACTICUS', kwargs_truncation_model_subhalos={},
+        infall_redshift_model='HYBRID_INFALL', kwargs_infall_model={},
+        truncation_model_fieldhalos='TRUNCATION_RN', kwargs_truncation_model_fieldhalos={},
+        shmf_log_slope=-1.9, cone_opening_angle_arcsec=6., log_m_host=13.3,  r_tidal=0.25,
+        LOS_normalization=1.0, geometry_type='DOUBLE_CONE', kwargs_cosmo=None):
+    """
+
+    Generates realizations of SIDM halos in terms of the core collapse timescale in different mass bins using the
+    parametric model by Yang et al.
+
+    :param z_lens: main deflector redshift
+    :param z_source: source redshift
+    :param log10_mass_ranges: a list of lists specifying halo mass ranges, e.g. [[6, 8], [8, 10]]
+    :param collapse_timescales: a list of core collapse timescales in Gyr for halos in each bin, e.g. [6.0, 1.0]
+    :param subhalo_time_scaling: a number that makes time pass quicker (>1) or lower (<1) for subhalos relative to
+    field halos
+    :param sigma_sub: normalization of the subhalo mass function
+    :param log10_sigma_sub: normalization of the subhalo mass function in log units (overwrites sigma_sub)
+    :param log_mlow: log base 10 of the minimum halo mass to render
+    :param log_mhigh: log base 10 of the maximum halo mass to render
+    :param concentration_model_subhalos: the concentration-mass relation applied to subhalos,
+    see concentration_models.py for a complete list of available models
+    :param kwargs_concentration_model_subhalos: keyword arguments for the subhalo MC relation
+    NOTE: keyword args returned by the load_concentration_model override these keywords with duplicate arguments
+    :param concentration_model_subhalos: the concentration-mass relation applied to field halos,
+    see concentration_models.py for a complete list of available models
+    :param kwargs_concentration_model_fieldhalos: keyword arguments for the field halo MC relation
+    NOTE: keyword args returned by the load_concentration_model override these keywords with duplicate arguments
+    :param truncation_model_subhalos: the truncation model applied to subhalos, see truncation_models for a complete list
+    :param kwargs_truncation_model_subhalos: keyword arguments for the truncation model applied to subhalos
+    :param truncation_model_fieldhalos: the truncation model applied to field halos, see truncation_models for a
+    complete list
+    :param kwargs_truncation_model_fieldhalos: keyword arguments for the truncation model applied to field halos
+    :param infall_redshift_model: a string that specifies that infall redshift sampling distribution, see the LensCosmo
+    class for details
+    :param kwargs_infall_model: keyword arguments for the infall redshift model
+    :param shmf_log_slope: the logarithmic slope of the subhalo mass function pivoting around 10^8 M_sun
+    :param cone_opening_angle_arcsec: the opening angle of the rendering volume in arcsec
+    :param log_m_host: log base 10 of the host halo mass [M_sun]
+    :param r_tidal: the core radius of the host halo in units of the host halo scale radius. Subhalos are distributed
+    in 3D with a cored NFW profile with this core radius
+    :param LOS_normalization: rescales the amplitude of the line-of-sight mass function
+    :param geometry_type: string that specifies the geometry of the rendering volume; options include
+    DOUBLE_CONE, CONE, CYLINDER
+    :param kwargs_cosmo: keyword arguments that specify the cosmology, see Cosmology/cosmology.py
+    :param collapsed_halo_profile: string that sets the density profile of core-collapsed halos
+    currently implemented models are SPL_CORE and GNFW (see example notebook)
+    :param kwargs_collapsed_profile: keyword arguments for the collapsed profile (see example notebook)
+    :return: a realization of dark matter structure in SIDM
+    """
+    two_halo_contribution = True
+    delta_power_law_index = 0.0
+    cdm = CDM(z_lens, z_source, sigma_sub, log_mlow, log_mhigh, log10_sigma_sub,
+              concentration_model_subhalos, kwargs_concentration_model_subhalos,
+              concentration_model_fieldhalos, kwargs_concentration_model_fieldhalos,
+              truncation_model_subhalos, kwargs_truncation_model_subhalos,
+              truncation_model_fieldhalos, kwargs_truncation_model_fieldhalos,
+              infall_redshift_model, kwargs_infall_model,
+              shmf_log_slope, cone_opening_angle_arcsec, log_m_host, r_tidal,
+              LOS_normalization, two_halo_contribution, delta_power_law_index,
+              geometry_type, kwargs_cosmo)
+    extension = RealizationExtensions(cdm)
+    sidm = extension.toSIDM(log10_mass_ranges, collapse_timescales, subhalo_time_scaling)
+    return sidm

pyHalo/preset_models.py

@@ -22,6 +22,9 @@ def preset_model_from_name(name):
     elif name == 'SIDM_core_collapse':
         from pyHalo.PresetModels.sidm import SIDM_core_collapse
         return SIDM_core_collapse
+    elif name == 'SIDM_parametric':
+        from pyHalo.PresetModels.sidm import SIDM_parametric
+        return SIDM_parametric
     elif name == 'ULDM':
         from pyHalo.PresetModels.uldm import ULDM
         return ULDM