Source code notes on NFWProfile and NFWPhaseSpace

This section of the documentation provides background material and detailed implementation notes on the functions and methods of the NFWProfile and NFWPhaseSpace models.

Outline

We will start in Class structure of the NFWPhaseSpace model with an overview of the orthogonal mix-in class design. The Modeling the NFW Spatial Profile section covers the mathematics of NFW spatial profiles and their Monte Carlo realizations, including explicit derivations of the exact form of all equations are they are implemented in code. The section Modeling the NFW Velocity Profile goes into the same level of detail but for the velocity profiles. The tutorial concludes in Monte Carlo realizations of the NFW profile by describing how the NFWPhaseSpace class can be used in both stand-alone fashion and as part of the Halotools mock-making framework to generate Monte Carlo realizations of points in NFW phase space.

Class structure of the NFWPhaseSpace model

The NFWPhaseSpace model is a container class for two independently-defined sets of behaviors:

1. Analytical descriptions of the spatial and velocity distribution of points within the halo

2. Monte Carlo methods for generating random realizations of points in phase space

For the first set of behaviors, the spatial distributions are modeled by NFWProfile, while the NFWPhaseSpace sub-class controls the velocities. The second set of behaviors is controlled entirely by MonteCarloGalProf, regardless of the details of the profile model.

Modeling the NFW Spatial Profile

The spatial profile of an NFW halo is modeled with the NFWProfile class, which is itself a sub-class of AnalyticDensityProf

Most of the functionality of the NFWProfile class derives from the behavior defined in the AnalyticDensityProf super-class. Here we only document the functionality and implementation that is unique to the NFWProfile model, and defer discussion of all super-class-derived behavior to the Source code notes on AnalyticDensityProf.

Halotools implementation of the NFW mass density

Equation 3 of Navarro, Frenk and White (1995), arXiv:9508025, gives the original definition of the NFW profile:

\[ \begin{align}\begin{aligned}\begin{split}\rho_{\rm NFW}(r) \equiv \rho_{\rm crit}\frac{\delta_{\rm c}}{r/r_{\rm s}(1 + r/r_{s})^{2}}, \\\end{split}\\\delta_{\rm c} \equiv \frac{200}{3}\frac{c^{3}}{{\rm ln}(1+c) - c/(1+c)}\end{aligned}\end{align} \]

In the above equation, the factor of \(200\rho_{\rm crit}\) reflects the 200m halo mass definition convention adopted in NFW95, and \(c \equiv R_{200}/r_{s}\) is the concentration parameter. In the notation adopted in Halo Mass and Radius Definitions, \(200\times\rho_{\rm crit} \equiv \Delta_{\rm thresh}\times\rho_{\rm ref} \equiv \rho_{\rm thresh}\), so we can reformulate the above expression in the following more general form:

\[ \begin{align}\begin{aligned}\begin{split}\rho_{\rm NFW}(r) = \rho_{\rm thresh}\frac{c^{3}/3g(c)}{r/r_{s}(1 + r/r_{s})^{2}}, \\\end{split}\\g(c) \equiv {\rm ln}(1+c) - c/(1+c).\end{aligned}\end{align} \]

If we now substitute \(r/r_{s} = cr/R_{\Delta}\) and define the scaled radius \(\tilde{r}\equiv r/R_{\Delta}\) the above expression can be further rewritten as:

\[\rho_{\rm NFW}(\tilde{r})/\rho_{\rm thresh} \equiv \tilde{\rho}_{\rm NFW}(\tilde{r}) = \frac{c^{3}/3g(c)}{c\tilde{r}(1 + c\tilde{r})^{2}}.\]

The above expression is the exact equation implemented in the dimensionless_mass_density method of the NFWProfile class. The quantity \(\rho_{\rm thresh}\) is calculated in Halotools using the density_threshold function, and \(g(c)\) is computed using the g method of the NFWProfile class.

For any sub-class of AnalyticDensityProf, once the dimensionless_mass_density method is defined, in principle all subsequent behavior is derived. In practice, if the associated integrals and derivatives can be computed analytically it is more efficient and numerically stable to implement the analytical results as over-rides of the super-class-defined methods. The subsections below derive the analytical equations used in all over-rides implemented in the NFWProfile class.

Derivation of the NFW cumulative mass PDF

The cumulative mass PDF, \(P_{\rm prof}(<\tilde{r})\), the cumulative probability of finding a randomly selected particle at a scaled-radius position less than \(\tilde{r}\), is defined as:

\[P_{\rm NFW}(<\tilde{r}) \equiv M_{\rm NFW}(<\tilde{r}) / M_{\Delta}.\]

In the above expression,

\[M_{\rm NFW}(<\tilde{r}) \equiv 4\pi\rho_{\rm thresh}\int_{0}^{\tilde{r}}d\tilde{r}' \tilde{r}'^{2}\tilde{\rho}_{\rm NFW}(\tilde{r}')\]

and

\[M_{\Delta} \equiv 4\pi\rho_{\rm thresh}\int_{0}^{1}d\tilde{r}' \tilde{r}'^{2}\tilde{\rho}_{\rm NFW}(\tilde{r}').\]

Plugging in the definition of \(\tilde{\rho}_{\rm NFW}\) and canceling the common pre-factors of \(4\pi\rho_{\rm thresh}c^{3}/3g(c)\) gives:

\[P_{\rm NFW}(<\tilde{r}) = \frac{\int_{0}^{\tilde{r}}d\tilde{r}' \tilde{r}'^{2}1/c\tilde{r}'(1 + c\tilde{r}')^{2}}{\int_{0}^{1}d\tilde{r}' \tilde{r}'^{2}1/c\tilde{r}'(1 + c\tilde{r}')^{2}}\]

Now we change integration variables \(\tilde{r}'\rightarrow c\tilde{r}'=y\):

\[P_{\rm NFW}(<\tilde{r}) = \frac{\int_{0}^{c\tilde{r}}dy\frac{y}{(1 + y)^{2}}}{\int_{0}^{1}dy\frac{y}{(1 + y)^{2}}}\]

and use the definition of \(g(x) \equiv {\rm ln}(1+x) - x/(1+x) = \int_{0}^{x}dy\frac{y}{(1+y)^{2}}\) to write the above expression as

\[P_{\rm NFW}(<\tilde{r}) = g(c\tilde{r}) / g(c)\]

The above equation is the exact expression used to calculate \(P_{\rm NFW}(<\tilde{r})\) via the cumulative_mass_PDF function.

Monte Carlo realizations of the NFW profile

Halotools uses Inverse Transform Sampling, a standard Monte Carlo technique, to produce random realizations of halo profiles. The basic idea of this technique is to draw a random uniform number, u, and intrepret u as the probability \(u = P(<r)\) of finding a point tracing an NFW radial profile interior to position r. The mapping between u and r is already implemented via the cumulative_mass_PDF function, so we only need to use this function to provide the inverse mapping. This we do numerically by tabulating \(P_{\rm NFW}(<\tilde{r})\) at a set of control points \(0<\tilde{r}<1\) and then using the scipy function InterpolatedUnivariateSpline. This technique is used ubiquitously throughout the package, and the interpolation is actually implemented using the custom_spline function, which is just a wrapper that customizes the edge case behavior of InterpolatedUnivariateSpline.

The simplest place in the code base to see where Inverse Transform Sampling gives Monte Carlo realizations of the NFW profile is in the mc_generate_nfw_radial_positions source code. Here the implementation is basically straightforward. Because NFW profiles are power laws, the interpolation is more stable when it is done in log-space.

Modeling the NFW Velocity Profile

The NFWPhaseSpace model solves for the velocity profile of satellite galaxies by making the following assumptions:

1. satellites trace the same spatial profile as their underlying gravitational potential well,

2. satellites are in virial equilibrium with their potential, and

3. satellite orbits are isotropic.

In any such system, the first moment of the radial velocity distribution of the satellites is zero (there is no net infall or outflow), and the second moment \(\sigma^{2}_{r}(r)\) can be calculated analytically by solving the Jeans equation, which under these assumptions takes the following form:

\[\sigma^{2}_{r}(r) = \frac{1}{\rho_{\rm sat}(r)}\int_{r}^{\infty}{\rm d}r\rho_{\rm sat}(r)\frac{{\rm d}\Phi(r)}{{\rm d}r},\]

In the above equation, \(\rho_{\rm sat}\) is the number density profile of the satellite galaxies and \(\Phi\) is the gravitational potential. For a justification of this simplification of the Jeans equation, see Jeans Equation Derivations.

The NFWPhaseSpace model assumes that \(\rho_{\rm sat} = \rho_{\rm NFW}\). So we can rewrite the above equation using the dimensionless quantities \(\tilde{r}\equiv r/R_{\Delta}\) and \(\tilde{\rho}_{\rm NFW}(\tilde{r}) \equiv \rho_{\rm NFW}(\tilde{r})/\rho_{\rm thresh}\) and canceling the common factors of \(\rho_{\rm thresh}\):

\[\sigma^{2}_{r}(\tilde{r}) = \frac{1}{\tilde{\rho}_{\rm NFW}(\tilde{r})}\int_{\tilde{r}}^{\infty}{\rm d}\tilde{r}\tilde{\rho}_{\rm NFW}(\tilde{r})\frac{{\rm d}\Phi(\tilde{r})}{{\rm d}\tilde{r}},\]

For any spherically symmetric gravitational potential,

\[ \begin{align}\begin{aligned}\begin{split}\Phi(x) = \frac{-GM_{\Delta}(<x)}{x} \\\end{split}\\\Rightarrow \frac{{\rm d}\Phi(\tilde{r})}{{\rm d}\tilde{r}} = \frac{GM_{\Delta}(<\tilde{r})}{\tilde{r}^{2}} \equiv \frac{V_{\rm circ}^{2}(\tilde{r})}{\tilde{r}} = V_{\rm vir}^{2}\frac{P_{\rm NFW}(<\tilde{r})}{\tilde{r}^{2}},\end{aligned}\end{align} \]

where in the second-to-last equality we have used the definition of \(V^{2}_{\rm circ}\), and in the last equality we have used the derivation provided in the Circular velocity profile section of the Source code notes on AnalyticDensityProf.

From the Derivation of the NFW cumulative mass PDF we have that \(P_{\rm NFW}(<\tilde{r}) = g(c\tilde{r}) / g(c)\), and using the expression for \(\tilde{\rho}_{\rm NFW}\) given in Halotools implementation of the NFW mass density, we can plug in these expressions to the above equation:

\[\sigma^{2}_{r}(\tilde{r}) = \frac{c\tilde{r}(1 + c\tilde{r})^{2}}{c^{3}/3g(c)}\int_{\tilde{r}}^{\infty}{\rm d}\tilde{r}\frac{c^{3}/3g(c)}{c\tilde{r}(1 + c\tilde{r})^{2}}V_{\rm vir}^{2}\frac{g(c\tilde{r})}{g(c)\tilde{r}^{2}}.\]

Canceling common factors of \(c^{3}/3g(c)\) and rearranging terms gives us:

\[\Rightarrow \sigma^{2}_{r}(\tilde{r}) = V_{\rm vir}^{2}\frac{c\tilde{r}(1 + c\tilde{r})^{2}}{g(c)}\int_{\tilde{r}}^{\infty}{\rm d}\tilde{r}\frac{g(c\tilde{r})}{c\tilde{r}^{3}(1 + c\tilde{r})^{2}}\]

Finally, we change integration variables \(\tilde{r}\rightarrow c\tilde{r}=y\) to give:

\[\Rightarrow \sigma^{2}_{r}(\tilde{r}) = V_{\rm vir}^{2}\frac{c^{2}\tilde{r}(1 + c\tilde{r})^{2}}{g(c)}\int_{c\tilde{r}}^{\infty}{\rm d}y\frac{g(y)}{y^{3}(1 + y)^{2}}\]

Defining the dimensionless radial velocity dispersion \(\tilde{\sigma}_{r}\equiv\sigma_{r}/V_{\rm vir}\), the above equation is the exact expression used in the dimensionless_radial_velocity_dispersion method of the NFWPhaseSpace class. The above expression is also the same expression appearing in Eq. 24 of More et al. (2008), arXiv:0807.4529, with the only differences being of notation: \(g(c) \leftrightarrow \mu(c)\) and \(c\tilde{r} \leftrightarrow r/r_{\rm s}\).

Monte Carlo realizations of the NFW profile

The NFWPhaseSpace model can either be used as a stand-alone class to generate an arbitrary number of points in NFW phase space, or as part of a composite galaxy-halo model that generates full-scale mock galaxy catalogs. We document each of these options in turn.

Stand-alone Monte Carlo realizations of NFW phase space

The mc_generate_nfw_phase_space_points method can be used to create an Astropy Table storing a collection of points in NFW phase space.

>>> from halotools.empirical_models import NFWPhaseSpace
>>> nfw = NFWPhaseSpace()
>>> data = nfw.mc_generate_nfw_phase_space_points(Ngals = 100, mass = 1e13, conc = 10)  

In the source code, the generation of these points happens in two steps. First, x, y, z points are drawn using the mc_halo_centric_pos method defined in the MonteCarloGalProf orthogonal mix-in class. Following the same methodology described in Monte Carlo realizations of the NFW profile, the mc_halo_centric_pos method uses inverse transform sampling together with the cumulative_mass_PDF function to draw random realizations of dimensionless NFW profile radii; these dimensionless radii are then scaled according to the halo mass and radius definition selected by the keyword arguments passed to the NFWPhaseSpace constructor. See the Monte Carlo Generation of Spatial Positions section of the Source code notes on MonteCarloGalProf for a detailed explanation of how this method works.

Once dimensionless radial positions have been drawn, the mc_generate_nfw_phase_space_points method passes these positions to the mc_radial_velocity method of the MonteCarloGalProf orthogonal mix-in class. This method works differently than the mc_halo_centric_pos method: the mc_radial_velocity method does not use inverse transform sampling. This is because radial velocity distributions are assumed to be Gaussian, and there is an optimized function numpy.random.normal in the Numpy code base for directly drawing from a Gaussian distribution. A Gaussian distribution is specified by its first two moments: the first moment is centered at the velocity of the host halo, and the second moment is calculated using the dimensionless_radial_velocity_dispersion method described in the previous section.

In practice, for performance reasons the mc_radial_velocity method actually uses a lookup table tabulated from dimensionless_radial_velocity_dispersion rather than the actual dimensionless_radial_velocity_dispersion method itself. For further information concerning this detail, see the Source code notes on MonteCarloGalProf.

At this point, random positions and velocities have been drawn for the satellites and the mc_generate_nfw_phase_space_points method bundles these arrays into an Astropy Table and returns the result.

Making mocks with NFWPhaseSpace satellites

There are a small number of boilerplate lines of code that must go into the constructor of any class in order for the class instance to be integrated into the factory design pattern of Halotools composite models. In this final section of the tutorial, we will look closely at the NFWPhaseSpace constructor to see how the analytical functions and Monte Carlo methods described above get incorporated into the Halotools framework.

As NFWPhaseSpace is primarily a container class for externally-defined behavior, the primary task of its constructor is to call the constructors of its super-classes. There are two super-classes whose behavior is being composed, NFWProfile and MonteCarloGalProf, whose roles we will now describe in turn.

Constructor of the NFWProfile class

The NFWProfile class is a sub-class of AnalyticDensityProf. Calling the constructor of this super-class is the first thing the NFWProfile constructor does. See the Constructor of the AnalyticDensityProf class section for what this accomplishes.

Any NFW profile is specified entirely by mass and concentration, and there exist many different calibrated models for the mapping between halo mass and halo concentration. The conc_mass_model keyword argument allows you to control the behavior of this relation. This keyword allows you to select between different nicknames for conc-mass relations in Halotools, or alternatively you can pass in any callable function that accepts a table keyword argument and returns an array of floats the same length as the Table object bound to the keyword.

Below we show a simple example of how to use a custom concentration-mass relation with NFWPhaseSpace:

def custom_conc_mass(table):
        mass = table['halo_mvir']
        return np.zeros_like(mass) + 5.

nfw_model = NFWPhaseSpace(conc_mass_model=custom_conc_mass)

When using the NFWProfile class in stand-alone fashion, all the analytical functions bound to the instance require that the halo concentration be supplied as an independent argument, and so the behavior inherited by the conc_mass_model keyword argument is irrelevant in such cases. However, when the NFWProfile class is used as a component model in the mock-making framework, the mapping between a halo catalog and the halo concentration must be provided. For such a use-case, the conc_mass_model keyword argument provides the user with the direct_from_halo_catalog model option to simply use the concentration in the halo_table itself.

As described in the Explanation of the prof_param_keys mechanism, one of the boilerplate standardizations of halo profile models is that all sub-classes of the AnalyticalDensityProf class must have a bound method with the same name as every element in the prof_param_keys list of strings bound to the instance. In the case of the NFWProfile, there is only a single profile parameter: conc_NFWmodel. Accordingly, there is a NFWProfile.conc_NFWmodel method; the behavior of this method derives entirely from the model provided by the conc_mass_model keyword argument.

Constructor of the MonteCarloGalProf class

The final super-class constructor called is MonteCarloGalProf.__init__, which performs four functions:

1. A python dictionary called new_haloprop_func_dict is created and bound to the instance.

As described in The new_haloprop_func_dict mechanism, the purpose of new_haloprop_func_dict is to create a new column in the halo_table before mock-population begins, and to automatically guarantee that all galaxies in the galaxy_table created during mock-population will inherit whatever new halo property is created. The keys of new_haloprop_func_dict are the names of the to-be-inherited property, the values are function objects that operate on the original halo_table and return the value of the newly created halo property. If component models require access to a halo property that is not already in the halo_table and this quantity is expensive to calculate, this mechanism can save considerable runtime during mock-population as it can be computed in advance.

In the case of our NFWPhaseSpace model, we calculate the conc_NFWmodel property with new_haloprop_func_dict. The newly created halo_table key will be called conc_NFWmodel, and the value bound to this key will be whatever result is returned by the NFWPhaseSpace.conc_NFWmodel function.

2. A numpy.dtype object called _galprop_dtypes_to_allocate is created and bound to the instance.

As described in The galprop_dtypes_to_allocate mechanism, the purpose of _galprop_dtypes_to_allocate is to inform the HodModelFactory the name and data type of the galaxy attributes that will be created by the component model, so that the appropriate memory can be pre-allocated without any hard-coding in the HodModelFactory. For the case of our NFWPhaseSpace model, we require x, y, z, vx, vy, vz, and we also allocate host_centric_distance as this is an interesting physical characteristic of satellite galaxies upon which other properties defined elsewhere may depend.

3. Build lookup tables for the spatial and velocity profiles using MonteCarloGalProf.setup_prof_lookup_tables.

The purpose of these lookup tables is to improve performance of the Monte Carlo generation of mock galaxy spatial positions and velocities. See Building lookup tables for the spatial and velocity profiles for details.

4. A python list called _mock_generation_calling_sequence is created and bound to the instance.

This list determines which bound methods of NFWPhaseSpace will be called during mock-population, and in which order. The NFWPhaseSpace model only has a single such method, assign_phase_space, which itself simply calls the mc_pos and mc_vel methods in turn. See The mock_generation_calling_sequence mechanism for details.