Add analytical expression for average star-forming mass evolved per binary

[avivajpeyi]

We sample from an inverse IMF CDF to compute the average star-forming mass evolved per binary. This can be done analytically and much faster (requiring less code!)

@ilyamandel shared the following

Mtot=(1/0.08)*(1/0.5)*(0.08^1.7-0.01^1.7)/1.7 + (1/0.5)*(0.5^0.7-0.08^0.7)/(0.7) + (-0.5^(-0.3))/(-0.3);
Ntot=(-(0.5)^(-1.3)+Mmax^(-1.3))/(-1.3)+(0.5^(-0.3)-0.08^(-0.3))/(-0.3)*(1/0.5)+(0.08^0.7-0.01^0.7)/0.7*(1/0.5)*(1/0.08);
Nint=(Mmax^(-1.3)-Mmin^(-1.3))/(-1.3);
Mtot*1.5/Nint

(where Mmax, Mmin are the COMPAS run's primary mass limits)

This does not account for the compas_m2_min = 0.1 Msun, and needs to be worked in to get the analytical function that we want.

[ilyamandel]

See https://github.com/TeamCOMPAS/COMPAS/issues/378 ;-) But yes, I'm on the hook for adding M2,min correctly to the above -- will do.

[ilyamandel]

@TomWagg -- it should have been 3 lines (Ntot isn't actually used). :)

Seriously, here's a slightly longer version -- with the min on m2 in place.

%p(M) = alpha * M^-2.3 for M between m3 and m4; 
%p(M) \propto M^-1.3 for M between m2 and m3; 
%p(M) \propto M^-0.3 for M between m1 and m2
m1=0.01; m2=0.08; m3=0.5; m4=200; 
Mmax = 150; Mmin = 5; M2min = 0.1; %thresholds on masses of M1 and M2 (otherwise, q=M2/M1 drawn from U[0,1]
alpha = (-(m4^(-1.3)-m3^(-1.3))/1.3 - (m3^(-0.3)-m2^(-0.3))/(m3*0.3) + (m2^0.7-m1^0.7)/(m2*m3*0.7))^(-1)
%average mass of stars (average mass of all binaries is a factor of 1.5 larger)
Mave = alpha * (-(m4^(-0.3)-m3^(-0.3))/0.3 + (m3^0.7-m2^0.7)/(m3*0.7) + (m2^1.7-m1^1.7)/(m2*m3*1.7))
%fraction of binaries that COMPAS simulates
fint = -alpha/1.3 * (Mmax^(-1.3) - Mmin^(-1.3)) + alpha* M2min/2.3 * (Mmax^(-2.3) - Mmin^(-2.3)) 
%mass represented by each binary simulated by COMPAS
Mrep = (1/fint) * Mave * 1.5

[avivajpeyi]

My apologies @ilyamandel -- I dont quite follow this 😅 In the Mave LHS, where does the 1.7 come from?

[avivajpeyi]

from compas_python_utils.cosmic_integration.ClassCOMPAS import COMPASData

def ilya_find_star_forming_mass_per_binary(Mmax, Mmin, M2min, m1=0.01, m2=0.08, m3=0.5, m4=200, a1=1.3, a2=0.3, a3=0.7):
    alpha = (-(m4**(-1.3)-m3**(-1.3))/1.3 - (m3**(-0.3)-m2**(-0.3))/(m3*0.3) + (m2**0.7-m1**0.7)/(m2*m3*0.7))**(-1)
    # average mass of stars (average mass of all binaries is a factor of 1.5 larger)
    Mave = alpha * (-(m4**(-0.3)-m3**(-0.3))/0.3 + (m3**0.7-m2**0.7)/(m3*0.7) + (m2**1.7-m1**1.7)/(m2*m3*1.7))
    # fraction of binaries that COMPAS simulates
    fint = -alpha/1.3 * (Mmax**(-1.3) - Mmin**(-1.3)) + alpha* M2min/2.3 * (Mmax**(-2.3) - Mmin**(-2.3))
    # mass represented by each binary simulated by COMPAS
    Mrep = (1/fint) * Mave * 1.5
    return Mrep

if __name__ == '__main__':
    import astropy.units as u
    Mmax = 150
    Mmin = 5
    M2min = 0.1
    data = COMPASData()
    data.Mupper = Mmax * u.M_sun
    data.Mlower = Mmin* u.M_sun
    data.m2_min = M2min* u.M_sun
    data.binaryFraction = 1
    data.find_star_forming_mass_per_binary_sampling()

    print(f"Ilya: {ilya_find_star_forming_mass_per_binary(Mmax=Mmax, Mmin=Mmin, M2min=M2min)}")
    print(f"COMPASData.find_star_forming_mass_per_binary_sampling: {data.mass_evolved_per_binary}")

Ilya: 79.11015055714887
COMPASData.find_star_forming_mass_per_binary_sampling: 79.01797410263515 solMass

@ilyamandel -- is it possible to incorporate an fraction_binary for m1 and m2 drawn from the IMF?

[ilyamandel]

My apologies @ilyamandel -- I dont quite follow this 😅 In the Mave LHS, where does the 1.7 come from?

1.7 comes from integrating M p(M) dM when p(M) \propto M^{-0.3}.

[ilyamandel]

is it possible to incorporate an fraction_binary

I generally assume that single stars are just wide binaries, and COMPAS therefore already accounts for them. :)

But if you want to account for fbin \neq 1, you can just replace Mrep = (1/fint) * Mave * 1.5 with Mrep = (1/fint) * Mave * (1.5 + (1-fbin)/fbin)

[avivajpeyi]

I generally assume that single stars are just wide binaries, and COMPAS therefore already accounts for them. :)

ah! it seems like the ComicIntegrator code requires us to input the fraction of binaries... Is that not required?

It does change Mrep quite a bit (and it seems like fbin is set to 0.7 by default)

https://github.com/TeamCOMPAS/COMPAS/blob/bc97c0c683977d9b29a4bc3df834b8871f9c30d6/compas_python_utils/cosmic_integration/FastCosmicIntegration.py#L779