Evolving a single binary¶

Below is the process to initialize and evolve a binary that could have formed a GW150914-like binary. First, import the modules in COSMIC that initialize and evolve the binary.

In [1]: from cosmic.sample.initialbinarytable import InitialBinaryTable

In [2]: from cosmic.evolve import Evolve

To initialize a single binary, populate the InitialBinaries method in the InitialBinaryTable class. Each initialized binary requires the following parameters:

m1 : ZAMS mass of the primary star in \(M_{\odot}\)
m2 : ZAMS mass of the secondary star in \(M_{\odot}\)
porb : initial orbital period in days
ecc : initial eccentricity
tphysf : total evolution time of the binary in Myr
kstar1 : initial primary stellar type, following the BSE convention
kstar2 : initial secondary stellar type, following the BSE convention
metallicity : metallicity of the population (e.g. \(Z_{\odot}=0.014\))

In [3]: single_binary = InitialBinaryTable.InitialBinaries(m1=85.543645, m2=84.99784, porb=446.795757,
   ...:                                                    ecc=0.448872, tphysf=13700.0,
   ...:                                                    kstar1=1, kstar2=1, metallicity=0.002)
   ...: 

In [4]: print(single_binary)
   kstar_1  kstar_2     mass_1    mass_2  ...  bhspin_1  bhspin_2  tphys  binfrac
0      1.0      1.0  85.543645  84.99784  ...       0.0       0.0    0.0      1.0

[1 rows x 38 columns]

The flags for the various binary evolution prescriptions used in BSE also need to be set. Each flag is saved in the BSEDict dictionary. Note that the BSEDict only needs to be specified the first time a binary is evolved with COSMIC or if you need to change the binary evolution prescriptions.

If you are unfamiliar with these prescriptions, it is highly advised to either run the defaults from the COSMIC install which are consistent with Breivik+2020

In [5]: BSEDict = {'xi': 1.0, 'bhflag': 1, 'neta': 0.5, 'windflag': 3, 'wdflag': 1, 'alpha1': 1.0, 'pts1': 0.001, 'pts3': 0.02, 'pts2': 0.01, 'epsnov': 0.001, 'hewind': 0.5, 'ck': 1000, 'bwind': 0.0, 'lambdaf': 0.0, 'mxns': 3.0, 'beta': -1.0, 'tflag': 1, 'acc2': 1.5, 'grflag' : 1, 'remnantflag': 4, 'ceflag': 0, 'eddfac': 1.0, 'ifflag': 0, 'bconst': 3000, 'sigma': 265.0, 'gamma': -2.0, 'pisn': 45.0, 'natal_kick_array' : [[-100.0,-100.0,-100.0,-100.0,0.0], [-100.0,-100.0,-100.0,-100.0,0.0]], 'bhsigmafrac' : 1.0, 'polar_kick_angle' : 90, 'qcrit_array' : [0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0], 'cekickflag' : 2, 'cehestarflag' : 0, 'cemergeflag' : 0, 'ecsn' : 2.25, 'ecsn_mlow' : 1.6, 'aic' : 1, 'ussn' : 0, 'sigmadiv' :-20.0, 'qcflag' : 1, 'eddlimflag' : 0, 'fprimc_array' : [2.0/21.0,2.0/21.0,2.0/21.0,2.0/21.0,2.0/21.0,2.0/21.0,2.0/21.0,2.0/21.0,2.0/21.0,2.0/21.0,2.0/21.0,2.0/21.0,2.0/21.0,2.0/21.0,2.0/21.0,2.0/21.0], 'bhspinflag' : 0, 'bhspinmag' : 0.0, 'rejuv_fac' : 1.0, 'rejuvflag' : 0, 'htpmb' : 1, 'ST_cr' : 1, 'ST_tide' : 1, 'bdecayfac' : 1, 'rembar_massloss' : 0.5, 'kickflag' : 1, 'zsun' : 0.014, 'bhms_coll_flag' : 0, 'don_lim' : -1, 'acc_lim' : -1, 'rtmsflag' : 0, 'wd_mass_lim': 1}

Once the binary is initialized and the BSE model is set, the system is evolved with the the Evolve class, which calls the evolv2.f subroutine in the BSE source code.

In [6]: bpp, bcm, initC, kick_info = Evolve.evolve(initialbinarytable=single_binary, BSEDict=BSEDict)

For every evolved binary system, BSE generates two arrays, which are stored as pandas DataFrames in COSMIC:

bpp - contains binary parameters at important stages in the binary’s evolution, including stellar evolutionary phase changes or mass transfer episodes.
bcm - contains several binary parameters at user specified time steps during the binary’s evolution. The default setting in COSMIC is to output the final stage of the binary at the evolution time specified by the user.

You can see the different parameters included in each DataFrame using the columns attribute of the DataFrame:

In [7]: print(bpp.columns)
Index(['tphys', 'mass_1', 'mass_2', 'kstar_1', 'kstar_2', 'sep', 'porb', 'ecc',
       'RRLO_1', 'RRLO_2', 'evol_type', 'aj_1', 'aj_2', 'tms_1', 'tms_2',
       'massc_1', 'massc_2', 'rad_1', 'rad_2', 'mass0_1', 'mass0_2', 'lum_1',
       'lum_2', 'teff_1', 'teff_2', 'radc_1', 'radc_2', 'menv_1', 'menv_2',
       'renv_1', 'renv_2', 'omega_spin_1', 'omega_spin_2', 'B_1', 'B_2',
       'bacc_1', 'bacc_2', 'tacc_1', 'tacc_2', 'epoch_1', 'epoch_2',
       'bhspin_1', 'bhspin_2', 'bin_num'],
      dtype='object')

In [8]: print(bcm.columns)
Index(['tphys', 'kstar_1', 'mass0_1', 'mass_1', 'lum_1', 'rad_1', 'teff_1',
       'massc_1', 'radc_1', 'menv_1', 'renv_1', 'epoch_1', 'omega_spin_1',
       'deltam_1', 'RRLO_1', 'kstar_2', 'mass0_2', 'mass_2', 'lum_2', 'rad_2',
       'teff_2', 'massc_2', 'radc_2', 'menv_2', 'renv_2', 'epoch_2',
       'omega_spin_2', 'deltam_2', 'RRLO_2', 'porb', 'sep', 'ecc', 'B_1',
       'B_2', 'SN_1', 'SN_2', 'bin_state', 'merger_type', 'bin_num'],
      dtype='object')

The units are broadly consistent with BSE and are described in Understanding COSMIC outputs.

The evol_type column in bpp indicates the evolutionary change that occurred for each line. The meaning of each number is described here, Evolve Type.

Each of the parameters in bpp or bcm can be accessed in the usual way for DataFrames:

In [9]: bpp.mass_1
Out[9]: 
   85.543645
   72.720332
   72.589218
   72.504133
   72.504133
  114.941973
  114.827105
   64.562553
   62.286136
   44.500000
   44.500000
Name: mass_1, dtype: float64

In [10]: bpp = bpp[['mass_1', 'mass_2', 'kstar_1', 'kstar_2', 'sep', 'evol_type']]

You can use the utils.convert_kstar_evol_type function to convert the kstar_1, kstar_2, and evol_type columns from integers to strings that describe each int:

In [11]: from cosmic.utils import convert_kstar_evol_type

In [12]: convert_kstar_evol_type(bpp)
Out[12]: 
       mass_1     mass_2  ...          sep                  evol_type
0   85.543645  84.997840  ...  1363.435508              initial state
0   72.720332  72.376321  ...  1602.531512               kstar change
0   72.589218  72.377845  ...   883.827396  begin Roche lobe overflow
0   72.504133  72.396628  ...   883.794364                    contact
0   72.504133  72.396628  ...   883.794364      begin common envelope
0  114.941973  72.396628  ...     0.000000        end common envelope
0  114.827105   0.000000  ...     0.000000               kstar change
0   64.562553   0.000000  ...     0.000000               kstar change
0   62.286136   0.000000  ...     0.000000       supernova of primary
0   44.500000   0.000000  ...     0.000000               kstar change
0   44.500000   0.000000  ...     0.000000         max evolution time

[11 rows x 6 columns]

Note that utils.convert_kstar_evol_type is only applicable to the bpp array.

You can also use the built-in plotting function to see how the system evolves:

In [13]: from cosmic.plotting import evolve_and_plot

In [14]: single_binary = InitialBinaryTable.InitialBinaries(m1=85.543645, m2=84.99784, porb=446.795757,
   ....:                                                    ecc=0.448872, tphysf=13700.0,
   ....:                                                    kstar1=1, kstar2=1, metallicity=0.002)
   ....: 

In [15]: BSEDict = {'xi': 1.0, 'bhflag': 1, 'neta': 0.5, 'windflag': 3, 'wdflag': 1, 'alpha1': 1.0, 'pts1': 0.001, 'pts3': 0.02, 'pts2': 0.01, 'epsnov': 0.001, 'hewind': 0.5, 'ck': 1000, 'bwind': 0.0, 'lambdaf': 0.0, 'mxns': 3.0, 'beta': -1.0, 'tflag': 1, 'acc2': 1.5, 'grflag' : 1, 'remnantflag': 4, 'ceflag': 0, 'eddfac': 1.0, 'ifflag': 0, 'bconst': 3000, 'sigma': 265.0, 'gamma': -2.0, 'pisn': 45.0, 'natal_kick_array' : [[-100.0,-100.0,-100.0,-100.0,0.0], [-100.0,-100.0,-100.0,-100.0,0.0]], 'bhsigmafrac' : 1.0, 'polar_kick_angle' : 90, 'qcrit_array' : [0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0], 'cekickflag' : 2, 'cehestarflag' : 0, 'cemergeflag' : 0, 'ecsn' : 2.25, 'ecsn_mlow' : 1.6, 'aic' : 1, 'ussn' : 0, 'sigmadiv' :-20.0, 'qcflag' : 1, 'eddlimflag' : 0, 'fprimc_array' : [2.0/21.0,2.0/21.0,2.0/21.0,2.0/21.0,2.0/21.0,2.0/21.0,2.0/21.0,2.0/21.0,2.0/21.0,2.0/21.0,2.0/21.0,2.0/21.0,2.0/21.0,2.0/21.0,2.0/21.0,2.0/21.0], 'bhspinflag' : 0, 'bhspinmag' : 0.0, 'rejuv_fac' : 1.0, 'rejuvflag' : 0, 'htpmb' : 1, 'ST_cr' : 1, 'ST_tide' : 1, 'bdecayfac' : 1, 'rembar_massloss' : 0.5, 'kickflag' : 1, 'zsun' : 0.014, 'bhms_coll_flag' : 0, 'don_lim' : -1, 'acc_lim' : -1, 'rtmsflag' : 0, 'wd_mass_lim': 1}

In [16]: fig = evolve_and_plot(single_binary, t_min=None, t_max=None, BSEDict=BSEDict, sys_obs={})

(Source code, png, hires.png, pdf)

In this case, all the action happens in the first few Myr, so let’s specify a t_max:

In [17]: fig = evolve_and_plot(initC, t_min=None, t_max=6.0, BSEDict={}, sys_obs={})

(Source code, png, hires.png, pdf)