Theory prediction for miscentered haloes (#622)

Theoretical prediction if Sigma and DeltaSigma profiles allowing for a user-defined miscentering distance.

* Miscentering theory notebooks to explore implementation options (in examples/for_dev/ folder)
* Add miscentering as option of existing user functions, in functional and OO-interface; update docstrings
* Allow backend-dependent calculation
* Update functional interface  and OO theory demo with miscentering option
* New tests for miscentering
* Make for_dev example folder to store notebook meant for developement only

---------

Co-authored-by: Hsin Fan <[email protected]>
Co-authored-by: m-aguena <[email protected]>
Co-authored-by: Anna Niemiec <[email protected]>

clmm/__init__.py

@@ -26,4 +26,4 @@
 )
 from . import support
 
-__version__ = "1.14.5"
+__version__ = "1.14.6"

clmm/theory/func_layer.py

@@ -1,6 +1,7 @@
 """@file func_layer.py
 Main functions to encapsule oo calls
 """
+
 # pylint: disable=too-many-lines
 # pylint: disable=invalid-name
 # Thin functonal layer on top of the class implementation of CLMModeling .
@@ -108,6 +109,8 @@ def compute_surface_density(
     halo_profile_model="nfw",
     massdef="mean",
     alpha_ein=None,
+    r_mis=None,
+    mis_from_backend=False,
     verbose=False,
     use_projected_quad=False,
     validate_input=True,
@@ -119,6 +122,14 @@ def compute_surface_density(
 
     where :math:`\rho(r)` is the 3d density profile.
 
+    If the `r_mis` keyword is specified, this function computes the miscentered surface 
+    density instead as
+
+    .. math::
+        \Sigma_{\rm mis}(R, R_{\rm mis}) = \frac{1}{2\pi}\int_0^{2\pi} d\theta \,
+        \Sigma\left(\sqrt{R^2 + R_{\rm mis}^2 - 2 R R_{\rm mis} \cos\theta} \right)\;,
+
+
     Parameters
     ----------
     r_proj : array_like
@@ -152,6 +163,13 @@ def compute_surface_density(
         Option only available for the NumCosmo and CCL backends.
         If None, use the default value of the backend. (0.25 for the NumCosmo backend and a
         cosmology-dependent value for the CCL backend.)
+
+    r_mis : float, optional
+        Projected miscenter distance in :math:`M\!pc`
+    mis_from_backend : bool, optional
+        If True, use the projected surface density from the backend for miscentering
+        calculations. If False, use the (faster) CLMM exact analytical 
+        implementation instead. (Default: False)
     verbose : boolean, optional
         If True, the Einasto slope (alpha_ein) is printed out. Only available for the NC and CCL
         backends.
@@ -181,10 +199,12 @@ def compute_surface_density(
     _modeling_object.set_mass(mdelta)
     if halo_profile_model == "einasto" or alpha_ein is not None:
         _modeling_object.set_einasto_alpha(alpha_ein)
-    if halo_profile_model == "einasto" and _modeling_object.backend=="ccl":
+    if halo_profile_model == "einasto" and _modeling_object.backend == "ccl":
         _modeling_object.set_projected_quad(use_projected_quad)
 
-    sigma = _modeling_object.eval_surface_density(r_proj, z_cl, verbose=verbose)
+    sigma = _modeling_object.eval_surface_density(
+        r_proj, z_cl, r_mis=r_mis, mis_from_backend=mis_from_backend, verbose=verbose
+    )
 
     _modeling_object.validate_input = True
     return sigma
@@ -200,6 +220,8 @@ def compute_mean_surface_density(
     halo_profile_model="nfw",
     massdef="mean",
     alpha_ein=None,
+    r_mis=None,
+    mis_from_backend=False,
     verbose=False,
     validate_input=True,
 ):
@@ -239,6 +261,12 @@ def compute_mean_surface_density(
     alpha_ein : float, optional
         If `halo_profile_model=='einasto'`, set the value of the Einasto slope. Option only
         available for the NumCosmo backend
+    r_mis : float, optional
+        Projected miscenter distance in :math:`M\!pc`
+    mis_from_backend : bool, optional
+        If True, use the projected surface density from the backend for miscentering
+        calculations. If False, use the (faster) CLMM exact analytical 
+        implementation instead. (Default: False)
     verbose : boolean, optional
         If True, the Einasto slope (alpha_ein) is printed out. Only available for the NC and CCL
         backends.
@@ -266,7 +294,9 @@ def compute_mean_surface_density(
     if alpha_ein is not None:
         _modeling_object.set_einasto_alpha(alpha_ein)
 
-    sigma_bar = _modeling_object.eval_mean_surface_density(r_proj, z_cl, verbose=verbose)
+    sigma_bar = _modeling_object.eval_mean_surface_density(
+        r_proj, z_cl, r_mis=r_mis, mis_from_backend=mis_from_backend, verbose=verbose
+    )
 
     _modeling_object.validate_input = True
     return sigma_bar
@@ -282,6 +312,8 @@ def compute_excess_surface_density(
     halo_profile_model="nfw",
     massdef="mean",
     alpha_ein=None,
+    r_mis=None,
+    mis_from_backend=False,
     verbose=False,
     validate_input=True,
 ):
@@ -323,6 +355,12 @@ def compute_excess_surface_density(
         Option only available for the NumCosmo and CCL backends.
         If None, use the default value of the backend. (0.25 for the NumCosmo backend and a
         cosmology-dependent value for the CCL backend.)
+    r_mis : float, optional
+        Projected miscenter distance in :math:`M\!pc`
+    mis_from_backend : bool, optional
+        If True, use the projected surface density from the backend for miscentering
+        calculations. If False, use the (faster) CLMM exact analytical 
+        implementation instead. (Default: False)
     verbose : boolean, optional
         If True, the Einasto slope (alpha_ein) is printed out. Only available for the NC and CCL
         backends.
@@ -344,7 +382,9 @@ def compute_excess_surface_density(
     if halo_profile_model == "einasto" or alpha_ein is not None:
         _modeling_object.set_einasto_alpha(alpha_ein)
 
-    deltasigma = _modeling_object.eval_excess_surface_density(r_proj, z_cl, verbose=verbose)
+    deltasigma = _modeling_object.eval_excess_surface_density(
+        r_proj, z_cl, r_mis=r_mis, mis_from_backend=mis_from_backend, verbose=verbose
+    )
 
     _modeling_object.validate_input = True
     return deltasigma

clmm/theory/miscentering.py

@@ -0,0 +1,251 @@
+"""@file miscentering.py
+Model functions with miscentering
+"""
+# Functions to model halo profiles
+
+import numpy as np
+from scipy.integrate import quad, dblquad, tplquad
+
+
+def integrand_surface_density_nfw(theta, r_proj, r_mis, r_s):
+    r"""Computes integrand for surface mass density with the NFW profile.
+
+    Parameters
+    ----------
+    theta : float
+        Angle of polar coordinates of the miscentering direction.
+    r_proj : array_like
+        Projected radial position from the cluster center in :math:`M\!pc`.
+    r_mis : float, optional
+        Projected miscenter distance in :math:`M\!pc`.
+    r_s : array_like
+        Scale radius
+
+    Returns
+    -------
+    numpy.ndarray, float
+        2D projected density in units of :math:`M_\odot\ Mpc^{-2}`.
+    """
+    r_norm = np.sqrt(r_proj**2.0 + r_mis**2.0 - 2.0 * r_proj * r_mis * np.cos(theta)) / r_s
+
+    r2m1 = r_norm**2.0 - 1.0
+    if r_norm < 1:
+        sqrt_r2m1 = np.sqrt(-r2m1)
+        res = np.arcsinh(sqrt_r2m1 / r_norm) / (-r2m1) ** (3.0 / 2.0) + 1.0 / r2m1
+    elif r_norm > 1:
+        sqrt_r2m1 = np.sqrt(r2m1)
+        res = -np.arcsin(sqrt_r2m1 / r_norm) / (r2m1) ** (3.0 / 2.0) + 1.0 / r2m1
+    else:
+        res = 1.0 / 3.0
+    return res
+
+
+def integrand_surface_density_einasto(r_par, theta, r_proj, r_mis, r_s, alpha_ein):
+    r"""Computes integrand for surface mass density with the Einasto profile.
+
+    Parameters
+    ----------
+    r_par : array_like
+        Parallel radial position from the cluster center in :math:`M\!pc`.
+    theta : float
+        Angle of polar coordinates of the miscentering direction.
+    r_proj : array_like
+        Projected radial position from the cluster center in :math:`M\!pc`.
+    r_mis : float, optional
+        Projected miscenter distance in :math:`M\!pc`.
+    r_s : array_like
+        Scale radius
+    alpha_ein : float
+        Einasto slope
+
+    Returns
+    -------
+    numpy.ndarray, float
+        2D projected density in units of :math:`M_\odot\ Mpc^{-2}`.
+    """
+    # Projected surface mass density element for numerical integration
+    r_norm = (
+        np.sqrt(r_par**2.0 + r_proj**2.0 + r_mis**2.0 - 2.0 * r_proj * r_mis * np.cos(theta))
+        / r_s
+    )
+
+    return np.exp(-2.0 * (r_norm**alpha_ein - 1.0) / alpha_ein)
+
+
+def integrand_surface_density_hernquist(theta, r_proj, r_mis, r_s):
+    r"""Computes integrand for surface mass density with the Hernquist profile.
+
+    Parameters
+    ----------
+    theta : float
+        Angle of polar coordinates of the miscentering direction.
+    r_proj : array_like
+        Projected radial position from the cluster center in :math:`M\!pc`.
+    r_mis : float, optional
+        Projected miscenter distance in :math:`M\!pc`.
+    r_s : array_like
+        Scale radius
+
+    Returns
+    -------
+    numpy.ndarray, float
+        2D projected density in units of :math:`M_\odot\ Mpc^{-2}`.
+    """
+    r_norm = np.sqrt(r_proj**2.0 + r_mis**2.0 - 2.0 * r_proj * r_mis * np.cos(theta)) / r_s
+
+    r2m1 = r_norm**2.0 - 1.0
+    if r_norm < 1:
+        res = -3 / r2m1**2 + (r2m1 + 3) * np.arcsinh(np.sqrt(-r2m1) / r_norm) / (-r2m1) ** 2.5
+    elif r_norm > 1:
+        res = -3 / r2m1**2 + (r2m1 + 3) * np.arcsin(np.sqrt(r2m1) / r_norm) / (r2m1) ** 2.5
+    else:
+        res = 4.0 / 15.0
+    return res
+
+
+def integrand_mean_surface_density_nfw(theta, r_proj, r_mis, r_s):
+    r"""Computes integrand for mean surface mass density with the NFW profile.
+
+    Parameters
+    ----------
+    theta : float
+        Angle of polar coordinates of the miscentering direction.
+    r_proj : array_like
+        Projected radial position from the cluster center in :math:`M\!pc`.
+    r_mis : float, optional
+        Projected miscenter distance in :math:`M\!pc`.
+    r_s : array_like
+        Scale radius
+
+    Returns
+    -------
+    numpy.ndarray, float
+        Mean surface density in units of :math:`M_\odot\ Mpc^{-2}`.
+    """
+    return r_proj * integrand_surface_density_nfw(theta, r_proj, r_mis, r_s)
+
+
+def integrand_mean_surface_density_einasto(r_par, theta, r_proj, r_mis, r_s, alpha_ein):
+    r"""Computes integrand for mean surface mass density with the Einasto profile.
+
+    Parameters
+    ----------
+    r_par : array_like
+        Parallel radial position from the cluster center in :math:`M\!pc`.
+    theta : float
+        Angle of polar coordinates of the miscentering direction.
+    r_proj : array_like
+        Projected radial position from the cluster center in :math:`M\!pc`.
+    r_mis : float, optional
+        Projected miscenter distance in :math:`M\!pc`.
+    r_s : array_like
+        Scale radius
+    alpha_ein : float
+        Einasto slope
+
+    Returns
+    -------
+    numpy.ndarray, float
+        Mean surface density in units of :math:`M_\odot\ Mpc^{-2}`.
+    """
+    return r_proj * integrand_surface_density_einasto(r_par, theta, r_proj, r_mis, r_s, alpha_ein)
+
+
+def integrand_mean_surface_density_hernquist(theta, r_proj, r_mis, r_s):
+    r"""Computes integrand for mean surface mass density with the Hernquist profile.
+
+    Parameters
+    ----------
+    theta : float
+        Angle of polar coordinates of the miscentering direction.
+    r_proj : array_like
+        Projected radial position from the cluster center in :math:`M\!pc`.
+    r_mis : float, optional
+        Projected miscenter distance in :math:`M\!pc`.
+    r_s : array_like
+        Scale radius
+
+    Returns
+    -------
+    numpy.ndarray, float
+        Mean surface density in units of :math:`M_\odot\ Mpc^{-2}`.
+    """
+    return r_proj * integrand_surface_density_hernquist(theta, r_proj, r_mis, r_s)
+
+
+def integrate_azimuthially_miscentered_surface_density(
+    r_proj, r_mis, integrand, norm, aux_args, extra_integral
+):
+    r"""Integrates azimuthally the miscentered surface mass density kernel.
+
+    Parameters
+    ----------
+    r_proj : array_like
+        Projected radial position from the cluster center in :math:`M\!pc`.
+    r_mis : float, optional
+        Projected miscenter distance in :math:`M\!pc`.
+    integrand : function
+        Function to be integrated
+    norm : float
+        Normalization value for integral
+    aux_args : list
+        Auxiliary arguments used in the integral
+    extra_integral : bool
+        Additional dimention for the integral
+
+    Returns
+    -------
+    numpy.ndarray, float
+        2D projected density in units of :math:`M_\odot\ Mpc^{-2}`.
+    """
+    args_list = [(r, r_mis, *aux_args) for r in r_proj]
+    if extra_integral:
+        res = [
+            dblquad(integrand, 0.0, np.pi, 0, np.inf, args=args, epsrel=1e-6)[0]
+            for args in args_list
+        ]
+    else:
+        res = [quad(integrand, 0.0, np.pi, args=args, epsrel=1e-6)[0] for args in args_list]
+
+    res = np.array(res) * norm / np.pi
+    return res
+
+
+def integrate_azimuthially_miscentered_mean_surface_density(
+    r_proj, r_mis, integrand, norm, aux_args, extra_integral
+):
+    r"""Integrates azimuthally the miscentered mean surface mass density kernel.
+
+    Parameters
+    ----------
+    theta : float
+        Angle of polar coordinates of the miscentering direction.
+    r_proj : array_like
+        Projected radial position from the cluster center in :math:`M\!pc`.
+    r_mis : float, optional
+        Projected miscenter distance in :math:`M\!pc`.
+    r_s : array_like
+        Scale radius
+
+    Returns
+    -------
+    numpy.ndarray, float
+        Mean surface density in units of :math:`M_\odot\ Mpc^{-2}`.
+    """
+
+    r_lower = np.zeros_like(r_proj)
+    r_lower[1:] = r_proj[:-1]
+    args = (r_mis, *aux_args)
+
+    if extra_integral:
+        res = [
+            tplquad(integrand, r_low, r_high, 0, np.pi, 0, np.inf, args=args, epsrel=1e-6)[0]
+            for r_low, r_high in zip(r_lower, r_proj)
+        ]
+    else:
+        res = [
+            dblquad(integrand, r_low, r_high, 0, np.pi, args=args, epsrel=1e-6)[0]
+            for r_low, r_high in zip(r_lower, r_proj)
+        ]
+
+    return np.cumsum(res) * norm * 2 / np.pi / r_proj**2