Diemer 2023 (truncated exponential) profile

This module implements the Diemer 2022b form of the density profile. Please see Halo Density Profiles for a general introduction to the Colossus density profile module.

Basics

The Diemer (2022b) profile was designed to fit the orbiting component of dark matter halos, even at large radii where the infalling component comes to dominate. The orbiting component has a sharp truncation that cannot be described by a power-law steepening term (as in the DK14 profile), but it can be described by the following “truncated exponential” form:

\[\rho(r) = \rho_{\rm s} \exp \left\{ -\frac{2}{\alpha} \left[ \left( \frac{r}{r_{\rm s}} \right)^\alpha - 1 \right] -\frac{1}{\beta} \left[ \left( \frac{r}{r_{\rm t}} \right)^\beta - \left( \frac{r_{\rm s}}{r_{\rm t}} \right)^\beta \right] \right\}\]

The meaning of this functional form is easiest to understand by considering its logarithmic slope:

\[\gamma(r) \equiv \frac{{\rm d} \ln \rho}{{\rm d} \ln r} = -2 \left( \frac{r}{r_{\rm s}} \right)^\alpha - \left( \frac{r}{r_{\rm t}} \right)^\beta\]

The first term is identical to an Einasto profile, and the second term causes a more or less sharp truncation. The formula has 5 free parameters with well-defined physical interpretations:

Param.	Symbol	Explanation
rhos	\(\rho_s\)	Density at the scale radius, in physical \(M_{\odot} h^2 / {\rm kpc}^3\)
rs	\(r_{\rm s}\)	The scale radius in physical kpc/h
alpha	\(\alpha\)	Determines how quickly the slope of the inner Einasto profile steepens
rt	\(r_{\rm t}\)	The radius where the profile steepens beyond the Einasto profile, in physical kpc/h
beta	\(\beta\)	Sharpness of the truncation

As with all profile models, the user can pass these fundamental parameters or mass and concentration to the constructor of the ModelAProfile class (the reason for the name will become apparent later). In the latter case, the user can also give additional information to create a more accurate profile model. In particular, the fitting function was calibrated for the median and mean profiles of halo samples selected by mass (selected_by = 'M') and selected by both mass and mass accretion rate (selected_by = 'Gamma'). The latter option results in a more accurate representation of the density profile, but the mass accretion rate must be known. See the deriveParameters() function for details.

Adding an infalling profile

In most real-world applications, we are interested in the total density rather than only that of orbiting matter. We thus want to add the overdensity of matter on a first infall (or “infalling profile”), as well as the mean density of the Universe. Such a composite orbiting+infalling model can easily be created with the compositeProfile() function:

from colossus.cosmology import cosmology
from colossus.halo import profile_composite

cosmology.setCosmology('planck18')
p = profile_composite.compositeProfile('diemer23', outer_names = ['mean', 'infalling'],
                        M = 1E12, c = 10.0, z = 0.0, mdef = 'vir', pl_delta_1 = 10.0, pl_s = 1.5)

With a single command, we have created a truncated exponential profile with two outer terms, the constant mean density and an infalling profile of the form

\[\rho_{\rm inf}(r) = \delta_1 \rho_{\rm m}(z) \left[ \left( \frac{\delta_1}{\delta_{\rm max}} \right)^{1/\zeta} + \left( \frac{r}{r_{\rm pivot}} \right)^{s/\zeta} \right]^{-\zeta} \]

where \(\delta_1\) is the overdensity normalization and \(s\) the slope. These parameters depend on the mass, accretion rate, and cosmology of the halo sample in question (see Diemer 2022b). The maximum overdensity at the center can safely be left to its default value unless the infalling profile is known in detail, as can \(\zeta = 0.5\). The pivot radius is, by default, set to \(R_{\rm 200m}\). This parameterization relies on the parameters of the inner profile, which is correctly handled by the constructor. When fitting, however, such an interdependence can create issues and it is recommended to set a fixed physical radius as a pivot. For more details, see the OuterTermInfalling class, as well as the code Tutorials.

Model variant with correction at scale radius

The orbiting profile model described above has one technically unaesthetic property: the logarithmic slope at the scale radius is no longer -2. Thus, Diemer 2022b also proposed a corrected variant called Model B, which can be created using the ModelBProfile class. Here, an extra term has been inserted into the slope to ensure that it remains -2 at \(r_{\rm s}\),

\[\gamma(r) \equiv \frac{{\rm d} \ln \rho}{{\rm d} \ln r} = -2 \left( \frac{r}{r_{\rm s}} \right)^\alpha - \left( \frac{r}{r_{\rm t}} \right)^\beta + \left( \frac{r_{\rm s}}{r_{\rm t}} \right)^\beta \left( \frac{r}{r_{\rm s}} \right)^\eta\]

where \(\eta = 0.1\) is a nuissance parameter that determines how quickly the correction term vanishes at small radii. The density function also becomes somewhat more complicated, but the user can ignore the underlying equations and use Model B exactly as Model A. The differences are so small that the parameters have virtually the same meaning. The main advantage of Model B is that it can be more stable in fits to profiles with a poorly defined scale radius, that is, profiles with a slope that is roughly -2 across a wide range of radii. Otherwise, we recommend using Model A for most applications. Given the Model A/B split is implemented via an abstract class and two specific derived classes:

GenericD22Profile (should never be instantiated by the user)
ModelAProfile (the default)
ModelBProfile (if correction at scale radius is needed)

Module reference

class halo.profile_diemer23.GenericD22Profile(selected_by='M', Gamma=None, **kwargs)

Base class for truncated exponential profiles.

Generic profile class for methods that are common to the Model A and B variants. This class should never be instantiated by the user.

Methods

`MDelta`(z, mdef)	The spherical overdensity mass of a given mass definition.
`RDelta`(z, mdef)	The spherical overdensity radius of a given mass definition.
`RMDelta`(z, mdef)	The spherical overdensity radius and mass of a given mass definition.
`Rsteepest`([search_range])	The radius where the logarithmic slope of the density profile is steepest.
`Vmax`()	The maximum circular velocity, and the radius where it occurs.
`circularVelocity`(r)	The circular velocity, \(v_c \equiv \sqrt{GM(<r)/r}\).
`cumulativePdf`(r[, Rmax, z, mdef])	The cumulative distribution function of the profile.
`deltaSigma`(r[, interpolate, ...])	The excess surface density at radius r.
`deltaSigmaInner`(r[, interpolate, ...])	The excess surface density at radius r due to the inner profile.
`deltaSigmaOuter`(r[, interpolate, ...])	The excess surface density at radius r due to the outer profile.
`density`(r)	Density as a function of radius.
`densityDerivativeLin`(r)	The linear derivative of density, \(d \rho / dr\).
`densityDerivativeLinInner`(r)	The linear derivative of the inner density, \(d \rho_{\rm inner} / dr\).
`densityDerivativeLinOuter`(r)	The linear derivative of the outer density, \(d \rho_{\rm outer} / dr\).
`densityDerivativeLog`(r)	The logarithmic derivative of density, \(d \log(\rho) / d \log(r)\).
`densityDerivativeLogInner`(r)	The logarithmic derivative of the inner density, \(d \log(\rho_{\rm inner}) / d \log(r)\).
`densityDerivativeLogOuter`(r)	The logarithmic derivative of the outer density, \(d \log(\rho_{\rm outer}) / d \log(r)\).
`densityInner`(r)	Density of the inner profile as a function of radius.
`densityOuter`(r)	Density of the outer profile as a function of radius.
`deriveParameters`(selected_by[, nu200m, z, Gamma])	Calibration of the parameters \(\alpha\), \(\beta\), and \(r_{\rm t}\).
`enclosedMass`(r[, accuracy])	The mass enclosed within radius r.
`enclosedMassInner`(r[, accuracy])	The mass enclosed within radius r due to the inner profile term.
`enclosedMassOuter`(r[, accuracy])	The mass enclosed within radius r due to the outer profile term.
`fit`(r, q, quantity[, q_err, q_cov, method, ...])	Fit the density, mass, surface density, or DeltaSigma profile to a given set of data points.
`getParameterArray`([mask])	Returns an array of the profile parameters.
`setNativeParameters`(M, c, z, mdef[, ...])	Set the native parameters from mass and concentration (and optionally others).
`setParameterArray`(pars[, mask])	Set the profile parameters from an array.
`surfaceDensity`(r[, interpolate, accuracy, ...])	The projected surface density at radius r.
`surfaceDensityInner`(r[, interpolate, ...])	The projected surface density at radius r due to the inner profile.
`surfaceDensityOuter`(r[, interpolate, ...])	The projected surface density at radius r due to the outer profile.
`update`()	Update the profile object after a change in parameters or cosmology.
`updateR200m`()	Update the internally stored R200m after a parameter change.

static deriveParameters(selected_by, nu200m=None, z=None, Gamma=None)

Calibration of the parameters \(\alpha\), \(\beta\), and \(r_{\rm t}\).

This function determines the values of those parameters in the Diemer22 profile that can be calibrated based on mass, and potentially mass accretion rate. The latter is the stronger determinant of the profile shape, but may not always be available (e.g., for mass-selected samples).

We set \(\alpha = 0.18\) and \(\beta = 3\), which are the default parameters for individual halo profiles. However, they are not necessarily optimal for any type of averaged sample, where the optimal values vary. We do not calibrate \(\alpha\) with mass as suggested by Gao et al. 2008 because we do not reproduce this relation in our data (Diemer 2022c).

The truncation ratius \(r_{\rm t}\) is calibrated as suggested by DK14 for Gamma-selected samples. If selected_by = 'M', we use a new parametrization because the meaning of rt differs for the different slope parameters of the DK14 profile. If selected_by = 'Gamma', \(r_{\rm t}\) is calibrated from Gamma and z. The DK14 calibrations are based on slightly different definitions of peak height (\(\nu = \nu_{\rm vir}\)), accretion rate, and for a different fitting function. However, the resulting \(r_{\rm t}\) values are very similar to the forthcoming analysis in Diemer 2022c.

Parameters

selected_by: str: The halo sample to which this profile refers can be selected mass M or by accretion rate Gamma.
nu200m: float: The peak height of the halo for which the parameters are to be calibrated. This parameter only needs to be passed if selected_by == 'M'.
z: float: Redshift
Gamma: float: The mass accretion rate over the past dynamical time, which is defined as the crossing time (see the dynamicalTime() function or Diemer 2017 for details). The definition in the DK14 profile is slightly different, but the definitions are close enough that they can be used interchangeably without great loss of accuracy. The Gamma parameter only needs to be passed if selected_by == 'Gamma'.

Returns

alpha: float: The Einasto steepening parameter.
beta: float: The steepening of the truncation term.
rt_R200m: float: The truncation radius in units of R200m.

setNativeParameters(M, c, z, mdef, selected_by=None, Gamma=None, **kwargs)

Set the native parameters from mass and concentration (and optionally others).

The truncated exponential profile has five free parameters, which are set by this function. The mass and concentration must be given as \(M_{\rm 200m}\) and \(c_{\rm 200m}\). Other mass definitions demand iteration, which can be achieved with the initialization routine in the parent class. This function ignores the presence of outer profiles.

Parameters

M: float: Spherical overdensity mass in \(M_{\odot}/h\).
c: float: The concentration, \(c = R / r_{\rm s}\), corresponding to the given halo mass and mass definition.
z: float: Redshift
mdef: str: The mass definition in which M and c are given. See Halo Mass Definitions for details.
selected_by: str: The halo sample to which this profile refers can be selected mass M or by accretion rate Gamma.
Gamma: float: The mass accretion rate as defined in DK14. This parameter only needs to be passed if selected_by == 'Gamma'. See comments in deriveParameters() function above.

densityDerivativeLinInner(r)

The linear derivative of the inner density, \(d \rho_{\rm inner} / dr\).

For the truncated exponential profile, the logarithmic derivative is much easier to evaluate. Thus, this function converts the logarithmic to the linear derivative.

Parameters

r: array_like: Radius in physical kpc/h; can be a number or a numpy array.

Returns

derivative: array_like: The linear derivative in physical \(M_{\odot} h / {\rm kpc}^2\); has the same dimensions as r.

MDelta(z, mdef)

The spherical overdensity mass of a given mass definition.

Parameters

z: float: Redshift
mdef: str: The mass definition for which the spherical overdensity mass is computed. See Halo Mass Definitions for details.

Returns

M: float: Spherical overdensity mass in \(M_{\odot} /h\).

See also

RDelta: The spherical overdensity radius of a given mass definition.
RMDelta: The spherical overdensity radius and mass of a given mass definition.

RDelta(z, mdef)

The spherical overdensity radius of a given mass definition.

Parameters

z: float: Redshift
mdef: str: The mass definition for which the spherical overdensity radius is computed. See Halo Mass Definitions for details.

Returns

R: float: Spherical overdensity radius in physical kpc/h.

See also

MDelta: The spherical overdensity mass of a given mass definition.
RMDelta: The spherical overdensity radius and mass of a given mass definition.

RMDelta(z, mdef)

The spherical overdensity radius and mass of a given mass definition.

This is a wrapper for the RDelta() and MDelta() functions which returns both radius and mass.

Parameters

z: float: Redshift
mdef: str: The mass definition for which the spherical overdensity mass is computed. See Halo Mass Definitions for details.

Returns

R: float: Spherical overdensity radius in physical kpc/h.
M: float: Spherical overdensity mass in \(M_{\odot} /h\).

See also

RDelta: The spherical overdensity radius of a given mass definition.
MDelta: The spherical overdensity mass of a given mass definition.

Rsteepest(search_range=10.0)

The radius where the logarithmic slope of the density profile is steepest.

This function finds the radius where the logarithmic slope of the profile is minimal, within some very generous bounds. The function makes sense only if at least one outer term has been added because the inner profile steepens with radius without ever become shallower again (for any reasonable functional form).

The radius of steepest slope is often taken as a proxy for the splashback radius, \(R_{\rm sp}\), but this correspondence is only approximate because the \(R_{\rm steep}\) is the result of a tradeoff between the orbiting and infalling profiles, whereas the splashback radius is determined by the dynamics of orbiting particles. See the Splashback radius section for a detailed description of the splashback radius.

Parameters

search_range: float: When searching for the radius of steepest slope, search within this factor of \(R_{\rm 200m}\), \(r_{\rm s}\), or another initial guess, which is determined automatically.

Returns

Rsteep: float: The radius where the slope is steepest, in physical kpc/h.

Vmax()

The maximum circular velocity, and the radius where it occurs.

Returns

vmax: float: The maximum circular velocity in km / s.
rmax: float: The radius where fmax occurs, in physical kpc/h.

See also

circularVelocity: The circular velocity, \(v_c \equiv \sqrt{GM(<r)/r}\).

circularVelocity(r)

The circular velocity, \(v_c \equiv \sqrt{GM(<r)/r}\).

Parameters

r: array_like: Radius in physical kpc/h; can be a number or a numpy array.

Returns

vc: float: The circular velocity in km / s; has the same dimensions as r.

See also

Vmax: The maximum circular velocity, and the radius where it occurs.

cumulativePdf(r, Rmax=None, z=None, mdef=None)

The cumulative distribution function of the profile.

Some density profiles do not converge to a finite mass at large radius, and the distribution thus needs to be cut off. The user can specify either a radius (in physical kpc/h) where the profile is cut off, or a mass definition and redshift to compute this radius (e.g., the virial radius \(R_{vir}\) at z = 0).

Parameters

r: array_like: Radius in physical kpc/h; can be a number or a numpy array.
Rmax: float: The radius where to cut off the profile in physical kpc/h.
z: float: Redshift
mdef: str: The radius definition for the cut-off radius. See Halo Mass Definitions for details.

Returns

pdf: array_like: The probability for mass to lie within radius r; has the same dimensions as r.

deltaSigma(r, interpolate=True, interpolate_surface_density=True, accuracy=0.0001, min_r_interpolate=1e-06, max_r_interpolate=100000000.0, max_r_integrate=1e+20)

The excess surface density at radius r.

This quantity is useful in weak lensing studies, and is defined as \(\Delta\Sigma(R) = \Sigma(<R)-\Sigma(R)\) where \(\Sigma(<R)\) is the averaged surface density within R weighted by area,

\[\Delta\Sigma(R) = \frac{1}{\pi R^2} \int_0^{R} 2 \pi r \Sigma(r) dr - \Sigma(R)\]

Parameters

r: array_like: Radius in physical kpc/h; can be a number or a numpy array.
interpolate: bool: Use an interpolation table for the surface density during the integration. This can speed up the evaluation significantly, as the surface density can be expensive to evaluate.
interpolate_surface_density: bool: Use an interpolation table for density during the computation of the surface density. This should make the evaluation somewhat faster, but can fail for some density terms which are negative at particular radii.
accuracy: float: The minimum accuracy of the integration (used both to compute the surface density and average it to get DeltaSigma).
min_r_interpolate: float: The minimum radius in physical kpc/h from which the surface density profile is averaged.
max_r_interpolate: float: The maximum radius in physical kpc/h to which the density profile is integrated when using interpolating density.
max_r_integrate: float: The maximum radius in physical kpc/h to which the density profile is integrated when using exact densities.

Returns

DeltaSigma: array_like: The excess surface density at radius r, in physical \(M_{\odot} h/{\rm kpc}^2\); has the same dimensions as r.

deltaSigmaInner(r, interpolate=True, interpolate_surface_density=True, accuracy=0.0001, min_r_interpolate=1e-06, max_r_interpolate=100000000.0, max_r_integrate=1e+20)

The excess surface density at radius r due to the inner profile.

Parameters

r: array_like: Radius in physical kpc/h; can be a number or a numpy array.
interpolate: bool: Use an interpolation table for the surface density during the integration. This can speed up the evaluation significantly, as the surface density can be expensive to evaluate.
interpolate_surface_density: bool: Use an interpolation table for density during the computation of the surface density. This should make the evaluation somewhat faster, but can fail for some density terms which are negative at particular radii.
accuracy: float: The minimum accuracy of the integration (used both to compute the surface density and average it to get DeltaSigma).
min_r_interpolate: float: The minimum radius in physical kpc/h from which the surface density profile is averaged.
max_r_interpolate: float: The maximum radius in physical kpc/h to which the density profile is integrated when using interpolating density.
max_r_integrate: float: The maximum radius in physical kpc/h to which the density profile is integrated when using exact densities.

Returns

DeltaSigma: array_like: The excess surface density at radius r, in physical \(M_{\odot} h/{\rm kpc}^2\); has the same dimensions as r.

deltaSigmaOuter(r, interpolate=True, interpolate_surface_density=True, accuracy=0.0001, min_r_interpolate=1e-06, max_r_interpolate=100000000.0, max_r_integrate=1e+20)

The excess surface density at radius r due to the outer profile.

Parameters

r: array_like: Radius in physical kpc/h; can be a number or a numpy array.
interpolate: bool: Use an interpolation table for the surface density during the integration. This can speed up the evaluation significantly, as the surface density can be expensive to evaluate.
interpolate_surface_density: bool: Use an interpolation table for density during the computation of the surface density. This should make the evaluation somewhat faster, but can fail for some density terms which are negative at particular radii.
accuracy: float: The minimum accuracy of the integration (used both to compute the surface density and average it to get DeltaSigma).
min_r_interpolate: float: The minimum radius in physical kpc/h from which the surface density profile is averaged.
max_r_interpolate: float: The maximum radius in physical kpc/h to which the density profile is integrated when using interpolating density.
max_r_integrate: float: The maximum radius in physical kpc/h to which the density profile is integrated when using exact densities.

Returns

DeltaSigma: array_like: The excess surface density at radius r, in physical \(M_{\odot} h/{\rm kpc}^2\); has the same dimensions as r.

density(r)

Density as a function of radius.

Parameters

r: array_like: Radius in physical kpc/h; can be a number or a numpy array.

Returns

density: array_like: Density in physical \(M_{\odot} h^2 / {\rm kpc}^3\); has the same dimensions as r.

densityDerivativeLin(r)

The linear derivative of density, \(d \rho / dr\).

This function should generally not be overwritten by child classes since it handles the general case of adding up the contributions from the inner and outer terms.

Parameters

r: array_like: Radius in physical kpc/h; can be a number or a numpy array.

Returns

derivative: array_like: The linear derivative in physical \(M_{\odot} h / {\rm kpc}^2\); has the same dimensions as r.

densityDerivativeLinOuter(r)

The linear derivative of the outer density, \(d \rho_{\rm outer} / dr\).

This function should generally not be overwritten by child classes since it handles the general case of adding up the contributions from all outer profile terms.

Parameters

r: array_like: Radius in physical kpc/h; can be a number or a numpy array.

Returns

derivative: array_like: The linear derivative in physical \(M_{\odot} h / {\rm kpc}^2\); has the same dimensions as r.

densityDerivativeLog(r)

The logarithmic derivative of density, \(d \log(\rho) / d \log(r)\).

This function should generally not be overwritten by child classes since it handles the general case of adding up the contributions from the inner and outer profile terms.

Parameters

r: array_like: Radius in physical kpc/h; can be a number or a numpy array.

Returns

derivative: array_like: The dimensionless logarithmic derivative; has the same dimensions as r.

densityDerivativeLogInner(r)

The logarithmic derivative of the inner density, \(d \log(\rho_{\rm inner}) / d \log(r)\).

This function evaluates the logarithmic derivative based on the linear derivative. If there is an analytic expression for the logarithmic derivative, child classes should overwrite this function.

Parameters

r: array_like: Radius in physical kpc/h; can be a number or a numpy array.

Returns

derivative: array_like: The dimensionless logarithmic derivative; has the same dimensions as r.

densityDerivativeLogOuter(r)

The logarithmic derivative of the outer density, \(d \log(\rho_{\rm outer}) / d \log(r)\).

This function should generally not be overwritten by child classes since it handles the general case of adding up the contributions from outer profile terms.

Parameters

r: array_like: Radius in physical kpc/h; can be a number or a numpy array.

Returns

derivative: array_like: The dimensionless logarithmic derivative; has the same dimensions as r.

abstract densityInner(r)

Density of the inner profile as a function of radius.

Abstract function which must be overwritten by child classes.

Parameters

r: array_like: Radius in physical kpc/h; can be a number or a numpy array.

Returns

density: array_like: Density in physical \(M_{\odot} h^2 / {\rm kpc}^3\); has the same dimensions as r.

densityOuter(r)

Density of the outer profile as a function of radius.

This function should generally not be overwritten by child classes since it handles the general case of adding up the contributions from all outer profile terms.

Parameters

r: array_like: Radius in physical kpc/h; can be a number or a numpy array.

Returns

density: array_like: Density in physical \(M_{\odot} h^2 / {\rm kpc}^3\); has the same dimensions as r.

enclosedMass(r, accuracy=1e-06)

The mass enclosed within radius r.

Parameters

r: array_like: Radius in physical kpc/h; can be a number or a numpy array.
accuracy: float: The minimum accuracy of the integration.

Returns

M: array_like: The mass enclosed within radius r, in \(M_{\odot}/h\); has the same dimensions as r.

enclosedMassInner(r, accuracy=1e-06)

The mass enclosed within radius r due to the inner profile term.

This function should be overwritten by child classes if an analytical expression exists.

Parameters

r: array_like: Radius in physical kpc/h; can be a number or a numpy array.
accuracy: float: The minimum accuracy of the integration.

Returns

M: array_like: The mass enclosed within radius r, in \(M_{\odot}/h\); has the same dimensions as r.

enclosedMassOuter(r, accuracy=1e-06)

The mass enclosed within radius r due to the outer profile term.

Parameters

r: array_like: Radius in physical kpc/h; can be a number or a numpy array.
accuracy: float: The minimum accuracy of the integration.

Returns

M: array_like: The mass enclosed within radius r, in \(M_{\odot}/h\); has the same dimensions as r.

fit(r, q, quantity, q_err=None, q_cov=None, method='leastsq', mask=None, verbose=True, tolerance=1e-05, maxfev=None, leastsq_algorithm='trf', use_legacy_leastsq=False, bounds=None, initial_step=0.01, nwalkers=100, random_seed=None, convergence_step=100, converged_GR=0.01, best_fit='median', output_every_n=100)

Fit the density, mass, surface density, or DeltaSigma profile to a given set of data points.

This function represents a general interface for finding the best-fit parameters of a halo density profile given a set of data points. These points can represent a number of different physical quantities: quantity can either be density, enclosed mass, surface density, or Delta Sigma (rho, M, Sigma, or DeltaSigma).

The data points q at radii r can optionally have error bars, and the user can pass a full covariance matrix. Please note that not passing any estimate of the uncertainty, i.e. q_err = None and q_cov = None, can lead to very poor fit results: the fitter will minimize the absolute difference between points, strongly favoring the high densities at the center.

The user can choose to vary only a sub-set of the profile parameters through the mask parameter. The current parameters of the profile instance serve as an initial guess.

By default, the parameters are transformed into log space during fitting to ensure positivity. Child classes can change this behavior by overwriting the _fitConvertParams() and _fitConvertParamsBack() functions.

There are two fundamental methods for performing the fit, a least-squares minimization (method = 'leastsq') and a Markov-Chain Monte Carlo (method = 'mcmc'). Both variants obey a few specific options (see below). The function returns a dictionary with outputs that somewhat depend on which method is chosen. After the function has completed, the profile instance represents the best-fit profile to the data points (i.e., its parameters are the best-fit parameters). Note that all output parameters are bundled into one dictionary. The explanations below refer to the entries in this dictionary.

Parameters

r: array_like: The radii of the data points, in physical kpc/h.
q: array_like: The data to fit; can either be density in physical \(M_{\odot} h^2 / {\rm kpc}^3\), enclosed mass in \(M_{\odot} /h\), or surface density in physical \(M_{\odot} h/{\rm kpc}^2\). Must have the same dimensions as r.
quantity: str: Indicates which quantity is given in as input in q, can be rho, M, Sigma, or DeltaSigma.
q_err: array_like: Optional; the uncertainty on the values in q in the same units. If method == 'mcmc', either q_err or q_cov must be passed. If method == 'leastsq' and neither q_err nor q_cov are passed, the absolute different between data points and fit is minimized. In this case, the returned chi2 is in units of absolute difference, meaning its value will depend on the units of q.
q_cov: array_like: Optional; the covariance matrix of the elements in q, as a 2-dimensional numpy array. This array must have dimensions of q**2 and be in units of the square of the units of q. If q_cov is passed, q_err is ignored since the diagonal elements of q_cov correspond to q_err**2.
method: str: The fitting method; can be leastsq for a least-squares minimization of mcmc for a Markov-Chain Monte Carlo.
mask: array_like: Optional; a numpy array of booleans that has the same length as the variables vector of the density profile class. Only variables where mask == True are varied in the fit, all others are kept constant. Important: this has to be a numpy array rather than a list.
verbose: bool / int: If true, output information about the fitting process. The flag can also be set as a number, where 1 has the same effect as True, and 2 outputs large amounts of information such as the fit parameters at each iteration.
tolerance: float: Only active when method == 'leastsq'. The accuracy to which the best-fit parameters are found.
maxfev: int: Only active when method == 'leastsq'. The maximum number of function evaluations before the fit is aborted. If None, the default value of the scipy least_squares function is used.
leastsq_algorithm: str: Only active when method == 'leastsq'. Can be any of the methods accepted by the scipy least_squares() function. Default is trf (trust region reflective), which works well for profile fits.
use_legacy_leastsq: bool: Only active when method == 'leastsq'. If True, this setting falls back to the old leastsq() in scipy rather than using the newer least_squares(). Should only be used for backward compatibility.
bounds: array_like: Only active when method == 'leastsq' and use_legacy_leastsq == False. If not None, this parameter must be an array of dimensions [2, N_par], giving two sets (lower and upper limits) of the fitted parameters (not the entire parameter array, if a mask is imposed). The limits must be given in linear space. If the parameters are fitted in log space (or some other transformation), that transformation is automatically applied to the bounds. For example, when fitting in log space (the default), the lower bounds must be positive, but the upper bounds can be np.inf.
initial_step: array_like: Only active when method == 'mcmc'. The MCMC samples (“walkers”) are initially distributed in a Gaussian around the initial guess. The width of the Gaussian is given by initial_step, either as an array of length N_par (giving the width of each Gaussian) or as a float number, in which case the width is set to initial_step times the initial value of the parameter.
nwalkers: int: Only active when method == 'mcmc'. The number of MCMC samplers that are run in parallel.
random_seed: int: Only active when method == 'mcmc'. If random_seed is not None, it is used to initialize the random number generator. This can be useful for reproducing results.
convergence_step: int: Only active when method == 'mcmc'. The convergence criteria are computed every convergence_step steps (and output is printed if verbose == True).
converged_GR: float: Only active when method == 'mcmc'. The maximum difference between different chains, according to the Gelman-Rubin criterion. Once the GR indicator is lower than this number in all parameters, the chain is ended. Setting this number too low leads to very long runtimes, but setting it too high can lead to inaccurate results.
best_fit: str: Only active when method == 'mcmc'. This parameter determines whether the mean or median value of the likelihood distribution is used as the output parameter set.
output_every_n: int: Only active when method == 'mcmc'. This parameter determines how frequently the MCMC chain outputs information. Only effective if verbose == True.

Returns

results: dict

A dictionary bundling the various fit results. Regardless of the fitting method, the dictionary always contains the following entries:

x: array_like: The best-fit result vector. If mask is passed, this vector only contains those variables that were varied in the fit.
q_fit: array_like: The fitted profile at the radii of the data points; has the same units as q and the same dimensions as r.
chi2: float: The chi^2 of the best-fit profile. If a covariance matrix was passed, the covariances are taken into account. If no uncertainty was passed at all, chi2 is in units of absolute difference, meaning its value will depend on the units of q.
ndof: int: The number of degrees of freedom, i.e. the number of fitted data points minus the number of free parameters.
chi2_ndof: float: The chi^2 per degree of freedom.

If method == 'leastsq', the dictionary additionally contains the entries returned by scipy.optimize.leastsq as well as the following:

nfev: int: The number of function calls used in the fit.
x_err: array_like: An array of dimensions [2, nparams] which contains an estimate of the lower and upper uncertainties on the fitted parameters. These uncertainties are computed from the covariance matrix estimated by the fitter. Please note that this estimate does not exactly correspond to a 68% likelihood. In order to get more statistically meaningful uncertainties, please use the MCMC samples instead of least-squares. In some cases, the fitter fails to return a covariance matrix, in which case x_err is None.

If method == 'mcmc', the dictionary contains the following entries:

x_initial: array_like: The initial positions of the walkers, in an array of dimensions [nwalkers, nparams].
chain_full: array_like: A numpy array of dimensions [n_independent_samples, nparams] with the parameters at each step in the chain. In this thin chain, only every nth step is output, where n is the auto-correlation time, meaning that the samples in this chain are truly independent.
chain_thin: array_like: Like the thin chain, but including all steps. Thus, the samples in this chain are not indepedent from each other. However, the full chain often gives better plotting results.
R: array_like: A numpy array containing the GR indicator at each step when it was saved.
x_mean: array_like: The mean of the chain for each parameter; has length nparams.
x_median: array_like: The median of the chain for each parameter; has length nparams.
x_stddev: array_like: The standard deviation of the chain for each parameter; has length nparams.
x_percentiles: array_like: The lower and upper values of each parameter that contain a certain percentile of the probability; has dimensions [n_percentages, 2, nparams] where the second dimension contains the lower/upper values.

getParameterArray(mask=None)

Returns an array of the profile parameters.

The profile parameters are internally stored in an ordered dictionary. For some applications (e.g., fitting), a simple array is more appropriate.

Parameters

mask: array_like: Optional; must be a numpy array (not a list) of booleans, with the same length as the parameter vector of the profile class (profile.N_par). Only those parameters that correspond to True values are returned.

Returns

par: array_like: A numpy array with the profile’s parameter values.

setParameterArray(pars, mask=None)

Set the profile parameters from an array.

The profile parameters are internally stored in an ordered dictionary. For some applications (e.g., fitting), setting them directly from an array might be necessary. If the profile contains values that depend on the parameters, the profile class must overwrite this function and update according to the new parameters.

Parameters

pars: array_like: The new parameter array.
mask: array_like: Optional; must be a numpy array (not a list) of booleans, with the same length as the parameter vector of the profile class (profile.N_par). If passed, only those parameters that correspond to True values are set (meaning the pars parameter must be shorter than profile.N_par).

surfaceDensity(r, interpolate=True, accuracy=0.0001, max_r_interpolate=100000000.0, max_r_integrate=1e+20)

The projected surface density at radius r.

The surface density is computed by projecting the 3D density along the line of sight,

\[\Sigma(R) = 2 \int_R^{\infty} \frac{r \rho(r)}{\sqrt{r^2-R^2}} dr\]

Parameters

r: array_like: Radius in physical kpc/h; can be a number or a numpy array.
interpolate: bool: Use an interpolation table for density during the integration. This should make the evaluation somewhat faster, depending on how large the radius array is.
accuracy: float: The minimum accuracy of the integration.
max_r_interpolate: float: The maximum radius in physical kpc/h to which the density profile is integrated when using interpolating density.
max_r_integrate: float: The maximum radius in physical kpc/h to which the density profile is integrated when using exact densities.

Returns

Sigma: array_like: The surface density at radius r, in physical \(M_{\odot} h/{\rm kpc}^2\); has the same dimensions as r.

surfaceDensityInner(r, interpolate=True, accuracy=0.0001, max_r_interpolate=100000000.0, max_r_integrate=1e+20)

The projected surface density at radius r due to the inner profile.

Parameters

r: array_like: Radius in physical kpc/h; can be a number or a numpy array.
interpolate: bool: Use an interpolation table for density during the integration. This should make the evaluation somewhat faster, depending on how large the radius array is.
accuracy: float: The minimum accuracy of the integration.
max_r_interpolate: float: The maximum radius in physical kpc/h to which the density profile is integrated when using interpolating density.
max_r_integrate: float: The maximum radius in physical kpc/h to which the density profile is integrated when using exact densities.

Returns

Sigma: array_like: The surface density at radius r, in physical \(M_{\odot} h/{\rm kpc}^2\); has the same dimensions as r.

surfaceDensityOuter(r, interpolate=True, accuracy=0.0001, max_r_interpolate=100000000.0, max_r_integrate=1e+20)

The projected surface density at radius r due to the outer profile.

This function checks whether there are explicit expressions for the surface density of the outer profile terms available, and uses them if possible. Note that there are some outer terms whose surface density integrates to infinity, such as the mean density of the universe which is constant to infinitely large radii.

Parameters

r: array_like: Radius in physical kpc/h; can be a number or a numpy array.
interpolate: bool: Use an interpolation table for density during the integration. This should make the evaluation somewhat faster, depending on how large the radius array is.
accuracy: float: The minimum accuracy of the integration.
max_r_interpolate: float: The maximum radius in physical kpc/h to which the density profile is integrated when using interpolating density.
max_r_integrate: float: The maximum radius in physical kpc/h to which the density profile is integrated when using exact densities.

Returns

Sigma: array_like: The surface density at radius r, in physical \(M_{\odot} h/{\rm kpc}^2\); has the same dimensions as r.

update()

Update the profile object after a change in parameters or cosmology.

If the parameters dictionary has been changed (e.g. by the user or during fitting), this function must be called to ensure consistency within the profile object. This involves deleting any pre-computed quantities (e.g., tabulated enclosed masses) and re-computing profile properties that depend on the parameters. The corresponding functions for outer terms are automatically called as well.

updateR200m()

Update the internally stored R200m after a parameter change.

If the profile has the internal option opt['R200m'] option, that does not stay in sync with the other profile parameters if they are changed (either inside or outside the constructor). This function adjusts \(R_{\rm 200m}\), in addition to whatever action is taken in the update function of the super class. Note that this adjustment needs to be done iteratively if any outer profiles rely on \(R_{\rm 200m}\).

class halo.profile_diemer23.ModelAProfile(selected_by='M', Gamma=None, **kwargs)

The Diemer 2023 (truncated exponential) density profile (default version).

The redshift must always be passed to this constructor, regardless of whether the fundamental parameters or a mass and concentration are given.

Parameters

rhos: float: The central scale density, in physical \(M_{\odot} h^2 / {\rm kpc}^3\).
rs: float: The scale radius in physical kpc/h.
rt: float: The radius where the profile steepens, in physical kpc/h.
alpha: float: Determines how quickly the slope of the inner Einasto profile steepens.
beta: float: Sharpness of the steepening.
M: float: Halo mass in \(M_{\odot}/h\).
c: float: Concentration in the same mass definition as M.
mdef: str: The mass definition to which M corresponds. See Halo Mass Definitions for details.
z: float: Redshift
selected_by: str: The halo sample to which this profile refers can be selected mass M or by accretion rate Gamma. This parameter influences how some of the fixed parameters in the profile are set, in particular those that describe the steepening term.
Gamma: float: The mass accretion rate over the past dynamical time, which is defined as the crossing time (see the dynamicalTime() function or Diemer 2017 for details). The definition in the DK14 profile is slightly different, but the definitions are close enough that they can be used interchangeably without great loss of accuracy. The Gamma parameter only needs to be passed if selected_by == 'Gamma'.

Methods

`MDelta`(z, mdef)	The spherical overdensity mass of a given mass definition.
`RDelta`(z, mdef)	The spherical overdensity radius of a given mass definition.
`RMDelta`(z, mdef)	The spherical overdensity radius and mass of a given mass definition.
`Rsteepest`([search_range])	The radius where the logarithmic slope of the density profile is steepest.
`Vmax`()	The maximum circular velocity, and the radius where it occurs.
`circularVelocity`(r)	The circular velocity, \(v_c \equiv \sqrt{GM(<r)/r}\).
`cumulativePdf`(r[, Rmax, z, mdef])	The cumulative distribution function of the profile.
`deltaSigma`(r[, interpolate, ...])	The excess surface density at radius r.
`deltaSigmaInner`(r[, interpolate, ...])	The excess surface density at radius r due to the inner profile.
`deltaSigmaOuter`(r[, interpolate, ...])	The excess surface density at radius r due to the outer profile.
`density`(r)	Density as a function of radius.
`densityDerivativeLin`(r)	The linear derivative of density, \(d \rho / dr\).
`densityDerivativeLinInner`(r)	The linear derivative of the inner density, \(d \rho_{\rm inner} / dr\).
`densityDerivativeLinOuter`(r)	The linear derivative of the outer density, \(d \rho_{\rm outer} / dr\).
`densityDerivativeLog`(r)	The logarithmic derivative of density, \(d \log(\rho) / d \log(r)\).
`densityDerivativeLogInner`(r)	The logarithmic derivative of the inner density, \(d \log(\rho_{\rm inner}) / d \log(r)\).
`densityDerivativeLogOuter`(r)	The logarithmic derivative of the outer density, \(d \log(\rho_{\rm outer}) / d \log(r)\).
`densityInner`(r)	Density of the inner profile as a function of radius.
`densityOuter`(r)	Density of the outer profile as a function of radius.
`deriveParameters`(selected_by[, nu200m, z, Gamma])	Calibration of the parameters \(\alpha\), \(\beta\), and \(r_{\rm t}\).
`enclosedMass`(r[, accuracy])	The mass enclosed within radius r.
`enclosedMassInner`(r[, accuracy])	The mass enclosed within radius r due to the inner profile term.
`enclosedMassOuter`(r[, accuracy])	The mass enclosed within radius r due to the outer profile term.
`fit`(r, q, quantity[, q_err, q_cov, method, ...])	Fit the density, mass, surface density, or DeltaSigma profile to a given set of data points.
`getParameterArray`([mask])	Returns an array of the profile parameters.
`setNativeParameters`(M, c, z, mdef[, ...])	Set the native parameters from mass and concentration (and optionally others).
`setParameterArray`(pars[, mask])	Set the profile parameters from an array.
`surfaceDensity`(r[, interpolate, accuracy, ...])	The projected surface density at radius r.
`surfaceDensityInner`(r[, interpolate, ...])	The projected surface density at radius r due to the inner profile.
`surfaceDensityOuter`(r[, interpolate, ...])	The projected surface density at radius r due to the outer profile.
`update`()	Update the profile object after a change in parameters or cosmology.
`updateR200m`()	Update the internally stored R200m after a parameter change.

densityInner(r)

Density of the inner profile as a function of radius.

Parameters

r: array_like: Radius in physical kpc/h; can be a number or a numpy array.

Returns

density: array_like: Density in physical \(M_{\odot} h^2 / {\rm kpc}^3\); has the same dimensions as r.

densityDerivativeLogInner(r)

The logarithmic derivative of the inner density, \(d \log(\rho_{\rm inner}) / d \log(r)\).

Parameters

r: array_like: Radius in physical kpc/h; can be a number or a numpy array.

Returns

derivative: array_like: The dimensionless logarithmic derivative; has the same dimensions as r.

MDelta(z, mdef)

The spherical overdensity mass of a given mass definition.

Parameters

z: float: Redshift
mdef: str: The mass definition for which the spherical overdensity mass is computed. See Halo Mass Definitions for details.

Returns

M: float: Spherical overdensity mass in \(M_{\odot} /h\).

See also

RDelta: The spherical overdensity radius of a given mass definition.
RMDelta: The spherical overdensity radius and mass of a given mass definition.

RDelta(z, mdef)

The spherical overdensity radius of a given mass definition.

Parameters

z: float: Redshift
mdef: str: The mass definition for which the spherical overdensity radius is computed. See Halo Mass Definitions for details.

Returns

R: float: Spherical overdensity radius in physical kpc/h.

See also

MDelta: The spherical overdensity mass of a given mass definition.
RMDelta: The spherical overdensity radius and mass of a given mass definition.

RMDelta(z, mdef)

The spherical overdensity radius and mass of a given mass definition.

This is a wrapper for the RDelta() and MDelta() functions which returns both radius and mass.

Parameters

z: float: Redshift
mdef: str: The mass definition for which the spherical overdensity mass is computed. See Halo Mass Definitions for details.

Returns

R: float: Spherical overdensity radius in physical kpc/h.
M: float: Spherical overdensity mass in \(M_{\odot} /h\).

See also

RDelta: The spherical overdensity radius of a given mass definition.
MDelta: The spherical overdensity mass of a given mass definition.

Rsteepest(search_range=10.0)

The radius where the logarithmic slope of the density profile is steepest.

Parameters

search_range: float: When searching for the radius of steepest slope, search within this factor of \(R_{\rm 200m}\), \(r_{\rm s}\), or another initial guess, which is determined automatically.

Returns

Rsteep: float: The radius where the slope is steepest, in physical kpc/h.

Vmax()

The maximum circular velocity, and the radius where it occurs.

Returns

vmax: float: The maximum circular velocity in km / s.
rmax: float: The radius where fmax occurs, in physical kpc/h.

See also

circularVelocity: The circular velocity, \(v_c \equiv \sqrt{GM(<r)/r}\).

circularVelocity(r)

The circular velocity, \(v_c \equiv \sqrt{GM(<r)/r}\).

Parameters

r: array_like: Radius in physical kpc/h; can be a number or a numpy array.

Returns

vc: float: The circular velocity in km / s; has the same dimensions as r.

See also

Vmax: The maximum circular velocity, and the radius where it occurs.

cumulativePdf(r, Rmax=None, z=None, mdef=None)

The cumulative distribution function of the profile.

Parameters

r: array_like: Radius in physical kpc/h; can be a number or a numpy array.
Rmax: float: The radius where to cut off the profile in physical kpc/h.
z: float: Redshift
mdef: str: The radius definition for the cut-off radius. See Halo Mass Definitions for details.

Returns

pdf: array_like: The probability for mass to lie within radius r; has the same dimensions as r.

deltaSigma(r, interpolate=True, interpolate_surface_density=True, accuracy=0.0001, min_r_interpolate=1e-06, max_r_interpolate=100000000.0, max_r_integrate=1e+20)

The excess surface density at radius r.

This quantity is useful in weak lensing studies, and is defined as \(\Delta\Sigma(R) = \Sigma(<R)-\Sigma(R)\) where \(\Sigma(<R)\) is the averaged surface density within R weighted by area,

\[\Delta\Sigma(R) = \frac{1}{\pi R^2} \int_0^{R} 2 \pi r \Sigma(r) dr - \Sigma(R)\]

Parameters

r: array_like: Radius in physical kpc/h; can be a number or a numpy array.
interpolate: bool: Use an interpolation table for the surface density during the integration. This can speed up the evaluation significantly, as the surface density can be expensive to evaluate.
interpolate_surface_density: bool: Use an interpolation table for density during the computation of the surface density. This should make the evaluation somewhat faster, but can fail for some density terms which are negative at particular radii.
accuracy: float: The minimum accuracy of the integration (used both to compute the surface density and average it to get DeltaSigma).
min_r_interpolate: float: The minimum radius in physical kpc/h from which the surface density profile is averaged.
max_r_interpolate: float: The maximum radius in physical kpc/h to which the density profile is integrated when using interpolating density.
max_r_integrate: float: The maximum radius in physical kpc/h to which the density profile is integrated when using exact densities.

Returns

DeltaSigma: array_like: The excess surface density at radius r, in physical \(M_{\odot} h/{\rm kpc}^2\); has the same dimensions as r.

deltaSigmaInner(r, interpolate=True, interpolate_surface_density=True, accuracy=0.0001, min_r_interpolate=1e-06, max_r_interpolate=100000000.0, max_r_integrate=1e+20)

The excess surface density at radius r due to the inner profile.

Parameters

r: array_like: Radius in physical kpc/h; can be a number or a numpy array.
interpolate: bool: Use an interpolation table for the surface density during the integration. This can speed up the evaluation significantly, as the surface density can be expensive to evaluate.
interpolate_surface_density: bool: Use an interpolation table for density during the computation of the surface density. This should make the evaluation somewhat faster, but can fail for some density terms which are negative at particular radii.
accuracy: float: The minimum accuracy of the integration (used both to compute the surface density and average it to get DeltaSigma).
min_r_interpolate: float: The minimum radius in physical kpc/h from which the surface density profile is averaged.
max_r_interpolate: float: The maximum radius in physical kpc/h to which the density profile is integrated when using interpolating density.
max_r_integrate: float: The maximum radius in physical kpc/h to which the density profile is integrated when using exact densities.

Returns

DeltaSigma: array_like: The excess surface density at radius r, in physical \(M_{\odot} h/{\rm kpc}^2\); has the same dimensions as r.

deltaSigmaOuter(r, interpolate=True, interpolate_surface_density=True, accuracy=0.0001, min_r_interpolate=1e-06, max_r_interpolate=100000000.0, max_r_integrate=1e+20)

The excess surface density at radius r due to the outer profile.

Parameters

r: array_like: Radius in physical kpc/h; can be a number or a numpy array.
interpolate: bool: Use an interpolation table for the surface density during the integration. This can speed up the evaluation significantly, as the surface density can be expensive to evaluate.
interpolate_surface_density: bool: Use an interpolation table for density during the computation of the surface density. This should make the evaluation somewhat faster, but can fail for some density terms which are negative at particular radii.
accuracy: float: The minimum accuracy of the integration (used both to compute the surface density and average it to get DeltaSigma).
min_r_interpolate: float: The minimum radius in physical kpc/h from which the surface density profile is averaged.
max_r_interpolate: float: The maximum radius in physical kpc/h to which the density profile is integrated when using interpolating density.
max_r_integrate: float: The maximum radius in physical kpc/h to which the density profile is integrated when using exact densities.

Returns

DeltaSigma: array_like: The excess surface density at radius r, in physical \(M_{\odot} h/{\rm kpc}^2\); has the same dimensions as r.

density(r)

Density as a function of radius.

Parameters

r: array_like: Radius in physical kpc/h; can be a number or a numpy array.

Returns

density: array_like: Density in physical \(M_{\odot} h^2 / {\rm kpc}^3\); has the same dimensions as r.

densityDerivativeLin(r)

The linear derivative of density, \(d \rho / dr\).

This function should generally not be overwritten by child classes since it handles the general case of adding up the contributions from the inner and outer terms.

Parameters

r: array_like: Radius in physical kpc/h; can be a number or a numpy array.

Returns

derivative: array_like: The linear derivative in physical \(M_{\odot} h / {\rm kpc}^2\); has the same dimensions as r.

densityDerivativeLinInner(r)

The linear derivative of the inner density, \(d \rho_{\rm inner} / dr\).

For the truncated exponential profile, the logarithmic derivative is much easier to evaluate. Thus, this function converts the logarithmic to the linear derivative.

Parameters

r: array_like: Radius in physical kpc/h; can be a number or a numpy array.

Returns

derivative: array_like: The linear derivative in physical \(M_{\odot} h / {\rm kpc}^2\); has the same dimensions as r.

densityDerivativeLinOuter(r)

The linear derivative of the outer density, \(d \rho_{\rm outer} / dr\).

This function should generally not be overwritten by child classes since it handles the general case of adding up the contributions from all outer profile terms.

Parameters

r: array_like: Radius in physical kpc/h; can be a number or a numpy array.

Returns

derivative: array_like: The linear derivative in physical \(M_{\odot} h / {\rm kpc}^2\); has the same dimensions as r.

densityDerivativeLog(r)

The logarithmic derivative of density, \(d \log(\rho) / d \log(r)\).

This function should generally not be overwritten by child classes since it handles the general case of adding up the contributions from the inner and outer profile terms.

Parameters

r: array_like: Radius in physical kpc/h; can be a number or a numpy array.

Returns

derivative: array_like: The dimensionless logarithmic derivative; has the same dimensions as r.

densityDerivativeLogOuter(r)

The logarithmic derivative of the outer density, \(d \log(\rho_{\rm outer}) / d \log(r)\).

This function should generally not be overwritten by child classes since it handles the general case of adding up the contributions from outer profile terms.

Parameters

r: array_like: Radius in physical kpc/h; can be a number or a numpy array.

Returns

derivative: array_like: The dimensionless logarithmic derivative; has the same dimensions as r.

densityOuter(r)

Density of the outer profile as a function of radius.

This function should generally not be overwritten by child classes since it handles the general case of adding up the contributions from all outer profile terms.

Parameters

r: array_like: Radius in physical kpc/h; can be a number or a numpy array.

Returns

density: array_like: Density in physical \(M_{\odot} h^2 / {\rm kpc}^3\); has the same dimensions as r.

static deriveParameters(selected_by, nu200m=None, z=None, Gamma=None)

Calibration of the parameters \(\alpha\), \(\beta\), and \(r_{\rm t}\).

Parameters

selected_by: str: The halo sample to which this profile refers can be selected mass M or by accretion rate Gamma.
nu200m: float: The peak height of the halo for which the parameters are to be calibrated. This parameter only needs to be passed if selected_by == 'M'.
z: float: Redshift
Gamma: float: The mass accretion rate over the past dynamical time, which is defined as the crossing time (see the dynamicalTime() function or Diemer 2017 for details). The definition in the DK14 profile is slightly different, but the definitions are close enough that they can be used interchangeably without great loss of accuracy. The Gamma parameter only needs to be passed if selected_by == 'Gamma'.

Returns

alpha: float: The Einasto steepening parameter.
beta: float: The steepening of the truncation term.
rt_R200m: float: The truncation radius in units of R200m.

enclosedMass(r, accuracy=1e-06)

The mass enclosed within radius r.

Parameters

r: array_like: Radius in physical kpc/h; can be a number or a numpy array.
accuracy: float: The minimum accuracy of the integration.

Returns

M: array_like: The mass enclosed within radius r, in \(M_{\odot}/h\); has the same dimensions as r.

enclosedMassInner(r, accuracy=1e-06)

The mass enclosed within radius r due to the inner profile term.

This function should be overwritten by child classes if an analytical expression exists.

Parameters

r: array_like: Radius in physical kpc/h; can be a number or a numpy array.
accuracy: float: The minimum accuracy of the integration.

Returns

M: array_like: The mass enclosed within radius r, in \(M_{\odot}/h\); has the same dimensions as r.

enclosedMassOuter(r, accuracy=1e-06)

The mass enclosed within radius r due to the outer profile term.

Parameters

r: array_like: Radius in physical kpc/h; can be a number or a numpy array.
accuracy: float: The minimum accuracy of the integration.

Returns

M: array_like: The mass enclosed within radius r, in \(M_{\odot}/h\); has the same dimensions as r.

Fit the density, mass, surface density, or DeltaSigma profile to a given set of data points.

The user can choose to vary only a sub-set of the profile parameters through the mask parameter. The current parameters of the profile instance serve as an initial guess.

Parameters

r: array_like: The radii of the data points, in physical kpc/h.
q: array_like: The data to fit; can either be density in physical \(M_{\odot} h^2 / {\rm kpc}^3\), enclosed mass in \(M_{\odot} /h\), or surface density in physical \(M_{\odot} h/{\rm kpc}^2\). Must have the same dimensions as r.
quantity: str: Indicates which quantity is given in as input in q, can be rho, M, Sigma, or DeltaSigma.
q_err: array_like: Optional; the uncertainty on the values in q in the same units. If method == 'mcmc', either q_err or q_cov must be passed. If method == 'leastsq' and neither q_err nor q_cov are passed, the absolute different between data points and fit is minimized. In this case, the returned chi2 is in units of absolute difference, meaning its value will depend on the units of q.
q_cov: array_like: Optional; the covariance matrix of the elements in q, as a 2-dimensional numpy array. This array must have dimensions of q**2 and be in units of the square of the units of q. If q_cov is passed, q_err is ignored since the diagonal elements of q_cov correspond to q_err**2.
method: str: The fitting method; can be leastsq for a least-squares minimization of mcmc for a Markov-Chain Monte Carlo.
mask: array_like: Optional; a numpy array of booleans that has the same length as the variables vector of the density profile class. Only variables where mask == True are varied in the fit, all others are kept constant. Important: this has to be a numpy array rather than a list.
verbose: bool / int: If true, output information about the fitting process. The flag can also be set as a number, where 1 has the same effect as True, and 2 outputs large amounts of information such as the fit parameters at each iteration.
tolerance: float: Only active when method == 'leastsq'. The accuracy to which the best-fit parameters are found.
maxfev: int: Only active when method == 'leastsq'. The maximum number of function evaluations before the fit is aborted. If None, the default value of the scipy least_squares function is used.
leastsq_algorithm: str: Only active when method == 'leastsq'. Can be any of the methods accepted by the scipy least_squares() function. Default is trf (trust region reflective), which works well for profile fits.
use_legacy_leastsq: bool: Only active when method == 'leastsq'. If True, this setting falls back to the old leastsq() in scipy rather than using the newer least_squares(). Should only be used for backward compatibility.
bounds: array_like: Only active when method == 'leastsq' and use_legacy_leastsq == False. If not None, this parameter must be an array of dimensions [2, N_par], giving two sets (lower and upper limits) of the fitted parameters (not the entire parameter array, if a mask is imposed). The limits must be given in linear space. If the parameters are fitted in log space (or some other transformation), that transformation is automatically applied to the bounds. For example, when fitting in log space (the default), the lower bounds must be positive, but the upper bounds can be np.inf.
initial_step: array_like: Only active when method == 'mcmc'. The MCMC samples (“walkers”) are initially distributed in a Gaussian around the initial guess. The width of the Gaussian is given by initial_step, either as an array of length N_par (giving the width of each Gaussian) or as a float number, in which case the width is set to initial_step times the initial value of the parameter.
nwalkers: int: Only active when method == 'mcmc'. The number of MCMC samplers that are run in parallel.
random_seed: int: Only active when method == 'mcmc'. If random_seed is not None, it is used to initialize the random number generator. This can be useful for reproducing results.
convergence_step: int: Only active when method == 'mcmc'. The convergence criteria are computed every convergence_step steps (and output is printed if verbose == True).
converged_GR: float: Only active when method == 'mcmc'. The maximum difference between different chains, according to the Gelman-Rubin criterion. Once the GR indicator is lower than this number in all parameters, the chain is ended. Setting this number too low leads to very long runtimes, but setting it too high can lead to inaccurate results.
best_fit: str: Only active when method == 'mcmc'. This parameter determines whether the mean or median value of the likelihood distribution is used as the output parameter set.
output_every_n: int: Only active when method == 'mcmc'. This parameter determines how frequently the MCMC chain outputs information. Only effective if verbose == True.

Returns

results: dict

A dictionary bundling the various fit results. Regardless of the fitting method, the dictionary always contains the following entries:

x: array_like: The best-fit result vector. If mask is passed, this vector only contains those variables that were varied in the fit.
q_fit: array_like: The fitted profile at the radii of the data points; has the same units as q and the same dimensions as r.
chi2: float: The chi^2 of the best-fit profile. If a covariance matrix was passed, the covariances are taken into account. If no uncertainty was passed at all, chi2 is in units of absolute difference, meaning its value will depend on the units of q.
ndof: int: The number of degrees of freedom, i.e. the number of fitted data points minus the number of free parameters.
chi2_ndof: float: The chi^2 per degree of freedom.

If method == 'leastsq', the dictionary additionally contains the entries returned by scipy.optimize.leastsq as well as the following:

nfev: int: The number of function calls used in the fit.
x_err: array_like: An array of dimensions [2, nparams] which contains an estimate of the lower and upper uncertainties on the fitted parameters. These uncertainties are computed from the covariance matrix estimated by the fitter. Please note that this estimate does not exactly correspond to a 68% likelihood. In order to get more statistically meaningful uncertainties, please use the MCMC samples instead of least-squares. In some cases, the fitter fails to return a covariance matrix, in which case x_err is None.

If method == 'mcmc', the dictionary contains the following entries:

x_initial: array_like: The initial positions of the walkers, in an array of dimensions [nwalkers, nparams].
chain_full: array_like: A numpy array of dimensions [n_independent_samples, nparams] with the parameters at each step in the chain. In this thin chain, only every nth step is output, where n is the auto-correlation time, meaning that the samples in this chain are truly independent.
chain_thin: array_like: Like the thin chain, but including all steps. Thus, the samples in this chain are not indepedent from each other. However, the full chain often gives better plotting results.
R: array_like: A numpy array containing the GR indicator at each step when it was saved.
x_mean: array_like: The mean of the chain for each parameter; has length nparams.
x_median: array_like: The median of the chain for each parameter; has length nparams.
x_stddev: array_like: The standard deviation of the chain for each parameter; has length nparams.
x_percentiles: array_like: The lower and upper values of each parameter that contain a certain percentile of the probability; has dimensions [n_percentages, 2, nparams] where the second dimension contains the lower/upper values.

getParameterArray(mask=None)

Returns an array of the profile parameters.

The profile parameters are internally stored in an ordered dictionary. For some applications (e.g., fitting), a simple array is more appropriate.

Parameters

mask: array_like: Optional; must be a numpy array (not a list) of booleans, with the same length as the parameter vector of the profile class (profile.N_par). Only those parameters that correspond to True values are returned.

Returns

par: array_like: A numpy array with the profile’s parameter values.

setNativeParameters(M, c, z, mdef, selected_by=None, Gamma=None, **kwargs)

Set the native parameters from mass and concentration (and optionally others).

Parameters

M: float: Spherical overdensity mass in \(M_{\odot}/h\).
c: float: The concentration, \(c = R / r_{\rm s}\), corresponding to the given halo mass and mass definition.
z: float: Redshift
mdef: str: The mass definition in which M and c are given. See Halo Mass Definitions for details.
selected_by: str: The halo sample to which this profile refers can be selected mass M or by accretion rate Gamma.
Gamma: float: The mass accretion rate as defined in DK14. This parameter only needs to be passed if selected_by == 'Gamma'. See comments in deriveParameters() function above.

setParameterArray(pars, mask=None)

Set the profile parameters from an array.

Parameters

pars: array_like: The new parameter array.
mask: array_like: Optional; must be a numpy array (not a list) of booleans, with the same length as the parameter vector of the profile class (profile.N_par). If passed, only those parameters that correspond to True values are set (meaning the pars parameter must be shorter than profile.N_par).

surfaceDensity(r, interpolate=True, accuracy=0.0001, max_r_interpolate=100000000.0, max_r_integrate=1e+20)

The projected surface density at radius r.

The surface density is computed by projecting the 3D density along the line of sight,

\[\Sigma(R) = 2 \int_R^{\infty} \frac{r \rho(r)}{\sqrt{r^2-R^2}} dr\]

Parameters

r: array_like: Radius in physical kpc/h; can be a number or a numpy array.
interpolate: bool: Use an interpolation table for density during the integration. This should make the evaluation somewhat faster, depending on how large the radius array is.
accuracy: float: The minimum accuracy of the integration.
max_r_interpolate: float: The maximum radius in physical kpc/h to which the density profile is integrated when using interpolating density.
max_r_integrate: float: The maximum radius in physical kpc/h to which the density profile is integrated when using exact densities.

Returns

Sigma: array_like: The surface density at radius r, in physical \(M_{\odot} h/{\rm kpc}^2\); has the same dimensions as r.

surfaceDensityInner(r, interpolate=True, accuracy=0.0001, max_r_interpolate=100000000.0, max_r_integrate=1e+20)

The projected surface density at radius r due to the inner profile.

Parameters

r: array_like: Radius in physical kpc/h; can be a number or a numpy array.
interpolate: bool: Use an interpolation table for density during the integration. This should make the evaluation somewhat faster, depending on how large the radius array is.
accuracy: float: The minimum accuracy of the integration.
max_r_interpolate: float: The maximum radius in physical kpc/h to which the density profile is integrated when using interpolating density.
max_r_integrate: float: The maximum radius in physical kpc/h to which the density profile is integrated when using exact densities.

Returns

Sigma: array_like: The surface density at radius r, in physical \(M_{\odot} h/{\rm kpc}^2\); has the same dimensions as r.

surfaceDensityOuter(r, interpolate=True, accuracy=0.0001, max_r_interpolate=100000000.0, max_r_integrate=1e+20)

The projected surface density at radius r due to the outer profile.

Parameters

r: array_like: Radius in physical kpc/h; can be a number or a numpy array.
interpolate: bool: Use an interpolation table for density during the integration. This should make the evaluation somewhat faster, depending on how large the radius array is.
accuracy: float: The minimum accuracy of the integration.
max_r_interpolate: float: The maximum radius in physical kpc/h to which the density profile is integrated when using interpolating density.
max_r_integrate: float: The maximum radius in physical kpc/h to which the density profile is integrated when using exact densities.

Returns

Sigma: array_like: The surface density at radius r, in physical \(M_{\odot} h/{\rm kpc}^2\); has the same dimensions as r.

update()

Update the profile object after a change in parameters or cosmology.

updateR200m()

Update the internally stored R200m after a parameter change.

class halo.profile_diemer23.ModelBProfile(selected_by='M', Gamma=None, eta=0.1, **kwargs)

The Diemer 2023 (truncated exponential) density profile (Model B).

This version corrects a minor flaw in the default model: the logarithmic slope at the scale radius is not -2 in the default Model A. In this Model B, this condition is enforced at the cost of an extra term, which gradually adjusts the slope between the center (where it is still zero) and the scale radius, where it offsets the effect of the truncation term. However, this correction is usually very small (except for extreme values of beta or rt). Thus, Model A and Model B profiles are virtually the same for almost all parameters. Model B can be a little more stable in fits to profiles without a clear scale radius. The nuissance parameter is set to \(\eta = 0.1\) by default; it is not recommended to change this parameter or to adjust it in fits.

The redshift must always be passed to this constructor, regardless of whether the fundamental parameters or a mass and concentration are given.

Parameters

rhos: float: The central scale density, in physical \(M_{\odot} h^2 / {\rm kpc}^3\).
rs: float: The scale radius in physical kpc/h.
rt: float: The radius where the profile steepens, in physical kpc/h.
alpha: float: Determines how quickly the slope of the inner Einasto profile steepens.
beta: float: Sharpness of the steepening.
eta: float: Nuissance parameter that determines how quickly the slope approaches zero at the halo center.
M: float: Halo mass in \(M_{\odot}/h\).
c: float: Concentration in the same mass definition as M.
mdef: str: The mass definition to which M corresponds. See Halo Mass Definitions for details.
z: float: Redshift
selected_by: str: The halo sample to which this profile refers can be selected mass M or by accretion rate Gamma. This parameter influences how some of the fixed parameters in the profile are set, in particular those that describe the steepening term.
Gamma: float: The mass accretion rate over the past dynamical time, which is defined as the crossing time (see the dynamicalTime() function or Diemer 2017 for details). The definition in the DK14 profile is slightly different, but the definitions are close enough that they can be used interchangeably without great loss of accuracy. The Gamma parameter only needs to be passed if selected_by == 'Gamma'.

Methods

`MDelta`(z, mdef)	The spherical overdensity mass of a given mass definition.
`RDelta`(z, mdef)	The spherical overdensity radius of a given mass definition.
`RMDelta`(z, mdef)	The spherical overdensity radius and mass of a given mass definition.
`Rsteepest`([search_range])	The radius where the logarithmic slope of the density profile is steepest.
`Vmax`()	The maximum circular velocity, and the radius where it occurs.
`circularVelocity`(r)	The circular velocity, \(v_c \equiv \sqrt{GM(<r)/r}\).
`cumulativePdf`(r[, Rmax, z, mdef])	The cumulative distribution function of the profile.
`deltaSigma`(r[, interpolate, ...])	The excess surface density at radius r.
`deltaSigmaInner`(r[, interpolate, ...])	The excess surface density at radius r due to the inner profile.
`deltaSigmaOuter`(r[, interpolate, ...])	The excess surface density at radius r due to the outer profile.
`density`(r)	Density as a function of radius.
`densityDerivativeLin`(r)	The linear derivative of density, \(d \rho / dr\).
`densityDerivativeLinInner`(r)	The linear derivative of the inner density, \(d \rho_{\rm inner} / dr\).
`densityDerivativeLinOuter`(r)	The linear derivative of the outer density, \(d \rho_{\rm outer} / dr\).
`densityDerivativeLog`(r)	The logarithmic derivative of density, \(d \log(\rho) / d \log(r)\).
`densityDerivativeLogInner`(r)	The logarithmic derivative of the inner density, \(d \log(\rho_{\rm inner}) / d \log(r)\).
`densityDerivativeLogOuter`(r)	The logarithmic derivative of the outer density, \(d \log(\rho_{\rm outer}) / d \log(r)\).
`densityInner`(r)	Density of the inner profile as a function of radius.
`densityOuter`(r)	Density of the outer profile as a function of radius.
`deriveParameters`(selected_by[, nu200m, z, Gamma])	Calibration of the parameters \(\alpha\), \(\beta\), and \(r_{\rm t}\).
`enclosedMass`(r[, accuracy])	The mass enclosed within radius r.
`enclosedMassInner`(r[, accuracy])	The mass enclosed within radius r due to the inner profile term.
`enclosedMassOuter`(r[, accuracy])	The mass enclosed within radius r due to the outer profile term.
`fit`(r, q, quantity[, q_err, q_cov, method, ...])	Fit the density, mass, surface density, or DeltaSigma profile to a given set of data points.
`getParameterArray`([mask])	Returns an array of the profile parameters.
`setNativeParameters`(M, c, z, mdef[, eta, ...])	Set the native Diemer22 parameters from mass and concentration (and optionally others).
`setParameterArray`(pars[, mask])	Set the profile parameters from an array.
`surfaceDensity`(r[, interpolate, accuracy, ...])	The projected surface density at radius r.
`surfaceDensityInner`(r[, interpolate, ...])	The projected surface density at radius r due to the inner profile.
`surfaceDensityOuter`(r[, interpolate, ...])	The projected surface density at radius r due to the outer profile.
`update`()	Update the profile object after a change in parameters or cosmology.
`updateR200m`()	Update the internally stored R200m after a parameter change.

setNativeParameters(M, c, z, mdef, eta=0.1, selected_by=None, Gamma=None, **kwargs)

Set the native Diemer22 parameters from mass and concentration (and optionally others).

The D22 profile has five free parameters, which are set by this function. The mass and concentration must be given as \(M_{\rm 200m}\) and \(c_{\rm 200m}\). Other mass definitions demand iteration, which can be achieved with the initialization routine in the parent class. This function ignores the presence of outer profiles.

Parameters

M: float: Spherical overdensity mass in \(M_{\odot}/h\).
c: float: The concentration, \(c = R / r_{\rm s}\), corresponding to the given halo mass and mass definition.
z: float: Redshift
mdef: str: The mass definition in which M and c are given. See Halo Mass Definitions for details.
eta: float: Nuissance parameter that determines how quickly the slope approaches zero at the halo center.
selected_by: str: The halo sample to which this profile refers can be selected mass M or by accretion rate Gamma.
Gamma: float: The mass accretion rate as defined in DK14. This parameter only needs to be passed if selected_by == 'Gamma'.

densityInner(r)

Density of the inner profile as a function of radius.

Parameters

r: array_like: Radius in physical kpc/h; can be a number or a numpy array.

Returns

density: array_like: Density in physical \(M_{\odot} h^2 / {\rm kpc}^3\); has the same dimensions as r.

densityDerivativeLogInner(r)

The logarithmic derivative of the inner density, \(d \log(\rho_{\rm inner}) / d \log(r)\).

Parameters

r: array_like: Radius in physical kpc/h; can be a number or a numpy array.

Returns

derivative: array_like: The dimensionless logarithmic derivative; has the same dimensions as r.

MDelta(z, mdef)

The spherical overdensity mass of a given mass definition.

Parameters

z: float: Redshift
mdef: str: The mass definition for which the spherical overdensity mass is computed. See Halo Mass Definitions for details.

Returns

M: float: Spherical overdensity mass in \(M_{\odot} /h\).

See also

RDelta: The spherical overdensity radius of a given mass definition.
RMDelta: The spherical overdensity radius and mass of a given mass definition.

RDelta(z, mdef)

The spherical overdensity radius of a given mass definition.

Parameters

z: float: Redshift
mdef: str: The mass definition for which the spherical overdensity radius is computed. See Halo Mass Definitions for details.

Returns

R: float: Spherical overdensity radius in physical kpc/h.

See also

MDelta: The spherical overdensity mass of a given mass definition.
RMDelta: The spherical overdensity radius and mass of a given mass definition.

RMDelta(z, mdef)

The spherical overdensity radius and mass of a given mass definition.

This is a wrapper for the RDelta() and MDelta() functions which returns both radius and mass.

Parameters

z: float: Redshift
mdef: str: The mass definition for which the spherical overdensity mass is computed. See Halo Mass Definitions for details.

Returns

R: float: Spherical overdensity radius in physical kpc/h.
M: float: Spherical overdensity mass in \(M_{\odot} /h\).

See also

RDelta: The spherical overdensity radius of a given mass definition.
MDelta: The spherical overdensity mass of a given mass definition.

Rsteepest(search_range=10.0)

The radius where the logarithmic slope of the density profile is steepest.

Parameters

search_range: float: When searching for the radius of steepest slope, search within this factor of \(R_{\rm 200m}\), \(r_{\rm s}\), or another initial guess, which is determined automatically.

Returns

Rsteep: float: The radius where the slope is steepest, in physical kpc/h.

Vmax()

The maximum circular velocity, and the radius where it occurs.

Returns

vmax: float: The maximum circular velocity in km / s.
rmax: float: The radius where fmax occurs, in physical kpc/h.

See also

circularVelocity: The circular velocity, \(v_c \equiv \sqrt{GM(<r)/r}\).

circularVelocity(r)

The circular velocity, \(v_c \equiv \sqrt{GM(<r)/r}\).

Parameters

r: array_like: Radius in physical kpc/h; can be a number or a numpy array.

Returns

vc: float: The circular velocity in km / s; has the same dimensions as r.

See also

Vmax: The maximum circular velocity, and the radius where it occurs.

cumulativePdf(r, Rmax=None, z=None, mdef=None)

The cumulative distribution function of the profile.

Parameters

r: array_like: Radius in physical kpc/h; can be a number or a numpy array.
Rmax: float: The radius where to cut off the profile in physical kpc/h.
z: float: Redshift
mdef: str: The radius definition for the cut-off radius. See Halo Mass Definitions for details.

Returns

pdf: array_like: The probability for mass to lie within radius r; has the same dimensions as r.

deltaSigma(r, interpolate=True, interpolate_surface_density=True, accuracy=0.0001, min_r_interpolate=1e-06, max_r_interpolate=100000000.0, max_r_integrate=1e+20)

The excess surface density at radius r.

This quantity is useful in weak lensing studies, and is defined as \(\Delta\Sigma(R) = \Sigma(<R)-\Sigma(R)\) where \(\Sigma(<R)\) is the averaged surface density within R weighted by area,

\[\Delta\Sigma(R) = \frac{1}{\pi R^2} \int_0^{R} 2 \pi r \Sigma(r) dr - \Sigma(R)\]

Parameters

r: array_like: Radius in physical kpc/h; can be a number or a numpy array.
interpolate: bool: Use an interpolation table for the surface density during the integration. This can speed up the evaluation significantly, as the surface density can be expensive to evaluate.
interpolate_surface_density: bool: Use an interpolation table for density during the computation of the surface density. This should make the evaluation somewhat faster, but can fail for some density terms which are negative at particular radii.
accuracy: float: The minimum accuracy of the integration (used both to compute the surface density and average it to get DeltaSigma).
min_r_interpolate: float: The minimum radius in physical kpc/h from which the surface density profile is averaged.
max_r_interpolate: float: The maximum radius in physical kpc/h to which the density profile is integrated when using interpolating density.
max_r_integrate: float: The maximum radius in physical kpc/h to which the density profile is integrated when using exact densities.

Returns

DeltaSigma: array_like: The excess surface density at radius r, in physical \(M_{\odot} h/{\rm kpc}^2\); has the same dimensions as r.

deltaSigmaInner(r, interpolate=True, interpolate_surface_density=True, accuracy=0.0001, min_r_interpolate=1e-06, max_r_interpolate=100000000.0, max_r_integrate=1e+20)

The excess surface density at radius r due to the inner profile.

Parameters

r: array_like: Radius in physical kpc/h; can be a number or a numpy array.
interpolate: bool: Use an interpolation table for the surface density during the integration. This can speed up the evaluation significantly, as the surface density can be expensive to evaluate.
interpolate_surface_density: bool: Use an interpolation table for density during the computation of the surface density. This should make the evaluation somewhat faster, but can fail for some density terms which are negative at particular radii.
accuracy: float: The minimum accuracy of the integration (used both to compute the surface density and average it to get DeltaSigma).
min_r_interpolate: float: The minimum radius in physical kpc/h from which the surface density profile is averaged.
max_r_interpolate: float: The maximum radius in physical kpc/h to which the density profile is integrated when using interpolating density.
max_r_integrate: float: The maximum radius in physical kpc/h to which the density profile is integrated when using exact densities.

Returns

DeltaSigma: array_like: The excess surface density at radius r, in physical \(M_{\odot} h/{\rm kpc}^2\); has the same dimensions as r.

deltaSigmaOuter(r, interpolate=True, interpolate_surface_density=True, accuracy=0.0001, min_r_interpolate=1e-06, max_r_interpolate=100000000.0, max_r_integrate=1e+20)

The excess surface density at radius r due to the outer profile.

Parameters

r: array_like: Radius in physical kpc/h; can be a number or a numpy array.
interpolate: bool: Use an interpolation table for the surface density during the integration. This can speed up the evaluation significantly, as the surface density can be expensive to evaluate.
interpolate_surface_density: bool: Use an interpolation table for density during the computation of the surface density. This should make the evaluation somewhat faster, but can fail for some density terms which are negative at particular radii.
accuracy: float: The minimum accuracy of the integration (used both to compute the surface density and average it to get DeltaSigma).
min_r_interpolate: float: The minimum radius in physical kpc/h from which the surface density profile is averaged.
max_r_interpolate: float: The maximum radius in physical kpc/h to which the density profile is integrated when using interpolating density.
max_r_integrate: float: The maximum radius in physical kpc/h to which the density profile is integrated when using exact densities.

Returns

DeltaSigma: array_like: The excess surface density at radius r, in physical \(M_{\odot} h/{\rm kpc}^2\); has the same dimensions as r.

density(r)

Density as a function of radius.

Parameters

r: array_like: Radius in physical kpc/h; can be a number or a numpy array.

Returns

density: array_like: Density in physical \(M_{\odot} h^2 / {\rm kpc}^3\); has the same dimensions as r.

densityDerivativeLin(r)

The linear derivative of density, \(d \rho / dr\).

This function should generally not be overwritten by child classes since it handles the general case of adding up the contributions from the inner and outer terms.

Parameters

r: array_like: Radius in physical kpc/h; can be a number or a numpy array.

Returns

derivative: array_like: The linear derivative in physical \(M_{\odot} h / {\rm kpc}^2\); has the same dimensions as r.

densityDerivativeLinInner(r)

The linear derivative of the inner density, \(d \rho_{\rm inner} / dr\).

For the truncated exponential profile, the logarithmic derivative is much easier to evaluate. Thus, this function converts the logarithmic to the linear derivative.

Parameters

r: array_like: Radius in physical kpc/h; can be a number or a numpy array.

Returns

derivative: array_like: The linear derivative in physical \(M_{\odot} h / {\rm kpc}^2\); has the same dimensions as r.

densityDerivativeLinOuter(r)

The linear derivative of the outer density, \(d \rho_{\rm outer} / dr\).

This function should generally not be overwritten by child classes since it handles the general case of adding up the contributions from all outer profile terms.

Parameters

r: array_like: Radius in physical kpc/h; can be a number or a numpy array.

Returns

derivative: array_like: The linear derivative in physical \(M_{\odot} h / {\rm kpc}^2\); has the same dimensions as r.

densityDerivativeLog(r)

The logarithmic derivative of density, \(d \log(\rho) / d \log(r)\).

This function should generally not be overwritten by child classes since it handles the general case of adding up the contributions from the inner and outer profile terms.

Parameters

r: array_like: Radius in physical kpc/h; can be a number or a numpy array.

Returns

derivative: array_like: The dimensionless logarithmic derivative; has the same dimensions as r.

densityDerivativeLogOuter(r)

The logarithmic derivative of the outer density, \(d \log(\rho_{\rm outer}) / d \log(r)\).

This function should generally not be overwritten by child classes since it handles the general case of adding up the contributions from outer profile terms.

Parameters

r: array_like: Radius in physical kpc/h; can be a number or a numpy array.

Returns

derivative: array_like: The dimensionless logarithmic derivative; has the same dimensions as r.

densityOuter(r)

Density of the outer profile as a function of radius.

This function should generally not be overwritten by child classes since it handles the general case of adding up the contributions from all outer profile terms.

Parameters

r: array_like: Radius in physical kpc/h; can be a number or a numpy array.

Returns

density: array_like: Density in physical \(M_{\odot} h^2 / {\rm kpc}^3\); has the same dimensions as r.

static deriveParameters(selected_by, nu200m=None, z=None, Gamma=None)

Calibration of the parameters \(\alpha\), \(\beta\), and \(r_{\rm t}\).

Parameters

selected_by: str: The halo sample to which this profile refers can be selected mass M or by accretion rate Gamma.
nu200m: float: The peak height of the halo for which the parameters are to be calibrated. This parameter only needs to be passed if selected_by == 'M'.
z: float: Redshift
Gamma: float: The mass accretion rate over the past dynamical time, which is defined as the crossing time (see the dynamicalTime() function or Diemer 2017 for details). The definition in the DK14 profile is slightly different, but the definitions are close enough that they can be used interchangeably without great loss of accuracy. The Gamma parameter only needs to be passed if selected_by == 'Gamma'.

Returns

alpha: float: The Einasto steepening parameter.
beta: float: The steepening of the truncation term.
rt_R200m: float: The truncation radius in units of R200m.

enclosedMass(r, accuracy=1e-06)

The mass enclosed within radius r.

Parameters

r: array_like: Radius in physical kpc/h; can be a number or a numpy array.
accuracy: float: The minimum accuracy of the integration.

Returns

M: array_like: The mass enclosed within radius r, in \(M_{\odot}/h\); has the same dimensions as r.

enclosedMassInner(r, accuracy=1e-06)

The mass enclosed within radius r due to the inner profile term.

This function should be overwritten by child classes if an analytical expression exists.

Parameters

r: array_like: Radius in physical kpc/h; can be a number or a numpy array.
accuracy: float: The minimum accuracy of the integration.

Returns

M: array_like: The mass enclosed within radius r, in \(M_{\odot}/h\); has the same dimensions as r.

enclosedMassOuter(r, accuracy=1e-06)

The mass enclosed within radius r due to the outer profile term.

Parameters

r: array_like: Radius in physical kpc/h; can be a number or a numpy array.
accuracy: float: The minimum accuracy of the integration.

Returns

M: array_like: The mass enclosed within radius r, in \(M_{\odot}/h\); has the same dimensions as r.

Fit the density, mass, surface density, or DeltaSigma profile to a given set of data points.

The user can choose to vary only a sub-set of the profile parameters through the mask parameter. The current parameters of the profile instance serve as an initial guess.

Parameters

r: array_like: The radii of the data points, in physical kpc/h.
q: array_like: The data to fit; can either be density in physical \(M_{\odot} h^2 / {\rm kpc}^3\), enclosed mass in \(M_{\odot} /h\), or surface density in physical \(M_{\odot} h/{\rm kpc}^2\). Must have the same dimensions as r.
quantity: str: Indicates which quantity is given in as input in q, can be rho, M, Sigma, or DeltaSigma.
q_err: array_like: Optional; the uncertainty on the values in q in the same units. If method == 'mcmc', either q_err or q_cov must be passed. If method == 'leastsq' and neither q_err nor q_cov are passed, the absolute different between data points and fit is minimized. In this case, the returned chi2 is in units of absolute difference, meaning its value will depend on the units of q.
q_cov: array_like: Optional; the covariance matrix of the elements in q, as a 2-dimensional numpy array. This array must have dimensions of q**2 and be in units of the square of the units of q. If q_cov is passed, q_err is ignored since the diagonal elements of q_cov correspond to q_err**2.
method: str: The fitting method; can be leastsq for a least-squares minimization of mcmc for a Markov-Chain Monte Carlo.
mask: array_like: Optional; a numpy array of booleans that has the same length as the variables vector of the density profile class. Only variables where mask == True are varied in the fit, all others are kept constant. Important: this has to be a numpy array rather than a list.
verbose: bool / int: If true, output information about the fitting process. The flag can also be set as a number, where 1 has the same effect as True, and 2 outputs large amounts of information such as the fit parameters at each iteration.
tolerance: float: Only active when method == 'leastsq'. The accuracy to which the best-fit parameters are found.
maxfev: int: Only active when method == 'leastsq'. The maximum number of function evaluations before the fit is aborted. If None, the default value of the scipy least_squares function is used.
leastsq_algorithm: str: Only active when method == 'leastsq'. Can be any of the methods accepted by the scipy least_squares() function. Default is trf (trust region reflective), which works well for profile fits.
use_legacy_leastsq: bool: Only active when method == 'leastsq'. If True, this setting falls back to the old leastsq() in scipy rather than using the newer least_squares(). Should only be used for backward compatibility.
bounds: array_like: Only active when method == 'leastsq' and use_legacy_leastsq == False. If not None, this parameter must be an array of dimensions [2, N_par], giving two sets (lower and upper limits) of the fitted parameters (not the entire parameter array, if a mask is imposed). The limits must be given in linear space. If the parameters are fitted in log space (or some other transformation), that transformation is automatically applied to the bounds. For example, when fitting in log space (the default), the lower bounds must be positive, but the upper bounds can be np.inf.
initial_step: array_like: Only active when method == 'mcmc'. The MCMC samples (“walkers”) are initially distributed in a Gaussian around the initial guess. The width of the Gaussian is given by initial_step, either as an array of length N_par (giving the width of each Gaussian) or as a float number, in which case the width is set to initial_step times the initial value of the parameter.
nwalkers: int: Only active when method == 'mcmc'. The number of MCMC samplers that are run in parallel.
random_seed: int: Only active when method == 'mcmc'. If random_seed is not None, it is used to initialize the random number generator. This can be useful for reproducing results.
convergence_step: int: Only active when method == 'mcmc'. The convergence criteria are computed every convergence_step steps (and output is printed if verbose == True).
converged_GR: float: Only active when method == 'mcmc'. The maximum difference between different chains, according to the Gelman-Rubin criterion. Once the GR indicator is lower than this number in all parameters, the chain is ended. Setting this number too low leads to very long runtimes, but setting it too high can lead to inaccurate results.
best_fit: str: Only active when method == 'mcmc'. This parameter determines whether the mean or median value of the likelihood distribution is used as the output parameter set.
output_every_n: int: Only active when method == 'mcmc'. This parameter determines how frequently the MCMC chain outputs information. Only effective if verbose == True.

Returns

results: dict

A dictionary bundling the various fit results. Regardless of the fitting method, the dictionary always contains the following entries:

x: array_like: The best-fit result vector. If mask is passed, this vector only contains those variables that were varied in the fit.
q_fit: array_like: The fitted profile at the radii of the data points; has the same units as q and the same dimensions as r.
chi2: float: The chi^2 of the best-fit profile. If a covariance matrix was passed, the covariances are taken into account. If no uncertainty was passed at all, chi2 is in units of absolute difference, meaning its value will depend on the units of q.
ndof: int: The number of degrees of freedom, i.e. the number of fitted data points minus the number of free parameters.
chi2_ndof: float: The chi^2 per degree of freedom.

If method == 'leastsq', the dictionary additionally contains the entries returned by scipy.optimize.leastsq as well as the following:

nfev: int: The number of function calls used in the fit.
x_err: array_like: An array of dimensions [2, nparams] which contains an estimate of the lower and upper uncertainties on the fitted parameters. These uncertainties are computed from the covariance matrix estimated by the fitter. Please note that this estimate does not exactly correspond to a 68% likelihood. In order to get more statistically meaningful uncertainties, please use the MCMC samples instead of least-squares. In some cases, the fitter fails to return a covariance matrix, in which case x_err is None.

If method == 'mcmc', the dictionary contains the following entries:

x_initial: array_like: The initial positions of the walkers, in an array of dimensions [nwalkers, nparams].
chain_full: array_like: A numpy array of dimensions [n_independent_samples, nparams] with the parameters at each step in the chain. In this thin chain, only every nth step is output, where n is the auto-correlation time, meaning that the samples in this chain are truly independent.
chain_thin: array_like: Like the thin chain, but including all steps. Thus, the samples in this chain are not indepedent from each other. However, the full chain often gives better plotting results.
R: array_like: A numpy array containing the GR indicator at each step when it was saved.
x_mean: array_like: The mean of the chain for each parameter; has length nparams.
x_median: array_like: The median of the chain for each parameter; has length nparams.
x_stddev: array_like: The standard deviation of the chain for each parameter; has length nparams.
x_percentiles: array_like: The lower and upper values of each parameter that contain a certain percentile of the probability; has dimensions [n_percentages, 2, nparams] where the second dimension contains the lower/upper values.

getParameterArray(mask=None)

Returns an array of the profile parameters.

The profile parameters are internally stored in an ordered dictionary. For some applications (e.g., fitting), a simple array is more appropriate.

Parameters

mask: array_like: Optional; must be a numpy array (not a list) of booleans, with the same length as the parameter vector of the profile class (profile.N_par). Only those parameters that correspond to True values are returned.

Returns

par: array_like: A numpy array with the profile’s parameter values.

setParameterArray(pars, mask=None)

Set the profile parameters from an array.

Parameters

pars: array_like: The new parameter array.
mask: array_like: Optional; must be a numpy array (not a list) of booleans, with the same length as the parameter vector of the profile class (profile.N_par). If passed, only those parameters that correspond to True values are set (meaning the pars parameter must be shorter than profile.N_par).

surfaceDensity(r, interpolate=True, accuracy=0.0001, max_r_interpolate=100000000.0, max_r_integrate=1e+20)

The projected surface density at radius r.

The surface density is computed by projecting the 3D density along the line of sight,

\[\Sigma(R) = 2 \int_R^{\infty} \frac{r \rho(r)}{\sqrt{r^2-R^2}} dr\]

Parameters

r: array_like: Radius in physical kpc/h; can be a number or a numpy array.
interpolate: bool: Use an interpolation table for density during the integration. This should make the evaluation somewhat faster, depending on how large the radius array is.
accuracy: float: The minimum accuracy of the integration.
max_r_interpolate: float: The maximum radius in physical kpc/h to which the density profile is integrated when using interpolating density.
max_r_integrate: float: The maximum radius in physical kpc/h to which the density profile is integrated when using exact densities.

Returns

Sigma: array_like: The surface density at radius r, in physical \(M_{\odot} h/{\rm kpc}^2\); has the same dimensions as r.

surfaceDensityInner(r, interpolate=True, accuracy=0.0001, max_r_interpolate=100000000.0, max_r_integrate=1e+20)

The projected surface density at radius r due to the inner profile.

Parameters

r: array_like: Radius in physical kpc/h; can be a number or a numpy array.
interpolate: bool: Use an interpolation table for density during the integration. This should make the evaluation somewhat faster, depending on how large the radius array is.
accuracy: float: The minimum accuracy of the integration.
max_r_interpolate: float: The maximum radius in physical kpc/h to which the density profile is integrated when using interpolating density.
max_r_integrate: float: The maximum radius in physical kpc/h to which the density profile is integrated when using exact densities.

Returns

Sigma: array_like: The surface density at radius r, in physical \(M_{\odot} h/{\rm kpc}^2\); has the same dimensions as r.

surfaceDensityOuter(r, interpolate=True, accuracy=0.0001, max_r_interpolate=100000000.0, max_r_integrate=1e+20)

The projected surface density at radius r due to the outer profile.

Parameters

r: array_like: Radius in physical kpc/h; can be a number or a numpy array.
interpolate: bool: Use an interpolation table for density during the integration. This should make the evaluation somewhat faster, depending on how large the radius array is.
accuracy: float: The minimum accuracy of the integration.
max_r_interpolate: float: The maximum radius in physical kpc/h to which the density profile is integrated when using interpolating density.
max_r_integrate: float: The maximum radius in physical kpc/h to which the density profile is integrated when using exact densities.

Returns

Sigma: array_like: The surface density at radius r, in physical \(M_{\odot} h/{\rm kpc}^2\); has the same dimensions as r.

update()

Update the profile object after a change in parameters or cosmology.

updateR200m()

Update the internally stored R200m after a parameter change.

Diemer 2023 (truncated exponential) profile

Previous topic

Diemer & Kravtsov 2014 profile

Next topic

Concentration

Diemer 2023 (truncated exponential) profile

Basics

Adding an infalling profile

Model variant with correction at scale radius

Module reference

Table of Contents

Previous topic

Next topic