Generate a population the COSMIC way¶

Beyond providing a simple interface and several updates to BSE (Hurley+2002), COSMIC is designed to adapt the number of systems simulated to each binary evolution model a user selects. This is done by iteratively simulating binaries initialized with ZAMS binary parameters and a star formation history until the distributions of binary parameters specified by the user converge. This process is carried out by cosmic-pop. Once this population is simulated, an astrophysical population can be sampled from the output of cosmic-pop.

cosmic-pop¶

The fixed population is the output of cosmic-pop and is designed to fully describe the distribution of binary parameters that results from a user specified star formation history (SFH) and binary evolution model (BSEDict, specified in an inifile). The fixed population is only simulated once and only contains information about the intrinsic properties of the binary (e.g. masses, semimajor axes, metallicities, etc.) Information about the location and orientation of each binary is not contained in the fixed population. The arguments necessary to run the cosmic-pop executable can be found using the help command:

cosmic-pop -h

usage: cosmic-pop [-h] [--inifile FILE] --final-kstar1 FINAL_KSTAR1
                  [FINAL_KSTAR1 ...] --final-kstar2 FINAL_KSTAR2
                  [FINAL_KSTAR2 ...] [--Niter NITER] [--Nstep NSTEP]
                  [-n NPROC] [--binary_state BINARY_STATE [BINARY_STATE ...]]
                  [--sampling_method SAMPLING_METHOD]
                  [--primary_model PRIMARY_MODEL]
                  [--binfrac_model BINFRAC_MODEL] [--ecc_model ECC_MODEL]
                  [--porb_model PORB_MODEL] [--SF_start SF_START]
                  [--SF_duration SF_DURATION] [--metallicity METALLICITY]
                  [--convergence_params CONVERGENCE_PARAMS [CONVERGENCE_PARAMS ...]]
                  [--convergence_limits CONVERGENCE_LIMITS]
                  [--convergence_filter CONVERGENCE_FILTER] [--match MATCH]
                  [--bcm_bpp_initcond_filter [BCM_BPP_INITCOND_FILTER]]
                  [--seed SEED] [--verbose]

optional arguments:
  -h, --help            show this help message and exit
  --inifile FILE        Name of ini file of params
  --final-kstar1 FINAL_KSTAR1 [FINAL_KSTAR1 ...]
                        Specify the final condition of kstar1 , you want
                        systems to end at for your samples
  --final-kstar2 FINAL_KSTAR2 [FINAL_KSTAR2 ...]
                        Specify the final condition of kstar2, you want
                        systems to end at for your samples
  --Niter NITER         Number of iterations of binaries to try, will check
                        ever Nstep for convergence
  --Nstep NSTEP         Number of binaries to try before checking for
                        convergence, it will check ever Nstep binaries until
                        it reach Niter binaries
  -n NPROC, --nproc NPROC
                        number of processors
  --binary_state BINARY_STATE [BINARY_STATE ...]
  --sampling_method SAMPLING_METHOD
  --primary_model PRIMARY_MODEL
                        Chooses the initial primary mass function from:
                        salpeter55, kroupa93, kroupa01
  --binfrac_model BINFRAC_MODEL
                        Chooses the binary fraction model from: a float
                        between [0,1] and vanHaaften
  --ecc_model ECC_MODEL
                        Chooses the initial eccentricity distribution model
                        from: thermal, uniform, and sana12
  --porb_model PORB_MODEL
                        Chooses the initial orbital period distribution model
                        from: log_uniform and sana12
  --SF_start SF_START   Sets the time in the past when star formation
                        initiates in Myr
  --SF_duration SF_DURATION
                        Sets the duration of constant star formation in Myr
  --metallicity METALLICITY
  --convergence_params CONVERGENCE_PARAMS [CONVERGENCE_PARAMS ...]
  --convergence_limits CONVERGENCE_LIMITS
  --convergence_filter CONVERGENCE_FILTER
  --match MATCH
  --bcm_bpp_initcond_filter [BCM_BPP_INITCOND_FILTER]
  --seed SEED
  --verbose             Run in Verbose Mode

Inputs¶

Params.ini¶

PLEASE SEE Configuration files for COSMIC for detailed information about the inifile and how it is constructed

The inifile contains five subsections: filters, convergence, rand_seed, and bse.

The filters subsection allows you to specify how you would like to filter the binary population. See the inifile for a description of each flag.

The convergence subsection allows you to specify a particular region of parameter space where you would like the convergence of each binary parameter distribution to be tested. The only implemented convergence filter is for LISA binaries, where the convergence is computed only for binaries with orbital period less than 5000 s. This allows for low probability, but high signal to noise binaries with very short orbital periods to be fully accounted for.

The rand_seed subsections allows you to initialize the population with a random seed for reproduceability. Note that for each simulated binary, COSMIC also returns the initial conditions, including a random seed that will reproduce that exact binary. The random seed produced in the rand_seed subsection allows full populations to be reproduced. This is particularly useful when comparing two popuations, e.g. binary black holes and binary neutron stars, which should be simulated separately, but using the same rand_seed value.

The sampling subsection allows you to control how you sample and achieve the initial binaries that you evolve.

The bse subsection is where all the BSE flags are specified.

Sample command line¶

Below is an example command line call for cosmic-pop:

cosmic-pop --final-kstar1 13 14 --final-kstar2 13 14 --inifile Params.ini --Nstep 1000 --Niter 1000000000 -n 2

A breakdown of each argument follows:

--final-kstar1 13 14 --final-kstar2 13 14 : this tells COSMIC to keep all systems where the primary star is a BH or NS and the secondary star is also a BH or NS.
--inifile Params.ini : this tells COSMIC where to look for the BSE flags that set the binary evolution model; in this case, the inifile is assumed to be in the same directory that the command line call is made
--Nstep 5000 --Niter 10000000 -n 1 : this tells COSMIC to use 1 processor to evolve a maximum of 1e7 systems and check in every 5000 systems to track how the shape of the distributions of the parameters specified in convergence-params change

Stopping conditions¶

There are two stopping conditions for cosmic-pop:

The shape of the parameter distributions for the parameters specified in convergence_params section of the inifile file converge to a shape, regardless of adding new binaries to the evolved population. This is quantified by the match criteria detailed in Breivik+2019 in prep
The number of binaries sampled exceeds --Niter.

Output of cosmic-pop¶

PLEASE SEE Describing the output of COSMIC/BSE: Columns names/Values/Units for more information about the output data frames including what each column means and the units.

The output of cosmic-pop is the fixed population, an hdf5 file with a naming scheme that tells you the Galactic component and final kstars of the population; the data file created by the cosmic-pop call above is: dat_DeltaBurst_13_14_13_14.h5.

The fixed population contains several pandas DataFrames accessed by the following keys:

conv : The converged population whose parameters are sepcified by the convergence subsection of the inifile
bpp : The evolutionary history of binaries which satisfy the user-specified final kstars and filter in the convergence subsection
bcm : The final state of binaries in the bcm array which satisfy the user-specified final kstars and filter in the convergence subsection
kick_info : The magnitude and direction of natal kicks, three dimensional systemic velocity changes, total tilt of orbital plane, and azimuthal angle of orbital angular momentum axis with respect to spins
initCond : The initial conditions for each binary which satisfies the user-specified final kstars and filter in the convergence subsection
idx : An integer that keeps track of the total number of simulated binaries to maintain proper indexing across several runs of cosmic-pop
match : Tracks the convergence where match = Log ₁₀ (1-convergence)
mass_binaries : Tracks the total mass of binaries needed to create the fixed population
mass_singles : Tracks the total mass of single stars needed to create the fixed population; if the binary fraction is 100%, the mass in singles will be zero
mass_stars : Tracks the total mass of all stars, including binaries and singles, needed to create the fixed population
n_binaries : Tracks the total number of binaries needed to create the fixed population
n_singles : Tracks the total number of single stars needed to create the fixed population
n_stars : Tracks the total number of stars, where n_stars = n_singles + 2*n_binaries, needed to create the fixed population

The conv, bpp, bcm, and initCond DataFrames share a common column called bin_num which is used to index the population across the DataFrames.

scaling to an astrophysical population¶

Once the fixed population is simulated, you can scale the simulation to an astrophysical population by resampling the conv DataFrame.

First, we need to load the data which is saved in the same directory where` cosmic-pop is called:

In [1]: import pandas

In [2]: import numpy

In [3]: conv = pandas.read_hdf('fixedpop/dat_DeltaBurst_13_14_13_14.h5', key='conv')

In [4]: total_mass = pandas.read_hdf('fixedpop/dat_DeltaBurst_13_14_13_14.h5', key='mass_stars')

In [5]: N_stars = pandas.read_hdf('fixedpop/dat_DeltaBurst_13_14_13_14.h5', key='n_stars')

Note

The masses and numbers of stars/binaries is saved at each iteration, so you’ll need to take the maximum mass/number:

In [6]: print(N_stars)
           0
  9832994
 19664494
 29497317
 39329473
 49161732
 58992422
 68823776
 78655686
 88189369
 98023296
107858668
117692244
127524814
136830882
146137733
155970005
165803742
175636330
185467780
195300614
205131835
214964373
224798029
234631004
244464116

In [7]: total_mass = max(numpy.array(total_mass))[0]

In [8]: N_stars = max(numpy.array(N_stars))[0]

In [9]: print(total_mass, N_stars)
99434460.72687817 244464116

Since COSMIC tracks both the total number and total mass of stars formed, you can either scale by number or mass.

To scale by number, multiply the number of systems in the conv array by the ratio of the total number of stars in the astrophysical population to the total number of stars used to generate the population:

In [10]: N_astro = 1e10

In [11]: N_13_14_13_14_astro = int(len(conv)*N_astro/N_stars)

In [12]: print(N_13_14_13_14_astro)
64590

Instead, if you want to scale by mass, you can choose between supplying your own mass or using the built in COSMIC models. The process for a user-supplied mass is nearly identical to scaling by number:

In [13]: M_astro = 1e11

In [14]: N_13_14_13_14_astro = int(len(conv)*M_astro/total_mass)

In [15]: print(N_13_14_13_14_astro)
1587980

If you specified a Milky Way star formation history in the inifile for the cosmic-pop call (e.g. ThinDisk, Thick Disk, or Bulge), you can easily scale to a representative Milky Way population with the MC_samp module. This will assume masses according to the McMillan 2011 model:

In [16]: from cosmic import MC_samp

In [17]: N_13_14_13_14_ThinDisk = MC_samp.mass_weighted_number(dat=conv, total_sampled_mass=total_mass, component_mass=MC_samp.select_component_mass('ThinDisk'))

In [18]: N_13_14_13_14_Bulge = MC_samp.mass_weighted_number(dat=conv, total_sampled_mass=total_mass, component_mass=MC_samp.select_component_mass('Bulge'))

In [19]: N_13_14_13_14_ThickDisk = MC_samp.mass_weighted_number(dat=conv, total_sampled_mass=total_mass, component_mass=MC_samp.select_component_mass('ThickDisk'))

It is also possible to convolve a star formation history with a DeltaBurst population. As an example, let’s create a population of NSs and BHs in the Milky Way thin disk. First, use the number of NS/BHs for the thin disk calculated above, and resample the conv array:

In [20]: print('The number of NS/BH systems in the simulated Milky Way thin disk is: {}'.format(N_13_14_13_14_ThinDisk))
The number of NS/BH systems in the simulated Milky Way thin disk is: 686007

In [21]: thin_disk_sample = conv.sample(n=N_13_14_13_14_ThinDisk, replace=True)

In [22]: print(thin_disk_sample)
    bin_num      tphys    mass_1    mass_2  kstar_1  kstar_2        sep       porb  ...      aj_1  aj_2         tms_1     tms_2   massc_1   massc_2     rad_1     rad_2
47   127641  31.182366  1.285621  1.277515     13.0     13.0   0.520399   0.450615  ...  4.012187   0.0  1.000000e+10  1.621963  1.285621  1.277515  0.000014  1.999322
11    16863  38.495229  1.278987  1.277515     13.0     13.0   2.827305   5.285162  ...  4.886382   0.0  1.000000e+10  1.995004  1.278987  1.277515  0.000014  9.657490
1     93432   5.683924  3.526645  3.121816     14.0     14.0  18.668544   3.625578  ...  0.302864   0.0  1.000000e+10  0.788577  3.526645  3.121816  0.000015  0.000013
21   120470  49.436932  1.278788  1.277515     13.0     13.0   1.670417   4.146062  ...  8.962992   0.0  1.000000e+10  2.919944  1.278788  1.277515  0.000014  7.439394
45    91011   8.545978  3.906329  1.277515     14.0     13.0   3.222999   0.294535  ...  3.410703   0.0  1.000000e+10  1.624231  3.906329  1.277515  0.000017  0.000014
..      ...        ...       ...       ...      ...      ...        ...        ...  ...       ...   ...           ...       ...       ...       ...       ...       ...
34    39345   4.116547  6.824194  6.406259     14.0     14.0  75.104585  20.738841  ...  0.080446   0.0  1.000000e+10  0.612832  6.824194  6.406259  0.000029  0.000027
61    10052  38.677202  1.279007  1.277515     13.0     13.0   2.307976   5.038857  ...  6.066015   0.0  1.000000e+10  2.092167  1.279007  1.277515  0.000014  9.234511
47    30683  24.045566  1.285094  1.507461     13.0     13.0   0.514052   0.687210  ...  2.587362   0.0  1.000000e+10  1.327204  1.285094  1.507461  0.000014  2.828796
80    88094  32.361136  1.279025  1.277515     13.0     13.0   2.428715   5.277959  ...  8.076725   0.0  1.000000e+10  1.957982  1.279025  1.277515  0.000014  9.695451
27    54542  35.099795  1.279035  1.277515     13.0     13.0   2.938625   5.360727  ...  6.343329   0.0  1.000000e+10  1.913573  1.279035  1.277515  0.000014  9.859501

[686007 rows x 24 columns]

Now, we can assign a uniform birth time, assuming that the age of the thin disk is 10 Gyr:

In [23]: thin_disk_sample['tbirth'] = np.random.uniform(0, 10000, N_13_14_13_14_ThinDisk)

Since we are interested in NSs/BHs that form before the present, we can filter out anything that has a formation time after 10 Gyr:

In [24]: thin_disk_sample = thin_disk_sample.loc[thin_disk_sample.tbirth + thin_disk_sample.tphys < 10000]

This leaves us with a population of NSs/BHs at formation where the formation time is the sum of the birth time and tphys in the conv array.

Assigning positions for Milky Way populations¶

The MC_samp module also contains methods which assign positions in the Milky Way according to the McMillan 2011 model. As an example, let’s assign positions and orientations for the thin disk population:

In [25]: xGx, yGx, zGx, inc, omega, OMEGA = MC_samp.galactic_positions(gx_component='ThinDisk', model='McMillan', size=len(thin_disk_sample))

In [26]: thin_disk_sample['xGx'] = xGx

In [27]: print(thin_disk_sample)
    bin_num      tphys    mass_1    mass_2  kstar_1  kstar_2        sep       porb  ...         tms_1     tms_2   massc_1   massc_2     rad_1     rad_2       tbirth       xGx
47   127641  31.182366  1.285621  1.277515     13.0     13.0   0.520399   0.450615  ...  1.000000e+10  1.621963  1.285621  1.277515  0.000014  1.999322  4315.117233 -2.162754
11    16863  38.495229  1.278987  1.277515     13.0     13.0   2.827305   5.285162  ...  1.000000e+10  1.995004  1.278987  1.277515  0.000014  9.657490  6828.844770 -0.067389
1     93432   5.683924  3.526645  3.121816     14.0     14.0  18.668544   3.625578  ...  1.000000e+10  0.788577  3.526645  3.121816  0.000015  0.000013  8049.138617  2.385584
21   120470  49.436932  1.278788  1.277515     13.0     13.0   1.670417   4.146062  ...  1.000000e+10  2.919944  1.278788  1.277515  0.000014  7.439394  6065.981724 -1.003903
45    91011   8.545978  3.906329  1.277515     14.0     13.0   3.222999   0.294535  ...  1.000000e+10  1.624231  3.906329  1.277515  0.000017  0.000014    33.667666 -6.402419
..      ...        ...       ...       ...      ...      ...        ...        ...  ...           ...       ...       ...       ...       ...       ...          ...       ...
34    39345   4.116547  6.824194  6.406259     14.0     14.0  75.104585  20.738841  ...  1.000000e+10  0.612832  6.824194  6.406259  0.000029  0.000027  6307.201242 -2.232109
61    10052  38.677202  1.279007  1.277515     13.0     13.0   2.307976   5.038857  ...  1.000000e+10  2.092167  1.279007  1.277515  0.000014  9.234511  5702.411183 -2.194809
47    30683  24.045566  1.285094  1.507461     13.0     13.0   0.514052   0.687210  ...  1.000000e+10  1.327204  1.285094  1.507461  0.000014  2.828796  4248.894961 -3.613706
80    88094  32.361136  1.279025  1.277515     13.0     13.0   2.428715   5.277959  ...  1.000000e+10  1.957982  1.279025  1.277515  0.000014  9.695451  1635.290559 -2.225438
27    54542  35.099795  1.279035  1.277515     13.0     13.0   2.938625   5.360727  ...  1.000000e+10  1.913573  1.279035  1.277515  0.000014  9.859501  5958.757172 -0.494579

[684043 rows x 26 columns]

Now we have a representative Milky Way population of NS/BH binaries in the thin disk! TADA!!