Using COSMIC to evolve binaries¶
COSMIC can evolve binaries for several different use cases. Below you’ll find examples to run a single binary system, multiple binary systems or a grid of binaries.
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
m2 : ZAMS mass of the secondary star in
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.
)
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 porb ecc metallicity tphysf mass0_1 ... tacc_2 epoch_1 epoch_2 tms_1 tms_2 bhspin_1 bhspin_2 tphys binfrac
0 1.0 1.0 85.543645 84.99784 446.795757 0.448872 0.002 13700.0 85.543645 ... 0.0 0.0 0.0 0.0 0.0 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.5, 'mxns': 2.5, 'beta': 0.125, 'tflag': 1, 'acc2': 1.5, 'remnantflag': 3, 'ceflag': 0, 'eddfac': 1.0, 'ifflag': 0, 'bconst': 3000, 'sigma': 265.0, 'gamma': -1.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.5, 'ecsn_mlow' : 1.4, 'aic' : 1, 'ussn' : 0, 'sigmadiv' :-20.0, 'qcflag' : 2, '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' : 0, 'bdecayfac' : 1, 'rembar_massloss' : 0.5, 'kickflag' : 0, 'zsun' : 0.014}
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 Describing the output of COSMIC/BSE: Columns names/Values/Units.
The evol_type column in bpp indicates the evolutionary change that occured 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]:
0 85.543645
0 77.610588
0 77.488343
0 59.516006
0 35.234327
0 35.041304
0 27.280996
0 26.780996
0 26.799354
0 26.801771
0 26.809532
0 26.810556
0 26.811851
0 26.812899
0 26.812899
0 26.812899
Name: mass_1, dtype: float64
In [10]: bpp[['mass_1', 'mass_2', 'kstar_1', 'kstar_2', 'sep', 'evol_type']]
Out[10]:
mass_1 mass_2 kstar_1 kstar_2 sep evol_type
0 85.543645 84.997840 1.0 1.0 1363.435508 1.0
0 77.610588 77.204563 2.0 1.0 1500.029741 2.0
0 77.488343 77.208726 2.0 1.0 827.546964 3.0
0 59.516006 95.001610 4.0 1.0 902.261909 2.0
0 35.234327 118.500023 4.0 1.0 1523.445555 4.0
0 35.041304 118.563115 7.0 1.0 1522.926443 2.0
0 27.280996 116.368926 7.0 1.0 1628.394867 15.0
0 26.780996 116.368926 14.0 1.0 1634.102529 2.0
0 26.799354 88.573419 14.0 2.0 2024.417588 2.0
0 26.801771 88.370688 14.0 2.0 2020.237661 3.0
0 26.809532 41.973522 14.0 2.0 476.120930 4.0
0 26.810556 41.958887 14.0 4.0 476.178665 2.0
0 26.811851 41.885153 14.0 7.0 476.635279 2.0
0 26.812899 31.733602 14.0 7.0 559.250334 16.0
0 26.812899 31.233602 14.0 14.0 564.109451 2.0
0 26.812899 31.233602 14.0 14.0 563.781534 10.0
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]:
tphys mass_1 mass_2 kstar_1 kstar_2 sep porb ... tacc_1 tacc_2 epoch_1 epoch_2 bhspin_1 bhspin_2 bin_num
0 0.000000 85.543645 84.997840 MS, > 0.7 M⊙ MS, > 0.7 M⊙ 1363.435508 446.795757 ... 0.0 0.0 0.000000 0.000000 0.0 0.0 0
0 3.690858 77.610588 77.204563 Hertzsprung Gap MS, > 0.7 M⊙ 1500.029741 541.146974 ... 0.0 0.0 -0.069029 -0.068690 0.0 0.0 0
0 3.692008 77.488343 77.208726 Hertzsprung Gap MS, > 0.7 M⊙ 827.546964 221.829904 ... 0.0 0.0 -0.069029 -0.068620 0.0 0.0 0
0 3.693628 59.516006 95.001610 Core Helium Burning MS, > 0.7 M⊙ 902.261909 252.686532 ... 0.0 0.0 -0.069029 0.819837 0.0 0.0 0
0 3.704436 35.234327 118.500023 Core Helium Burning MS, > 0.7 M⊙ 1523.445555 555.810963 ... 0.0 0.0 -0.069029 1.488651 0.0 0.0 0
0 3.705723 35.041304 118.563115 Naked Helium Star MS MS, > 0.7 M⊙ 1522.926443 555.761805 ... 0.0 0.0 3.695480 1.488886 0.0 0.0 0
0 4.030606 27.280996 116.368926 Naked Helium Star MS MS, > 0.7 M⊙ 1628.394867 635.417321 ... 0.0 0.0 3.670998 1.479874 0.0 0.0 0
0 4.030606 26.780996 116.368926 Black Hole MS, > 0.7 M⊙ 1634.102529 639.875607 ... 0.0 0.0 4.030606 1.479874 0.0 0.0 0
0 4.865319 26.799354 88.573419 Black Hole Hertzsprung Gap 2024.417588 982.812441 ... 0.0 0.0 4.030606 1.258575 0.0 0.0 0
0 4.866836 26.801771 88.370688 Black Hole Hertzsprung Gap 2020.237661 980.621779 ... 0.0 0.0 4.030606 1.258575 0.0 0.0 0
0 4.867668 26.809532 41.973522 Black Hole Hertzsprung Gap 476.120930 145.180184 ... 0.0 0.0 4.030606 1.258575 0.0 0.0 0
0 4.867798 26.810556 41.958887 Black Hole Core Helium Burning 476.178665 145.220960 ... 0.0 0.0 4.030606 1.258575 0.0 0.0 0
0 4.881760 26.811851 41.885153 Black Hole Naked Helium Star MS 476.635279 145.506548 ... 0.0 0.0 4.030606 4.870512 0.0 0.0 0
0 5.177980 26.812899 31.733602 Black Hole Naked Helium Star MS 559.250334 200.323115 ... 0.0 0.0 4.030606 4.845725 0.0 0.0 0
0 5.177980 26.812899 31.233602 Black Hole Black Hole 564.109451 203.811742 ... 0.0 0.0 4.030606 5.177980 0.0 0.0 0
0 13700.000000 26.812899 31.233602 Black Hole Black Hole 563.781534 203.634054 ... 0.0 0.0 4.030606 5.177980 0.0 0.0 0
[16 rows x 44 columns]
Note that utils.convert_kstar_evol_type is only applicable to the bpp
array.
multiple binaries¶
Multiple systems can also be initialized and evolved; below is an example for systems that could form GW150914 and GW170817 - like binaries.
In [13]: binary_set = InitialBinaryTable.InitialBinaries(m1=[85.543645, 11.171469], m2=[84.99784, 6.67305], porb=[446.795757, 170.758343], ecc=[0.448872, 0.370], tphysf=[13700.0, 13700.0], kstar1=[1, 1], kstar2=[1, 1], metallicity=[0.002, 0.02])
In [14]: print(binary_set)
kstar_1 kstar_2 mass_1 mass_2 porb ecc metallicity tphysf mass0_1 ... tacc_2 epoch_1 epoch_2 tms_1 tms_2 bhspin_1 bhspin_2 tphys binfrac
0 1.0 1.0 85.543645 84.99784 446.795757 0.448872 0.002 13700.0 85.543645 ... 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0
1 1.0 1.0 11.171469 6.67305 170.758343 0.370000 0.020 13700.0 11.171469 ... 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0
[2 rows x 38 columns]
In [15]: import numpy as np
In [16]: np.random.seed(5)
In [17]: bpp, bcm, initC, kick_info = Evolve.evolve(initialbinarytable=binary_set, BSEDict=BSEDict)
Note that the BSEDict did not be reinitialized since the BSE model did not change.
As before, bpp, bcm, and initC are returned as pandas DataFrames which assign an index to each binary system we evolve. We can access each binary as follows:
In [18]: print(bpp.loc[0])
tphys mass_1 mass_2 kstar_1 kstar_2 sep porb ecc ... bacc_2 tacc_1 tacc_2 epoch_1 epoch_2 bhspin_1 bhspin_2 bin_num
0 0.000000 85.543645 84.997840 1.0 1.0 1363.435508 446.795757 0.448872 ... 0.0 0.0 0.0 0.000000 0.000000 0.0 0.0 0
0 3.690858 77.610588 77.204563 2.0 1.0 1500.029741 541.146974 0.448645 ... 0.0 0.0 0.0 -0.069029 -0.068690 0.0 0.0 0
0 3.692008 77.488343 77.208726 2.0 1.0 827.546964 221.829904 0.000000 ... 0.0 0.0 0.0 -0.069029 -0.068620 0.0 0.0 0
0 3.693628 59.516006 95.001610 4.0 1.0 902.261909 252.686532 0.000000 ... 0.0 0.0 0.0 -0.069029 0.819837 0.0 0.0 0
0 3.704436 35.234327 118.500023 4.0 1.0 1523.445555 555.810963 0.000000 ... 0.0 0.0 0.0 -0.069029 1.488651 0.0 0.0 0
0 3.705723 35.041304 118.563115 7.0 1.0 1522.926443 555.761805 0.000000 ... 0.0 0.0 0.0 3.695480 1.488886 0.0 0.0 0
0 4.030606 27.280996 116.368926 7.0 1.0 1628.394867 635.417321 0.000000 ... 0.0 0.0 0.0 3.670998 1.479874 0.0 0.0 0
0 4.030606 26.780996 116.368926 14.0 1.0 1634.102529 639.875607 0.003493 ... 0.0 0.0 0.0 4.030606 1.479874 0.0 0.0 0
0 4.865319 26.799354 88.573419 14.0 2.0 2024.417588 982.812441 0.003491 ... 0.0 0.0 0.0 4.030606 1.258575 0.0 0.0 0
0 4.866836 26.801771 88.370688 14.0 2.0 2020.237661 980.621779 0.000000 ... 0.0 0.0 0.0 4.030606 1.258575 0.0 0.0 0
0 4.867668 26.809532 41.973522 14.0 2.0 476.120930 145.180184 0.000000 ... 0.0 0.0 0.0 4.030606 1.258575 0.0 0.0 0
0 4.867798 26.810556 41.958887 14.0 4.0 476.178665 145.220960 0.000000 ... 0.0 0.0 0.0 4.030606 1.258575 0.0 0.0 0
0 4.881760 26.811851 41.885153 14.0 7.0 476.635279 145.506548 0.000000 ... 0.0 0.0 0.0 4.030606 4.870512 0.0 0.0 0
0 5.177980 26.812899 31.733602 14.0 7.0 559.250334 200.323115 0.000000 ... 0.0 0.0 0.0 4.030606 4.845725 0.0 0.0 0
0 5.177980 26.812899 31.233602 14.0 14.0 564.109451 203.811742 0.008614 ... 0.0 0.0 0.0 4.030606 5.177980 0.0 0.0 0
0 13700.000000 26.812899 31.233602 14.0 14.0 563.781534 203.634054 0.008614 ... 0.0 0.0 0.0 4.030606 5.177980 0.0 0.0 0
[16 rows x 44 columns]
In [19]: print(bcm.loc[0])
tphys kstar_1 mass0_1 mass_1 lum_1 rad_1 teff_1 massc_1 ... ecc B_1 B_2 SN_1 SN_2 bin_state merger_type bin_num
0 0.0 1.0 85.543645 85.543645 9.416137e+05 11.060285 54306.854457 0.000000 ... 0.448872 0.0 0.0 0.0 0.0 0 -001 0
0 13700.0 14.0 27.280996 26.812899 1.000000e-10 0.000114 1719.549945 26.812899 ... 0.008614 0.0 0.0 1.0 1.0 0 -001 0
[2 rows x 39 columns]
In [20]: print(initC.loc[0])
kstar_1 1.000000
kstar_2 1.000000
mass_1 85.543645
mass_2 84.997840
porb 446.795757
...
fprimc_11 0.095238
fprimc_12 0.095238
fprimc_13 0.095238
fprimc_14 0.095238
fprimc_15 0.095238
Name: 0, Length: 132, dtype: float64
In [21]: print(bpp.loc[1])
tphys mass_1 mass_2 kstar_1 kstar_2 sep porb ecc ... bacc_2 tacc_1 tacc_2 epoch_1 epoch_2 bhspin_1 bhspin_2 bin_num
1 0.000000 11.171469 6.673050 1.0 1.0 338.356712 170.758343 0.370000 ... 0.0 0.000000 0.0 0.000000 0.000000 0.0 0.0 1
1 19.429462 10.767530 6.665578 2.0 1.0 346.228856 178.825577 0.369953 ... 0.0 0.000000 0.0 -0.968174 -0.024232 0.0 0.0 1
1 19.463620 10.765398 6.665623 2.0 1.0 218.163264 89.450584 0.000000 ... 0.0 0.000000 0.0 -0.975415 -0.023942 0.0 0.0 1
1 19.479182 4.872815 8.954552 3.0 1.0 291.471876 155.094861 0.000000 ... 0.0 0.000000 0.0 -1.129457 11.666434 0.0 0.0 1
1 19.486113 3.483800 10.343160 4.0 1.0 428.057952 276.034051 0.000000 ... 0.0 0.000000 0.0 -1.129457 14.337154 0.0 0.0 1
1 19.496870 2.431140 11.395311 4.0 1.0 720.157365 602.362478 0.000000 ... 0.0 0.000000 0.0 -1.129457 15.553354 0.0 0.0 1
1 19.526722 2.424998 11.400067 7.0 1.0 719.624887 601.724693 0.000000 ... 0.0 0.000000 0.0 19.471002 15.556206 0.0 0.0 1
1 22.918450 2.235745 11.395898 8.0 1.0 729.834943 618.920750 0.000000 ... 0.0 0.000000 0.0 19.160262 15.552576 0.0 0.0 1
1 23.250117 2.156953 11.402081 8.0 1.0 732.890348 624.476774 0.000000 ... 0.0 0.000000 0.0 19.160262 15.559809 0.0 0.0 1
1 23.250117 1.427227 11.402081 13.0 1.0 1010.152196 1038.845024 0.514914 ... 0.0 0.000000 0.0 23.250117 15.559809 0.0 0.0 1
1 34.259925 1.427233 10.996009 13.0 2.0 1042.564056 1106.901082 0.514912 ... 0.0 11.009808 0.0 23.250117 14.610607 0.0 0.0 1
1 34.307074 1.427239 10.991812 13.0 3.0 1013.350998 1060.883756 0.500529 ... 0.0 11.056957 0.0 23.250117 14.597304 0.0 0.0 1
1 34.308349 1.427241 10.991568 13.0 3.0 506.171779 374.523517 0.000000 ... 0.0 11.058309 0.0 23.250117 14.597304 0.0 0.0 1
1 34.308349 1.427241 10.991568 13.0 3.0 506.171779 374.523517 0.000000 ... 0.0 11.058309 0.0 23.250117 14.597304 0.0 0.0 1
1 34.308349 1.427241 2.495020 13.0 7.0 6.064825 0.874040 0.000000 ... 0.0 11.058309 0.0 23.250117 14.597304 0.0 0.0 1
1 34.308349 1.427241 2.495020 13.0 7.0 6.064825 0.874040 0.000000 ... 0.0 11.058309 0.0 23.250117 34.308349 0.0 0.0 1
1 37.598485 1.429600 2.290705 13.0 8.0 6.372670 0.966641 0.000000 ... 0.0 14.348444 0.0 23.250117 34.003825 0.0 0.0 1
1 37.836980 1.431201 2.249341 13.0 8.0 6.426085 0.984092 0.000000 ... 0.0 14.587139 0.0 23.250117 34.003825 0.0 0.0 1
1 37.841823 1.431480 1.716722 13.0 9.0 6.146256 0.995306 0.000000 ... 0.0 14.591982 0.0 23.250117 34.003825 0.0 0.0 1
1 37.841823 1.431480 1.716722 13.0 9.0 6.146256 0.995306 0.000000 ... 0.0 14.591982 0.0 23.250117 34.003825 0.0 0.0 1
1 37.841823 1.431480 1.333588 13.0 12.0 5.773101 0.966792 0.000000 ... 0.0 14.591982 0.0 23.250117 34.003825 0.0 0.0 1
1 37.841823 1.431480 1.333588 13.0 12.0 5.773101 0.966792 0.000000 ... 0.0 14.591982 0.0 23.250117 37.841823 0.0 0.0 1
1 13700.000000 1.431480 1.333588 13.0 12.0 5.020490 0.784040 0.000000 ... 0.0 14.591982 0.0 23.250117 37.841823 0.0 0.0 1
[23 rows x 44 columns]
grid of binaries¶
Sometimes it is helpful to run a grid of initial binaries to explore how changing a single paramter affects the evolved binary. Here we evolve the same system that produces a GW150914-like binary, but run over several initial orbital periods spaced evenly in log space.
In [22]: n_grid = 10
In [23]: binary_grid = InitialBinaryTable.InitialBinaries(m1=np.ones(n_grid)*100.0, m2=np.ones(n_grid)*85.0, porb=np.logspace(3,5,n_grid), ecc=np.ones(n_grid)*0.65, tphysf=np.ones(n_grid)*13700.0, kstar1=np.ones(n_grid), kstar2=np.ones(n_grid), metallicity=np.ones(n_grid)*0.005)
In [24]: print(binary_grid)
kstar_1 kstar_2 mass_1 mass_2 porb ecc metallicity tphysf mass0_1 ... tacc_2 epoch_1 epoch_2 tms_1 tms_2 bhspin_1 bhspin_2 tphys binfrac
0 1.0 1.0 100.0 85.0 1000.000000 0.65 0.005 13700.0 100.0 ... 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0
1 1.0 1.0 100.0 85.0 1668.100537 0.65 0.005 13700.0 100.0 ... 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0
2 1.0 1.0 100.0 85.0 2782.559402 0.65 0.005 13700.0 100.0 ... 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0
3 1.0 1.0 100.0 85.0 4641.588834 0.65 0.005 13700.0 100.0 ... 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0
4 1.0 1.0 100.0 85.0 7742.636827 0.65 0.005 13700.0 100.0 ... 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0
5 1.0 1.0 100.0 85.0 12915.496650 0.65 0.005 13700.0 100.0 ... 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0
6 1.0 1.0 100.0 85.0 21544.346900 0.65 0.005 13700.0 100.0 ... 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0
7 1.0 1.0 100.0 85.0 35938.136638 0.65 0.005 13700.0 100.0 ... 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0
8 1.0 1.0 100.0 85.0 59948.425032 0.65 0.005 13700.0 100.0 ... 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0
9 1.0 1.0 100.0 85.0 100000.000000 0.65 0.005 13700.0 100.0 ... 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0
[10 rows x 38 columns]
In [25]: bpp, bcm, initC, kick_info = Evolve.evolve(initialbinarytable=binary_grid, BSEDict=BSEDict)
In [26]: print(bpp)
tphys mass_1 mass_2 kstar_1 kstar_2 sep porb ecc ... bacc_2 tacc_1 tacc_2 epoch_1 epoch_2 bhspin_1 bhspin_2 bin_num
0 0.000000 100.000000 85.000000 1.0 1.0 2397.042879 1.000000e+03 0.650000 ... 0.0 0.0 0.0 0.000000 0.000000 0.0 0.0 0
0 3.532986 75.202022 72.179052 2.0 1.0 2998.714278 1.567667e+03 0.649328 ... 0.0 0.0 0.0 -0.207963 -0.112228 0.0 0.0 0
0 3.534807 74.972259 72.182104 2.0 1.0 1053.002253 3.264597e+02 0.000000 ... 0.0 0.0 0.0 -0.207963 -0.112173 0.0 0.0 0
0 3.536082 61.736810 85.266246 4.0 1.0 1088.444387 3.432562e+02 0.000000 ... 0.0 0.0 0.0 -0.207963 0.606249 0.0 0.0 0
0 3.586119 34.343506 108.570167 4.0 1.0 1533.848480 5.823834e+02 0.000000 ... 0.0 0.0 0.0 -0.207963 1.297939 0.0 0.0 0
.. ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ...
9 3.900117 28.549314 38.700770 14.0 4.0 138349.274758 7.272623e+05 0.644507 ... 0.0 0.0 0.0 3.900117 -0.132272 0.0 0.0 9
9 3.938406 28.556114 33.951555 14.0 7.0 148739.036119 8.408973e+05 0.644363 ... 0.0 0.0 0.0 3.900117 3.730108 0.0 0.0 9
9 4.057944 28.556114 26.971840 14.0 7.0 167447.745004 1.065701e+06 0.644363 ... 0.0 0.0 0.0 3.900117 3.696169 0.0 0.0 9
9 4.057944 28.556114 26.471840 14.0 14.0 169117.761402 1.086587e+06 0.644762 ... 0.0 0.0 0.0 3.900117 4.057944 0.0 0.0 9
9 13700.000000 28.556114 26.471840 14.0 14.0 169021.306377 1.085657e+06 0.644762 ... 0.0 0.0 0.0 3.900117 4.057944 0.0 0.0 9
[143 rows x 44 columns]
In [27]: print(bcm)
tphys kstar_1 mass0_1 mass_1 lum_1 rad_1 teff_1 massc_1 ... ecc B_1 B_2 SN_1 SN_2 bin_state merger_type bin_num
0 0.0 1.0 100.000000 100.000000 1.240906e+06 13.897608 51907.996879 0.000000 ... 0.650000 0.0 0.0 0.0 0.0 0 -001 0
0 13700.0 14.0 21.578953 21.107916 1.000000e-10 0.000089 1938.045651 21.107916 ... 0.012188 0.0 0.0 1.0 1.0 0 -001 0
1 0.0 1.0 100.000000 100.000000 1.240906e+06 13.897608 51907.996879 0.000000 ... 0.650000 0.0 0.0 0.0 0.0 0 -001 1
1 13700.0 14.0 22.433369 21.956872 1.000000e-10 0.000093 1900.209353 21.956872 ... 0.011969 0.0 0.0 1.0 1.0 0 -001 1
2 0.0 1.0 100.000000 100.000000 1.240906e+06 13.897608 51907.996879 0.000000 ... 0.650000 0.0 0.0 0.0 0.0 0 -001 2
2 13700.0 14.0 22.923370 22.442103 1.000000e-10 0.000095 1879.554414 22.442103 ... 0.011805 0.0 0.0 1.0 1.0 0 -001 2
3 0.0 1.0 100.000000 100.000000 1.240906e+06 13.897608 51907.996879 0.000000 ... 0.650000 0.0 0.0 0.0 0.0 0 -001 3
3 13700.0 14.0 23.935364 24.053205 1.000000e-10 0.000102 1815.516535 24.053205 ... 0.010906 0.0 0.0 1.0 1.0 0 -001 3
4 0.0 1.0 100.000000 100.000000 1.240906e+06 13.897608 51907.996879 0.000000 ... 0.650000 0.0 0.0 0.0 0.0 0 -001 4
4 13700.0 14.0 26.089665 27.091764 1.000000e-10 0.000115 1710.677084 27.091764 ... 0.009711 0.0 0.0 1.0 1.0 0 -001 4
5 0.0 1.0 100.000000 100.000000 1.240906e+06 13.897608 51907.996879 0.000000 ... 0.650000 0.0 0.0 0.0 0.0 0 -001 5
5 13700.0 14.0 29.700970 29.480986 1.000000e-10 0.000125 1639.893634 29.480986 ... 0.395195 0.0 0.0 1.0 1.0 0 -001 5
6 0.0 1.0 100.000000 100.000000 1.240906e+06 13.897608 51907.996879 0.000000 ... 0.650000 0.0 0.0 0.0 0.0 0 -001 6
6 13700.0 14.0 29.275585 28.846795 1.000000e-10 0.000122 1657.821990 28.846795 ... 0.595454 0.0 0.0 1.0 1.0 0 -001 6
7 0.0 1.0 100.000000 100.000000 1.240906e+06 13.897608 51907.996879 0.000000 ... 0.650000 0.0 0.0 0.0 0.0 0 -001 7
7 13700.0 14.0 29.144942 28.674756 1.000000e-10 0.000122 1662.787764 28.674756 ... 0.640203 0.0 0.0 1.0 1.0 0 -001 7
8 0.0 1.0 100.000000 100.000000 1.240906e+06 13.897608 51907.996879 0.000000 ... 0.650000 0.0 0.0 0.0 0.0 0 -001 8
8 13700.0 14.0 29.082199 28.596065 1.000000e-10 0.000121 1665.074018 28.596065 ... 0.636038 0.0 0.0 1.0 1.0 0 -001 8
9 0.0 1.0 100.000000 100.000000 1.240906e+06 13.897608 51907.996879 0.000000 ... 0.650000 0.0 0.0 0.0 0.0 0 -001 9
9 13700.0 14.0 29.049314 28.556114 1.000000e-10 0.000121 1666.238366 28.556114 ... 0.644762 0.0 0.0 1.0 1.0 0 -001 9
[20 rows x 39 columns]
dynamically set time resolution for bcm array¶
COSMIC has the ability to set time resolution of the bcm array depending on the current state of the evolution. Below we demonstrate three scenarios, setting dtp only during mass transfer, setting dtp to the same resolution for all of the evolution except for after the system merges or is disrupted, and finally an example of setting dtp only during the HMXRB stage of the evoltuion.
First, print all time steps during mass transfer
In [28]: single_binary = InitialBinaryTable.InitialBinaries(m1=7.806106, m2=5.381412, porb=2858.942021, ecc=0.601408, tphysf=13700.0, kstar1=1, kstar2=1, metallicity=0.02)
In [29]: 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.5, 'mxns': 3.0, 'beta': 0.125, 'tflag': 1, 'acc2': 1.5, 'remnantflag': 3, 'ceflag': 0, 'eddfac': 1.0, 'ifflag': 0, 'bconst': 3000, 'sigma': 265.0, 'gamma': -1.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' : [1000.0,1000.0,1000.0,1000.0,1000.0,1000.0,1000.0,1000.0,1000.0, 1000.0,1000.0,1000.0,1000.0,1000.0,1000.0,1000.0],'cekickflag' : 2, 'cehestarflag' : 0, 'cemergeflag' : 0, 'ecsn' : 2.5, 'ecsn_mlow' : 1.4, 'aic' : 1, 'ussn' : 0, 'sigmadiv' :-20.0, 'qcflag' : 2, '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' : 0, 'bdecayfac' : 1, 'rembar_massloss' : 0.5, 'kickflag' : 0, 'zsun' : 0.014}
In [30]: bpp, bcm, initC, kick_info = Evolve.evolve(initialbinarytable=single_binary, BSEDict=BSEDict, timestep_conditions =[['RRLO_1>=1', 'dtp=0.0'], ['RRLO_2>=1', 'dtp=0.0']])
In [31]: print(bcm[['tphys', 'kstar_1', 'kstar_2', 'mass_1', 'mass_2', 'RRLO_1', 'RRLO_2']])
tphys kstar_1 kstar_2 mass_1 mass_2 RRLO_1 RRLO_2
0 0.0 1.0 1.0 7.806106 5.381412 0.010953 0.010459
0 13700.0 12.0 13.0 1.358765 1.427227 0.000100 0.000100
Second, pick a certain resolution for the bcm array until the system mergers or is disrutped and then only print the final state
In [32]: bpp, bcm, initC, kick_info = Evolve.evolve(initialbinarytable=single_binary, BSEDict=BSEDict, timestep_conditions =[['binstate=0', 'dtp=1.0']])
In [33]: print(bcm[['tphys', 'kstar_1', 'kstar_2', 'mass_1', 'mass_2', 'bin_state']])
tphys kstar_1 kstar_2 mass_1 mass_2 bin_state
0 0.0 1.0 1.0 7.806106 5.381412 0
0 1.0 1.0 1.0 7.805049 5.381334 0
0 2.0 1.0 1.0 7.803970 5.381256 0
0 3.0 1.0 1.0 7.802864 5.381177 0
0 4.0 1.0 1.0 7.801731 5.381098 0
.. ... ... ... ... ... ...
0 60.0 12.0 4.0 1.355383 9.735468 0
0 61.0 12.0 4.0 1.356135 9.493343 0
0 62.0 12.0 4.0 1.356355 9.334358 0
0 63.0 12.0 13.0 1.358782 1.427227 2
0 13700.0 12.0 13.0 1.358782 1.427227 2
[65 rows x 6 columns]
Finally, we show how to print a fine resolution only during the HMXRB stage of the evolution.
In [34]: 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 [35]: 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.5, 'mxns': 2.5, 'beta': 0.125, 'tflag': 1, 'acc2': 1.5, 'remnantflag': 3, 'ceflag': 0, 'eddfac': 1.0, 'ifflag': 0, 'bconst': 3000, 'sigma': 265.0, 'gamma': -1.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.5, 'ecsn_mlow' : 1.4, 'aic' : 1, 'ussn' : 0, 'sigmadiv' :-20.0, 'qcflag' : 2, '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' : 0, 'bdecayfac' : 1, 'rembar_massloss' : 0.5, 'kickflag': 0, 'zsun' : 0.014}
In [36]: bpp, bcm, initC, kick_info = Evolve.evolve(initialbinarytable=single_binary, BSEDict=BSEDict, timestep_conditions =[['kstar_1=14', 'kstar_2<10','dtp=0.1'], ['kstar_2=14', 'kstar_1<10','dtp=0.1']])
In [37]: print(bcm[['tphys', 'kstar_1', 'kstar_2', 'mass_1', 'mass_2', 'bin_state']])
tphys kstar_1 kstar_2 mass_1 mass_2 bin_state
0 0.000000 1.0 1.0 85.543645 84.997840 0
0 4.030606 14.0 1.0 26.780996 116.368926 0
0 4.130606 14.0 1.0 26.781133 115.587281 0
0 4.230606 14.0 1.0 26.781303 114.762463 0
0 4.330606 14.0 1.0 26.781516 113.906980 0
0 4.430606 14.0 1.0 26.781783 113.042355 0
0 4.530606 14.0 1.0 26.782114 112.201659 0
0 4.630606 14.0 1.0 26.783100 110.458173 0
0 4.730606 14.0 1.0 26.790076 99.938912 0
0 4.830606 14.0 1.0 26.798307 89.741635 0
0 4.930606 14.0 7.0 26.812023 40.196096 0
0 5.030606 14.0 7.0 26.812388 36.693634 0
0 5.130606 14.0 7.0 26.812742 33.280292 0
0 5.230606 14.0 14.0 26.812895 31.233489 0
0 13700.000000 14.0 14.0 26.812895 31.233489 0
restarting a binary¶
COSMIC allows you to restart a binary from any point in its evolution from a COSMIC generated bpp array. Below we provide an example of the same evolutionary track started from the beginning and three different points in the evolution, once sometime between the beginning and the first object going supernova, once between the first and second supernova, and finally after both supernova.:
>>> single_binary = InitialBinaryTable.InitialBinaries(m1=25.543645, m2=20.99784, porb=446.795757, ecc=0.448872, tphysf=13700.0, kstar1=1, kstar2=1, metallicity=0.002)
>>> 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.5, 'mxns': 3.0, 'beta': 0.125, 'tflag': 1, 'acc2': 1.5, 'remnantflag': 3, 'ceflag': 0, 'eddfac': 1.0, 'ifflag': 0, 'bconst': 3000, 'sigma': 265.0, 'gamma': -1.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.5, 'ecsn_mlow' : 1.4, 'aic' : 1, 'ussn' : 0, 'sigmadiv' :-20.0, 'qcflag' : 2, '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' : 0, 'bdecayfac' : 1, 'randomseed' : -1235453, 'rembar_massloss' : 0.5, 'kickflag' : 0, 'zsun' : 0.014}
>>> for i in [3, 7, 11]:
>>> bpp, bcm, initC, kick_info = Evolve.evolve(initialbinarytable=single_binary, BSEDict=BSEDict)
>>> for column in bpp.columns:
>>> initC = initC.assign(**{column:bpp.iloc[i][column]})
>>> bpp_mid, bcm_mid, initC_mid, kick_info = Evolve.evolve(initialbinarytable=initC)
>>> if i == 3:
>>> print("From beginning")
>>> print(bpp)
>>> print("Started in middle at Index {0}".format(i))
>>> print(bpp_mid)