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
  71.959504
  32.352601
  32.176152
  25.488585
  24.988585
  24.988590
  24.989628
  24.990676
  24.990676
  24.990676
  24.990676
  25.001441
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
85.543645   84.997840  ...  1363.435508                     initial state
72.720332   72.376321  ...  1602.531512                      kstar change
72.589218   72.377845  ...   883.827396         begin Roche lobe overflow
71.959504   72.806393  ...   882.511619                      kstar change
32.352601  110.573279  ...  1614.251471            end Roche lobe overlow
32.176152  110.618802  ...  1614.063847                      kstar change
25.488585  106.891756  ...  1741.069066              supernova of primary
24.988585  106.891756  ...  1747.695130                      kstar change
24.988590   88.885767  ...  2023.984648                      kstar change
24.989628   88.681486  ...  2019.531513         begin Roche lobe overflow
24.990676   88.469784  ...  2018.508488                      kstar change
24.990676   88.469784  ...  2018.508488             begin common envelope
24.990676   41.981622  ...   320.677863               end common envelope
24.990676   41.981622  ...   320.677863            end Roche lobe overlow
25.001441   41.981622  ...          NaN  RLOF interpolation timeout error

[15 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)