mass function normalization

[damonge]

The mass function should be normalized so that one recovers the mean matter density when integrating over masses. I believe this is not the case in the current implementation.

Running the code below should reproduce the issue.

import numpy as np
import pyccl as ccl
from scipy.integrate import quad

cpar={'h':0.6774,
      'oc':0.1188,
      'ob':0.02230,
      'A_s':2.142E-9,
      'n_s': 0.9667,
      'delta':200.}
cpar['Om']=(cpar['oc']+cpar['ob'])/cpar['h']**2
cpar['Oc']= cpar['oc']/cpar['h']**2
cpar['Ob']= cpar['ob']/cpar['h']**2

cosmo=ccl.Cosmology(Omega_c=cpar['Oc'],Omega_b=cpar['Ob'],h=cpar['h'],A_s=cpar['A_s'],
                    n_s=cpar['n_s'])

def norm() :
    def integ(lm) :
        m=10.**lm
        nm=m*ccl.massfunc(cosmo,m,1.,odelta=cpar['delta'])/(2.744948E11*cpar['Om']*cpar['h']**2)
        return nm
    return quad(integ,4,17)[0]

print "This should be 1 : %.3lF"%(norm())

[tmcclintock]

Can you post the output as well?

[alexander-mead]

What are your integration limits? It looks like 10^4 to 10^17 solar masses. I guess there is this standard problem that the mass function tends to slowly diverge at low masses and so getting correct limit is quite hard in practice, but it’s usually just a numerical problem…

For example, with the Sheth & Torman function (which I know a bit better) if you integrate between 100 solar mass haloes and 10^17 solar mass haloes you only recover ~half of the mass at z=0. This is telling you that, if you take the mass function literally, you should expect ~half of the mass in the universe to reside in haloes smaller than 100 Solar masses!

Whether the mass function is correctly normalised (if you integrate correctly from 0 to infinity in mass) depends on which version of the Tinker function you use. I think the fitting parameters from those in the 2010 paper are constrained to be correctly normalised. Also the Tinker bias fit is constrained to be correctly normalised such that the mean bias of matter = 1. Sheth & Tormen is correctly normalised, as is the peak-background split derived bias from Sheth, Mo & Tormen…

[damonge]

Right, this shouldn't be a problem with the normalization, since our benchmarks at M>1E8 pass, and I agree it must have to do with whatever happens on low masses.

For most purposes this may be fine, but the halo model involves these types of integrals, so it may be that it matters for #339 .

@tmcclintock : the output is "This should be 1 : 0.726"

[damonge]

I think the integrals I'm talking about are probably much better behaved when expressed as a function of nu, so that may be worth bearing in mind.

[m-aguena]

If I may contribute to the conversation, I agree with @damonge. This integral usually behave better as a function of nu (or sigma). One of the main issues of doing it on mass is the conversion of mass->sigma, that is usually difficult to have at low masses.

How is the sigma(M) being computed? is it an interpolation of some kind?

[alexander-mead]

I just checked my own implementation of Tinker, which has no splines, and integrating between 10^3 and 10^17 Solar masses gets 0.555, rather than 1, for a vanilla LCDM cosmology at z=0. I think this is correct given what I said above. I think you need to integrate down to ~Jupiter mass before this starts to converge.

[damonge]

OK, great. Let's close this for the time being then. I guess @alexander-mead and @villarrealas are being smart at this with their halo model work!

[alexander-mead]

The Sheth & Tormen function diverges like nu^-0.6 at low nu, so it is still quite slow to converge when doing the integral in terms of nu.

What I have done in the past to deal with this type of issue is to basically compute the integral down to a lower limit (say 10^2 solar masses) and then check how far away from unity I am. I then alter the mass function to be whatever one I was using between by integration limits, and then stick the rest of the mass in as a delta function at the lower integration limit. This works fine for computing one- and two-halo terms in halo-model approaches: these low mass haloes don't contribute to the one-halo term anyway (so, actually you don't need to worry about this at all for a one-halo term) and for the two-halo term you just need to ensure that the normalisation is correct so that linear theory is recovered at large scales.

What you definitely should not do is to multiply the mass function by an overall amplitude to correct for the fact that you are missing mass, as this amounts to saying that there are no low-mass haloes and ~twice as many high-mass haloes. This is crazy, but you see it done in papers sometimes.

[damonge]

This makes a lot of sense, actually.

[damonge]

Actually, one question @alexander-mead : does this trick work also to ensure that the bias of dark matter is 1? (assuming you use the trick to ensure that you recover rho_M)

[alexander-mead]

I think so, in fact, the way I apply this is at the level of assuming a delta function in 'mass function x bias', rather than in 'mass function'. It should not matter though because for most (all?) bias functions the low-mass haloes have a constant anti-bias, so making the delta function in either 'mass function x bias' or 'mass function' should not matter.