[general issue] `mass_binaries` != sum(mass_1 + mass_2)?

[TomWagg]

I imagine I'm missing something here but I wanted to ask about this case in which I thought the total mass of the sampled binaries should be the same as mass_binaries. Here's the call:

initial_binaries, _, mass_binaries, _, _ = InitialBinaryTable.sampler('independent',
                                                            np.linspace(0, 15, 16), np.linspace(0, 15, 16),
                                                            binfrac_model=1.0,
                                                            primary_model='kroupa01', ecc_model='sana12',
                                                            porb_model='sana12', qmin=-1,
                                                            SF_start=particle["t_form"],
                                                            SF_duration=0.0, met=particle["Z"],
                                                            size=int(np.ceil(mass)))

I figured that since my kstars have everything and the binfrac=1.0 then it should be the case that approximately

    mass_binaries / (initial_binaries["mass_1"] + initial_binaries["mass_2"]).sum() == 1.0

since nothing would get masked. But I'm finding that ratio fluctuates around 1.08 (from ~1.075-1.085) (independent of metallicity/formation time/size)

Am I missing something that makes it so this isn't 1.0? Thanks!

[TomWagg]

Perhaps another clue - I checked the outputted n_binaries vs. the length of initial_binaries and it's always 50% higher.

[TomWagg]

Okay I believe I understand where the difference is coming from now - it's the sampled secondary stars!

I think this is because the minimum stellar mass that COSMIC/BSE can handle and the absolute minimum stellar mass are not the same. Although the primary stars are sampled with a minimum mass of 0.08, the secondaries are sampled with q and so can have a mass as low as 0.008 (side note: is that big enough for a star?) and any binary with a secondary mass < 0.08 gets thrown out.

Apparently on average this comes to ~1/3 of the sampled population and ~7% of the mass...which I guess would be obvious if I took the time to integrate the IMF 🤷🏻

[TomWagg]

Okay I was curious so I did the integration:

@np.vectorize
def IMF(m, m1=0.01, m2=0.08, m3=0.5, m4=200.0, a12=0.3, a23=1.3, a34=2.3):
    """ 
        Calculate the fraction of stellar mass between m and m + dm for a three part broken power law.
        Default values follow Kroupa (2001)
            zeta(m) ~ m^(-a_ij)
        
        Args:
            m       --> [float, list of floats] mass or masses at which to evaluate
            mi      --> [float]                 masses at which to transition the slope
            aij     --> [float]                 slope of the IMF between mi and mj
            
        Returns:
            zeta(m) --> [float, list of floats] value or values of the IMF at m
    """
    # calculate normalisation constants that ensure the IMF is continuous
    b1 = 1 / ( 
                (m2**(1 - a12) - m1**(1 - a12)) / (1 - a12) \
                + m2**(-(a12 - a23)) * (m3**(1 - a23) - m2**(1 - a23)) / (1 - a23) \
                + m2**(-(a12 - a23)) * m3**(-(a23 - a34)) * (m4**(1 - a34) - m3**(1 - a34)) / (1 - a34)
                )
    b2 = b1 * m2**(-(a12 - a23))
    b3 = b2 * m3**(-(a23 - a34))

    # evaluate IMF
    if m < m1:
        return 0
    elif m < m2:
        return b1 * m**(-a12)
    elif m < m3:
        return b2 * m**(-a23)
    elif m < m4:
        return b3 * m**(-a34)
    else:
        return 0

m_range = np.linspace(0.01, 0.08, 100000)
imf_vals = IMF(m_range, m1=0.01, m4=150)

print(np.trapz(imf_vals, m_range))
print(np.trapz(imf_vals * m_range, m_range))

This gives

0.3714647956325013
0.015493749559631911

Which checks out! About a 1/3 of the IMF is < 0.08 and though only ~1.5% of the mass is there, I guess the primary companions bring it up to 7% of the binary mass?