weird IMF assumed in multidim sampler

[kareemelbadry]

I've been using the multidim sampler to generate a binary population following Moe & DiStefano 2017:

from cosmic.sample.initialbinarytable import InitialBinaryTable
from cosmic.sample.sampler import multidim

InitialBinaries, mass_sin, mass_bin, n_sin, n_bin = InitialBinaryTable.sampler('multidim', final_kstar1=np.linspace(0, 14, 15), final_kstar2=np.linspace(0, 14, 15), rand_seed=0, nproc=1, SF_start=13700.0, SF_duration=0.0, met=0.02, size=10000)

I noticed that this results in a total mass distribution (primaries + secondaries + singles) that is significantly more top-heavy than a standard Kroupa IMF: default_multidim.pdf (this plot compares the generated population (red) to draws from a Kroupa IMF with the same total mass (black). Edit: the x-axis should be labeled log(M/Msun), not M/Msun.)

The reason for this seems to be that the code (lines 780-790 of multidim.py) assumes a rather top-heavy primary mass function, with alpha = 1.6 between 0.5 and 1 Msun (where Kroupa is 2.3) and alpha = 0.8 between 0.08 and 0.5 Msun (where Kroupa is 1.3).

As far as I can tell, this mass function comes from Eq 5 of https://arxiv.org/abs/1708.08248. That paper quotes Moe & DiStefano, but I think Moe & DiStefano don't actually say anything about the IMF or primary mass function (at least in the paper; maybe they did make these assumptions in their IDL code.)

I think it might be more appropriate to set the alphas to 2.3 and 1.3, which would give much better agreement of the total mass distribution with Kroupa: modified_multidim.pdf

If you care mostly about massive stars, this won't directly affect the massive binaries, but it will change the conversion between simulated binary population and total stellar population mass, by close to a factor of 2.

Or is there a known reason for assuming this more top-heavy mass function?

[katiebreivik]

Hey @kareemelbadry! Sorry this has taken me quite a long time to get to checking on.

The code in the multidim is indeed basically a one-to-one copy from IDL to Python as far as I'm aware. That said -- we could add in the option for different slopes. Is this something you're hoping for in the short term?

[kareemelbadry]

Thanks! I changed the slope in my local version, so I don't need changes to the public repo in the short term. I just wanted to point out the issue because a slope of -1.6 is kind of random (and is not used or justified in the Moe & Di Stefano paper), and I imagine most users aren't aware they're using it.