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]
                  [--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]
                  [--pop_select POP_SELECT] [--match MATCH]
                  [--apply_convergence_limits [APPLY_CONVERGENCE_LIMITS]]
                  [--seed SEED] [--verbose] [-n NPROC | --mpi]

options:
  -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
  --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 ...]
                        specifies the list of parameters for which you would
                        like to track the distribution shapes for convergence
  --convergence_limits CONVERGENCE_LIMITS
                        dictionary that can contain limits for convergence
                        params
  --pop_select POP_SELECT
                        Used in combination with the specified final_kstar1
                        and final_kstar2 values to select the subpopulation of
                        interest from the evolved population
  --match MATCH         provides the tolerance for the convergence calculation
  --apply_convergence_limits [APPLY_CONVERGENCE_LIMITS]
                        filters the evolved binary population to contain only
                        the binaries that satsify the convergence limits
  --seed SEED
  --verbose             Run in Verbose Mode
  -n NPROC, --nproc NPROC
                        number of processors
  --mpi                 Run with MPI.

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 with multiprocessing (i.e. shared memory on a single computer):

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

Below is an example command line call for cosmic-pop with MPI (i.e. across many computers):

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

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 2 : this tells COSMIC to use 2 processors (or in the MPI case, MPI tasks) 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+2020
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

Now we can generate the astrophysical population:

In [16]: pop_astro = conv.sample(N_13_14_13_14_astro, replace=True)

If you specified a star formation history in the inifile for the cosmic-pop call (i.e. by setting SF_start and SF_duration to follow a star formation history of your choice), the data in the resampled population will also be consistent with that star formation history. If you assigned the same SF_start to all binaries in the cosmic-pop call, you can assign birth times in post processing. As an example, we can assign a uniform birth time, assuming that the age of the population is 10 Gyr:

In [17]: pop_astro['tbirth'] = np.random.uniform(0, 10000, N_13_14_13_14_astro)

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 [18]: pop_astro = pop_astro.loc[pop_astro.tbirth + pop_astro.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.