Independent distributions¶
This guide will show you how to sample an initial binary population using the independent sampler in COSMIC.
First import the InitialBinaryTable class and the independent sampler
In [1]: from cosmic.sample.initialbinarytable import InitialBinaryTable
In [2]: from cosmic.sample.sampler import independent
The independent sampler contains multiple models for each binary parameter. You can access the available models using the independent sampler help call:
In [3]: help(independent.get_independent_sampler)
Help on function get_independent_sampler in module cosmic.sample.sampler.independent:
get_independent_sampler(final_kstar1, final_kstar2, primary_model, ecc_model, porb_model, SF_start, SF_duration, binfrac_model, met, size=None, total_mass=inf, sampling_target='size', trim_extra_samples=False, q_power_law=0, **kwargs)
Generates an initial binary sample according to user specified models
Parameters
----------
final_kstar1 : `int or list`
Int or list of final kstar1
final_kstar2 : `int or list`
Int or list of final kstar2
primary_model : `str`
Model to sample primary mass; choices include: kroupa93, kroupa01, salpeter55, custom
if 'custom' is selected, must also pass arguemts:
alphas : `array`
list of power law indicies
mcuts : `array`
breaks in the power laws.
e.g. alphas=[-1.3,-2.3,-2.3],mcuts=[0.08,0.5,1.0,150.] reproduces standard Kroupa2001 IMF
ecc_model : `str`
Model to sample eccentricity; choices include: thermal, uniform, sana12
porb_model : `str` or `dict`
Model to sample orbital period; choices include: log_uniform, sana12, raghavan10, moe19
or a custom power law distribution defined with a dictionary with keys "min", "max", and "slope"
(e.g. {"min": 0.15, "max": 0.55, "slope": -0.55}) would reproduce the Sana+2012 distribution
qmin : `float`
kwarg which sets the minimum mass ratio for sampling the secondary
where the mass ratio distribution is flat in q
if q > 0, qmin sets the minimum mass ratio
q = -1, this limits the minimum mass ratio to be set such that
the pre-MS lifetime of the secondary is not longer than the full
lifetime of the primary if it were to evolve as a single star
m_max : `float`
kwarg which sets the maximum primary and secondary mass for sampling
NOTE: this value changes the range of the IMF and should *not* be used
as a means of selecting certain kstar types!
m1_min : `float`
kwarg which sets the minimum primary mass for sampling
NOTE: this value changes the range of the IMF and should *not* be used
as a means of selecting certain kstar types!
m2_min : `float`
kwarg which sets the minimum secondary mass for sampling
the secondary as uniform in mass_2 between m2_min and mass_1
msort : `float`
Stars with M>msort can have different pairing and sampling of companions
qmin_msort : `float`
Same as qmin for M>msort
m2_min_msort : `float`
Same as m2_min for M>msort
SF_start : `float`
Time in the past when star formation initiates in Myr
SF_duration : `float`
Duration of constant star formation beginning from SF_Start in Myr
binfrac_model : `str or float`
Model for binary fraction; choices include: vanHaaften, offner22, or a fraction where 1.0 is 100% binaries
binfrac_model_msort : `str or float`
Same as binfrac_model for M>msort
met : `float`
Sets the metallicity of the binary population where solar metallicity is zsun
size : `int`
Size of the population to sample
total_mass : `float`
Total mass to use as a target for sampling
sampling_target : `str`
Which type of target to use for sampling (either "size" or "total_mass"), by default "size".
Note that `total_mass` must not be None when `sampling_target=="total_mass"`.
trim_extra_samples : `str`
Whether to trim the sampled population so that the total mass sampled is as close as possible to
`total_mass`. Ignored when `sampling_target==size`.
Note that given the discrete mass of stars, this could mean your sample is off by 300
solar masses in the worst case scenario (of a 150+150 binary being sampled). In reality the majority
of cases track the target total mass to within a solar mass.
zsun : `float`
optional kwarg for setting effective radii, default is 0.02
q_power_law : `float`
Exponent for the mass ratio distribution power law, default is 0 (flat in q). Note that
q_power_law cannot be exactly -1, as this would result in a divergent distribution.
Returns
-------
InitialBinaryTable : `pandas.DataFrame`
DataFrame in the format of the InitialBinaryTable
mass_singles : `float`
Total mass in single stars needed to generate population
mass_binaries : `float`
Total mass in binaries needed to generate population
n_singles : `int`
Number of single stars needed to generate a population
n_binaries : `int`
Number of binaries needed to generate a population
Targetting specific final kstar types¶
The final_kstar1 and final_kstar2 parameters are lists that contain the kstar types
that you would like the final population to contain.
The final kstar is the final state of the binary system we are interested in and is based on the BSE kstar naming convention, see Evolutionary State of the Star for more information.
Thus, if you want to generate a population containing double white dwarfs with CO and ONe WD primaries and He-WD secondaries, the final kstar inputs would be:
In [4]: final_kstar1 = [11, 12]
In [5]: final_kstar2 = [10]
Since we are interested in binaries, we only retain the binary systems that are likely to produce the user specified final kstar types. However, we also keep track of the total mass of the single and binary stars as well as the number of binary and single stars so that we can scale our results to larger populations. If you don’t want to filter the binaries, you can supply final kstars as
In [6]: final_kstars = np.linspace(0, 14, 15)
In [7]: InitialBinaries, mass_singles, mass_binaries, n_singles, n_binaries = InitialBinaryTable.sampler('independent', final_kstars, final_kstars, binfrac_model=0.5, primary_model='kroupa01', ecc_model='sana12', porb_model='sana12', qmin=-1, m2_min=0.08, SF_start=13700.0, SF_duration=0.0, met=0.02, size=10000)
Additionally if you are interested in single stars then you can specify keep_singles=True.
Understanding parameter sampling models¶
Similar to the help for the sampler, the different models that can be used for each parameter to be sampled can be accessed by the help function for the argument. The syntax for each parameter sample is always: sample_`parameter`. See the example for the star formation history (SFH) below:
In [8]: help(independent.Sample.sample_SFH)
Help on function sample_SFH in module cosmic.sample.sampler.independent:
sample_SFH(self, SF_start=13700.0, SF_duration=0.0, met=0.02, size=None)
Sample an evolution time for each binary based on a user-specified
time at the start of star formation and the duration of star formation.
The default is a burst of star formation 13,700 Myr in the past.
Parameters
----------
SF_start : float
Time in the past when star formation initiates in Myr
SF_duration : float
Duration of constant star formation beginning from SF_Start in Myr
met : float
metallicity of the population [Z_sun = 0.02]
Default: 0.02
size : int, optional
number of evolution times to sample
NOTE: this is set in cosmic-pop call as Nstep
Returns
-------
tphys : array
array of evolution times of size=size
metallicity : array
array of metallicities
Stopping conditions for sampling¶
In addition to telling COSMIC how to sample, you also need to tell it how much to sample. You can
do this in one of two ways: (1) by specifying the number of binaries to sample, or (2) by specifying the total mass of singles and binaries to sample. Let’s look at both of these in more detail.
Number of binaries¶
Using the final kstar inputs we mentioned above, the initial binary population can be sampled based on the desired number of binaries (10000 in this case) as follows
In [9]: InitialBinaries, mass_singles, mass_binaries, n_singles, n_binaries = InitialBinaryTable.sampler('independent', final_kstar1, final_kstar2, binfrac_model=0.5, primary_model='kroupa01', ecc_model='sana12', porb_model='sana12', qmin=-1, m2_min=0.08, SF_start=13700.0, SF_duration=0.0, met=0.02, size=10000)
In [10]: print(InitialBinaries)
kstar_1 kstar_2 mass_1 mass_2 ... bhspin_1 bhspin_2 tphys binfrac
0 1.0 1.0 1.715988 1.490909 ... 0.0 0.0 0.0 0.5
1 1.0 1.0 0.863649 0.803246 ... 0.0 0.0 0.0 0.5
2 1.0 0.0 0.922107 0.699386 ... 0.0 0.0 0.0 0.5
3 1.0 1.0 1.549598 0.993874 ... 0.0 0.0 0.0 0.5
4 1.0 1.0 1.287078 0.772630 ... 0.0 0.0 0.0 0.5
... ... ... ... ... ... ... ... ... ...
10039 1.0 1.0 2.560260 1.892963 ... 0.0 0.0 0.0 0.5
10040 1.0 1.0 0.983484 0.876515 ... 0.0 0.0 0.0 0.5
10041 1.0 0.0 1.138245 0.550008 ... 0.0 0.0 0.0 0.5
10042 1.0 1.0 1.269011 1.257934 ... 0.0 0.0 0.0 0.5
10043 1.0 1.0 2.698330 1.800682 ... 0.0 0.0 0.0 0.5
[10044 rows x 38 columns]
Note
The length of the initial binary data set, InitialBinaries, does not always match
the size parameter provided to sampler().
This is because of the various cuts that the sampler makes to the population (e.g. the binary fraction,
which is either a fraction between 0 and 1 or mass dependent following the
prescription in van Haaften+2013.) specified by the user.
Total mass sampled¶
Alternatively, we could do the same thing but now instead set our sampling_target to be the total mass and aim for 15000 solar masses. This is done by setting sampling_target="total_mass" and total_mass=15000.
In [11]: InitialBinaries, mass_singles, mass_binaries, n_singles, n_binaries = InitialBinaryTable.sampler('independent', final_kstar1, final_kstar2, binfrac_model=0.5, primary_model='kroupa01', ecc_model='sana12', porb_model='sana12', qmin=-1, m2_min=0.08, SF_start=13700.0, SF_duration=0.0, met=0.02, sampling_target="total_mass", total_mass=15000)
In [12]: print(InitialBinaries)
kstar_1 kstar_2 mass_1 mass_2 ... bhspin_1 bhspin_2 tphys binfrac
0 1.0 1.0 1.035583 0.835998 ... 0.0 0.0 0.0 0.5
1 1.0 1.0 0.951510 0.924857 ... 0.0 0.0 0.0 0.5
2 1.0 1.0 1.351939 0.727992 ... 0.0 0.0 0.0 0.5
3 1.0 0.0 3.239166 0.629694 ... 0.0 0.0 0.0 0.5
4 1.0 1.0 2.002929 0.921669 ... 0.0 0.0 0.0 0.5
... ... ... ... ... ... ... ... ... ...
1332 1.0 1.0 6.372072 0.769274 ... 0.0 0.0 0.0 0.5
1333 1.0 1.0 2.239466 0.884469 ... 0.0 0.0 0.0 0.5
1334 1.0 1.0 2.466188 2.383919 ... 0.0 0.0 0.0 0.5
1335 1.0 0.0 1.632684 0.550748 ... 0.0 0.0 0.0 0.5
1336 1.0 1.0 1.214975 1.155442 ... 0.0 0.0 0.0 0.5
[1337 rows x 38 columns]
And we can check what the total sampled mass was by looking at the sum of the mass_singles and mass_binaries variables
In [13]: print(mass_singles + mass_binaries)
22816.602858103794
Tip
If you’d like to avoid your sample overshooting your desired total_mass and instead get as close to this value as possible,
you can set trim_extra_samples=True. This will trim the sample to get a total mass as close as possible to your target.
In many cases, this will be within a solar mass, but could be as large as twice the maximum stellar mass (for the very rare case that
the final binary drawn is the most massive primary with an equal mass ratio).
Mass dependent binary fractions and mass pairings¶
If you want to have separate binary fractions and mass pairings for low and high mass stars, you can by supplying the msort kwarg to the sampler. This sets the mass above which an alternative mass pairing (specified by kwargs qmin_msort and m2_min_msort) and binary fraction model (specified by kwarg binfrac_model_msort) are used. This is handy if you want, for example, a higher binary fraction and more equal mass pairings for high mass stars.
Below we show the effect of different assumptions for the independent initial sampler. The standard assumptions are shown in purple, the assumptions of Sana et al. 2012 are shown in orange, and the assumptions of Moe et al. 2019 are shown in green.
(Source code, png, hires.png, pdf)
