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 ... 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' : 0, '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 Describing the output of COSMIC/BSE: Columns names/Values/Units.
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]:
0 85.543645
0 72.720332
0 72.589218
0 71.959504
0 32.352601
0 32.176152
0 25.488585
0 24.988585
0 24.988590
0 24.989628
0 24.990676
0 24.990676
0 24.990676
0 24.990676
0 24.990687
0 24.990687
0 24.990687
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 71.959504 72.806393 ... 882.511619 kstar change
0 32.352601 110.573279 ... 1614.251471 end Roche lobe overlow
0 32.176152 110.618802 ... 1614.063847 kstar change
0 25.488585 106.891756 ... 1741.069066 supernova of primary
0 24.988585 106.891756 ... 1747.695130 kstar change
0 24.988590 88.885767 ... 2023.280971 kstar change
0 24.989628 88.681512 ... 2018.828817 begin Roche lobe overflow
0 24.990676 88.469775 ... 2017.805217 kstar change
0 24.990676 88.469775 ... 2017.805217 begin common envelope
0 24.990676 41.981621 ... 320.660165 end common envelope
0 24.990676 41.981621 ... 320.660165 end Roche lobe overlow
0 24.990687 31.554344 ... 379.806762 supernova of secondary
0 24.990687 31.054344 ... 383.225670 kstar change
0 24.990687 31.054344 ... 382.989460 max evolution time
[17 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' : 0, '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)
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 [18]: 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 [19]: print(binary_set)
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 1.0 1.0 11.171469 6.67305 ... 0.0 0.0 0.0 1.0
[2 rows x 38 columns]
In [20]: import numpy as np
In [21]: np.random.seed(5)
In [22]: 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 [23]: print(bpp.loc[0])
tphys mass_1 mass_2 ... bhspin_1 bhspin_2 bin_num
0 0.000000 85.543645 84.997840 ... 0.0 0.0 0
0 3.717070 72.719783 72.376037 ... 0.0 0.0 0
0 3.718368 72.588666 72.377556 ... 0.0 0.0 0
0 3.720007 71.958903 72.806115 ... 0.0 0.0 0
0 3.739911 32.346795 110.583971 ... 0.0 0.0 0
0 3.741053 32.175552 110.600277 ... 0.0 0.0 0
0 4.071364 25.486398 106.874794 ... 0.0 0.0 0
0 4.071364 24.986398 106.874794 ... 0.0 0.0 0
0 4.894907 24.986403 88.691971 ... 0.0 0.0 0
0 4.896430 24.987444 88.487242 ... 0.0 0.0 0
0 4.897383 24.988491 88.276317 ... 0.0 0.0 0
0 4.897383 24.988491 88.276317 ... 0.0 0.0 0
0 4.897383 24.988491 41.858050 ... 0.0 0.0 0
0 4.897383 24.988491 41.858050 ... 0.0 0.0 0
0 5.206247 24.988501 31.476750 ... 0.0 0.0 0
0 5.206247 24.988501 30.976750 ... 0.0 0.0 0
0 205.190000 24.988501 30.976750 ... 0.0 0.0 0
[17 rows x 44 columns]
In [24]: print(bcm.loc[0])
tphys kstar_1 mass0_1 mass_1 ... SN_2 bin_state merger_type bin_num
0 0.00 1.0 85.543645 85.543645 ... 0.0 0 -001 0
0 0.01 1.0 85.527367 85.527367 ... 0.0 0 -001 0
0 0.02 1.0 85.511049 85.511049 ... 0.0 0 -001 0
0 0.03 1.0 85.494691 85.494691 ... 0.0 0 -001 0
0 0.04 1.0 85.478292 85.478292 ... 0.0 0 -001 0
.. ... ... ... ... ... ... ... ... ...
0 205.15 14.0 25.486398 24.988501 ... 1.0 0 -001 0
0 205.16 14.0 25.486398 24.988501 ... 1.0 0 -001 0
0 205.17 14.0 25.486398 24.988501 ... 1.0 0 -001 0
0 205.18 14.0 25.486398 24.988501 ... 1.0 0 -001 0
0 205.19 14.0 25.486398 24.988501 ... 1.0 0 -001 0
[20520 rows x 39 columns]
In [25]: 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: 138, dtype: float64
In [26]: print(bpp.loc[1])
tphys mass_1 mass_2 ... bhspin_1 bhspin_2 bin_num
1 0.000000 11.171469 6.673050 ... 0.0 0.0 1
1 19.427278 10.767909 6.664658 ... 0.0 0.0 1
1 19.461438 10.765777 6.664702 ... 0.0 0.0 1
1 19.476912 10.544417 6.884988 ... 0.0 0.0 1
1 19.476912 10.544417 6.884988 ... 0.0 0.0 1
1 19.476912 2.415442 6.884988 ... 0.0 0.0 1
1 19.476912 2.415442 6.884988 ... 0.0 0.0 1
1 22.952681 2.226294 6.887974 ... 0.0 0.0 1
1 23.231406 2.177806 6.896735 ... 0.0 0.0 1
1 23.254641 1.889676 7.181936 ... 0.0 0.0 1
1 23.256935 1.371085 7.700471 ... 0.0 0.0 1
1 23.258694 1.370344 7.700740 ... 0.0 0.0 1
1 47.641575 1.370365 7.628335 ... 0.0 0.0 1
1 47.673758 1.370368 7.628122 ... 0.0 0.0 1
1 47.673758 1.370368 7.628122 ... 0.0 0.0 1
1 47.673758 1.370344 8.188112 ... 0.0 0.0 1
1 47.785623 0.000000 7.998227 ... 0.0 0.0 1
1 47.785623 0.000000 1.277584 ... 0.0 0.0 1
1 247.770000 0.000000 1.277584 ... 0.0 0.0 1
[19 rows x 44 columns]
The plotting function can also take in multiple binaries. Let’s plot both the GW150914-like progenitor evolution and the GW170817-like progenitor evolutions. For the GW170817-like progenitor, we expect most of the evolution to take place in the first ~60 Myr.
In [27]: fig = evolve_and_plot(binary_set, t_min=None, t_max=[6.0, 60.0], BSEDict=BSEDict, sys_obs={})
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 [28]: n_grid = 10
In [29]: 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 [30]: 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' : 0, 'zsun' : 0.014, 'bhms_coll_flag' : 0, 'don_lim' : -1, 'acc_lim' : -1, 'rtmsflag' : 0, 'wd_mass_lim': 1}
In [31]: print(binary_grid)
kstar_1 kstar_2 mass_1 mass_2 ... bhspin_1 bhspin_2 tphys binfrac
0 1.0 1.0 100.0 85.0 ... 0.0 0.0 0.0 1.0
1 1.0 1.0 100.0 85.0 ... 0.0 0.0 0.0 1.0
2 1.0 1.0 100.0 85.0 ... 0.0 0.0 0.0 1.0
3 1.0 1.0 100.0 85.0 ... 0.0 0.0 0.0 1.0
4 1.0 1.0 100.0 85.0 ... 0.0 0.0 0.0 1.0
5 1.0 1.0 100.0 85.0 ... 0.0 0.0 0.0 1.0
6 1.0 1.0 100.0 85.0 ... 0.0 0.0 0.0 1.0
7 1.0 1.0 100.0 85.0 ... 0.0 0.0 0.0 1.0
8 1.0 1.0 100.0 85.0 ... 0.0 0.0 0.0 1.0
9 1.0 1.0 100.0 85.0 ... 0.0 0.0 0.0 1.0
[10 rows x 38 columns]
In [32]: bpp, bcm, initC, kick_info = Evolve.evolve(initialbinarytable=binary_grid, BSEDict=BSEDict)
In [33]: print(bpp)
tphys mass_1 mass_2 ... bhspin_1 bhspin_2 bin_num
0 0.000000 100.000000 85.000000 ... 0.0 0.0 0
0 3.561836 74.433283 67.900444 ... 0.0 0.0 0
0 3.563725 74.195874 67.901480 ... 0.0 0.0 0
0 3.564957 73.978352 67.965270 ... 0.0 0.0 0
0 3.656147 34.358643 98.037335 ... 0.0 0.0 0
.. ... ... ... ... ... ... ...
9 3.930319 28.153730 35.520384 ... 0.0 0.0 9
9 3.960762 28.159125 31.484483 ... 0.0 0.0 9
9 4.100892 28.159125 24.732495 ... 0.0 0.0 9
9 4.100892 28.159125 24.232495 ... 0.0 0.0 9
9 13700.000000 28.159125 24.232495 ... 0.0 0.0 9
[139 rows x 44 columns]
In [34]: print(bcm)
tphys kstar_1 mass0_1 ... bin_state merger_type bin_num
0 0.0 1.0 100.000000 ... 0 -001 0
0 13700.0 14.0 22.440552 ... 0 -001 0
1 0.0 1.0 100.000000 ... 0 -001 1
1 13700.0 14.0 23.651593 ... 0 -001 1
2 0.0 1.0 100.000000 ... 0 -001 2
2 13700.0 14.0 25.319576 ... 0 -001 2
3 0.0 1.0 100.000000 ... 0 -001 3
3 13700.0 14.0 27.066154 ... 0 -001 3
4 0.0 1.0 100.000000 ... 0 -001 4
4 13700.0 14.0 24.672076 ... 0 -001 4
5 0.0 1.0 100.000000 ... 0 -001 5
5 13700.0 14.0 29.517261 ... 0 -001 5
6 0.0 1.0 100.000000 ... 0 -001 6
6 13700.0 14.0 28.990549 ... 0 -001 6
7 0.0 1.0 100.000000 ... 0 -001 7
7 13700.0 14.0 28.809440 ... 0 -001 7
8 0.0 1.0 100.000000 ... 0 -001 8
8 13700.0 14.0 28.708663 ... 0 -001 8
9 0.0 1.0 100.000000 ... 0 -001 9
9 13700.0 14.0 28.653730 ... 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 HMB stage of the evolution.
First, print all time steps during mass transfer
In [35]: 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 [36]: 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' : 0, 'zsun' : 0.014, 'bhms_coll_flag' : 0, 'don_lim' : -1, 'acc_lim' : -1, 'rtmsflag' : 0, 'wd_mass_lim': 1}
In [37]: 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 [38]: 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.000000 1.0 1.0 7.806106 5.381412 0.010953 0.010459
0 43.565604 5.0 1.0 7.346069 5.396914 1.000222 0.008730
0 43.565604 8.0 1.0 2.037972 5.396914 0.520417 0.110930
0 43.579670 8.0 1.0 2.032299 5.399465 1.000951 0.110963
0 43.579674 8.0 1.0 2.032297 5.399466 1.001137 0.110963
.. ... ... ... ... ... ... ...
0 43.610746 9.0 1.0 1.323823 6.103425 1.704861 0.064712
0 43.611735 9.0 1.0 1.322435 6.104787 0.033876 0.064677
0 80.634129 12.0 2.0 1.321769 6.083292 0.000090 1.000175
0 80.634129 15.0 5.0 0.000000 5.225974 -1.000000 0.000100
0 13700.000000 15.0 13.0 0.000000 1.277584 -1.000000 0.000100
[107 rows x 7 columns]
Second, pick a certain resolution for the bcm array until the system mergers or is disrutped and then only print the final state
In [39]: bpp, bcm, initC, kick_info = Evolve.evolve(initialbinarytable=single_binary, BSEDict=BSEDict, timestep_conditions =[['binstate=0', 'dtp=1.0']])
In [40]: 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.381097 0
.. ... ... ... ... ... ...
0 78.0 12.0 1.0 1.321764 6.085939 0
0 79.0 12.0 1.0 1.321764 6.085078 0
0 80.0 12.0 1.0 1.321764 6.084196 0
0 81.0 15.0 13.0 0.000000 1.277584 1
0 13700.0 15.0 13.0 0.000000 1.277584 1
[83 rows x 6 columns]
Finally, we show how to print a fine resolution only during the HMXB stage of the evolution.
In [41]: 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 [42]: 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' : 0, 'zsun' : 0.014, 'bhms_coll_flag' : 0, 'don_lim' : -1, 'acc_lim' : -1, 'rtmsflag' : 0, 'wd_mass_lim': 1}
In [43]: 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 [44]: 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.071374 14.0 1.0 24.988585 106.891756 0
0 4.171374 14.0 1.0 24.988586 105.634889 0
0 4.271374 14.0 1.0 24.988586 104.335873 0
0 4.371374 14.0 1.0 24.988586 103.016130 0
0 4.471374 14.0 1.0 24.988586 101.707633 0
0 4.571374 14.0 1.0 24.988586 100.454871 0
0 4.671374 14.0 1.0 24.988586 99.314847 0
0 4.771374 14.0 1.0 24.988586 98.352198 0
0 4.871374 14.0 1.0 24.988590 89.040292 0
0 4.971374 14.0 7.0 24.990678 39.441122 0
0 5.071374 14.0 7.0 24.990681 35.988161 0
0 5.171374 14.0 7.0 24.990685 32.642016 0
0 5.271374 14.0 14.0 24.990687 31.054448 0
0 13700.000000 14.0 14.0 24.990687 31.054448 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
In [45]: 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)
In [46]: 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, 'remnantflag': 3, '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.5, 'ecsn_mlow' : 1.4, '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' : 0, 'bdecayfac' : 1, 'randomseed' : -1235453, 'grflag' : 1, 'rembar_massloss' : 0.5, 'kickflag' : 0, 'zsun' : 0.014, 'grflag' : 1, 'bhms_coll_flag' : 0, 'don_lim' : -1, 'acc_lim' : -1, 'rtmsflag' : 0, 'wd_mass_lim': 1}
In [47]: 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, BSEDict={})
....: if i == 3:
....: print("From beginning")
....: print(bpp)
....: print("Started in middle at Index {0}".format(i))
....: print(bpp_mid)
....:
From beginning
tphys mass_1 mass_2 ... bhspin_1 bhspin_2 bin_num
0 0.000000 25.543645 20.997840 ... 0.0 0.0 0
0 7.710210 24.912154 20.771614 ... 0.0 0.0 0
0 7.721745 24.907454 20.771066 ... 0.0 0.0 0
0 8.032373 24.637985 20.771096 ... 0.0 0.0 0
0 8.182493 8.914101 36.245340 ... 0.0 0.0 0
0 8.197721 8.892583 36.248877 ... 0.0 0.0 0
0 8.475201 8.580568 36.159241 ... 0.0 0.0 0
0 8.507067 8.527310 36.148466 ... 0.0 0.0 0
0 8.507067 8.027310 36.148466 ... 0.0 0.0 0
0 11.026338 8.027310 34.999378 ... 0.0 0.0 0
0 11.033157 8.027395 34.984228 ... 0.0 0.0 0
0 11.174781 8.045876 33.801486 ... 0.0 0.0 0
0 11.174781 8.045876 33.801486 ... 0.0 0.0 0
0 11.174781 8.045876 12.596798 ... 0.0 0.0 0
0 11.174781 8.045876 12.596798 ... 0.0 0.0 0
0 11.734073 8.073471 11.349314 ... 0.0 0.0 0
0 11.755369 8.075304 11.284171 ... 0.0 0.0 0
0 11.755369 8.075304 7.508719 ... 0.0 0.0 0
0 286.935007 8.075304 7.508719 ... 0.0 0.0 0
0 286.935007 0.000000 15.584023 ... 0.0 0.0 0
0 13700.000000 0.000000 15.584023 ... 0.0 0.0 0
[21 rows x 44 columns]
Started in middle at Index 3
tphys mass_1 mass_2 ... bhspin_1 bhspin_2 bin_num
0 8.032373 24.637985 20.771096 ... 0.0 0.0 0
0 8.182706 8.913234 36.246063 ... 0.0 0.0 0
0 8.197454 8.891835 36.248618 ... 0.0 0.0 0
0 8.475196 8.579641 36.158904 ... 0.0 0.0 0
0 8.507067 8.526389 36.148128 ... 0.0 0.0 0
0 8.507067 8.026389 36.148128 ... 0.0 0.0 0
0 11.026380 8.026389 34.999061 ... 0.0 0.0 0
0 11.033199 8.026474 34.983912 ... 0.0 0.0 0
0 11.174862 8.044952 33.800764 ... 0.0 0.0 0
0 11.174862 8.044952 33.800764 ... 0.0 0.0 0
0 11.174862 8.044952 12.596827 ... 0.0 0.0 0
0 11.174862 8.044952 12.596827 ... 0.0 0.0 0
0 11.734149 8.072531 11.349343 ... 0.0 0.0 0
0 11.755445 8.074363 11.284200 ... 0.0 0.0 0
0 11.755445 8.074363 7.508757 ... 0.0 0.0 0
0 287.246686 8.074363 7.508757 ... 0.0 0.0 0
0 287.246686 0.000000 15.583120 ... 0.0 0.0 0
0 13708.032373 0.000000 15.583120 ... 0.0 0.0 0
[18 rows x 44 columns]
Started in middle at Index 7
tphys mass_1 mass_2 ... bhspin_1 bhspin_2 bin_num
0 11.026338 8.027310 34.999378 ... 0.0 0.0 0
0 11.033157 8.027385 34.984228 ... 0.0 0.0 0
0 11.181821 8.046183 33.708746 ... 0.0 0.0 0
0 11.181821 8.046183 33.708746 ... 0.0 0.0 0
0 11.181821 8.046183 12.630280 ... 0.0 0.0 0
0 11.181821 8.046183 12.630280 ... 0.0 0.0 0
0 11.739568 8.070481 11.376838 ... 0.0 0.0 0
0 11.760721 8.072084 11.311761 ... 0.0 0.0 0
0 11.760721 8.072084 7.545724 ... 0.0 0.0 0
0 2459.988840 8.072084 7.545724 ... 0.0 0.0 0
0 2459.988840 0.000000 15.617808 ... 0.0 0.0 0
0 13708.507067 0.000000 15.617808 ... 0.0 0.0 0
[12 rows x 44 columns]
Started in middle at Index 11
tphys mass_1 mass_2 ... bhspin_1 bhspin_2 bin_num
0 11.174781 8.045876 33.801486 ... 0.0 0.0 0
0 11.174781 8.045876 33.801486 ... 0.0 0.0 0
0 11.174781 8.045876 12.596798 ... 0.0 0.0 0
0 11.174781 8.045876 12.596798 ... 0.0 0.0 0
0 11.734073 8.073471 11.349314 ... 0.0 0.0 0
0 11.755369 8.075304 11.284171 ... 0.0 0.0 0
0 11.755369 8.075304 7.508719 ... 0.0 0.0 0
0 286.935007 8.075304 7.508719 ... 0.0 0.0 0
0 286.935007 0.000000 15.584023 ... 0.0 0.0 0
0 13711.174781 0.000000 15.584023 ... 0.0 0.0 0
[10 rows x 44 columns]
One example of where restarting a binary can be extremely helpful is studying how natal kicks affect a binary independently of its previous evolution. This is particularly relevant for Gaia BH1 <https://ui.adsabs.harvard.edu/abs/2023MNRAS.518.1057E/abstract>_ and Gaia BH2 <https://ui.adsabs.harvard.edu/abs/2023MNRAS.521.4323E/abstract>_ which are difficult to produce through the standard common envelope channels. We can still study the effect of natal kicks on these binaries if we restart the evolution after the mass transfer would occur. We can do this by using a binary which gets us to the right masses given the metallicity, then overwrite some of the initial conditions to resample the natal kicks and pre-explosion separation.
In [48]: from cosmic import utils
....: import pandas as pd
....:
In [49]: single_binary = InitialBinaryTable.InitialBinaries(m1=65.0, m2=0.93, porb=4500, ecc=0.448872,
....: tphysf=13700.0, kstar1=1, kstar2=1, metallicity=0.014*0.6)
....:
In [50]: bpp, bcm, initC, kick_info = Evolve.evolve(initialbinarytable=single_binary, BSEDict=BSEDict)
In [51]: for column in bpp.columns:
....: initC = initC.assign(**{column:bpp.iloc[6][column]})
....:
In [52]: initC = pd.concat([initC]*1000)
....: initC['natal_kick_1'] = np.random.uniform(0, 100, 1000)
....: initC['phi_1'] = np.random.uniform(-90, 90, 1000)
....: initC['theta_1'] = np.random.uniform(0, 360, 1000)
....: initC['mean_anomaly_1'] = np.random.uniform(0, 360, 1000)
....: initC['porb'] = np.random.uniform(50, 190, 1000)
....: initC['sep'] = utils.a_from_p(p=initC.porb.values, m1=initC.mass_1.values, m2=initC.mass_2.values)
....: initC['bin_num'] = np.linspace(0, 1000, 1000)
....:
In [53]: bpp_restart, bcm_restart, initC_restart, kick_info_restart = Evolve.evolve(initialbinarytable=initC, BSEDict={})
In [54]: bpp_BH = bpp_restart.loc[(bpp_restart.kstar_1 == 14) & (bpp_restart.kstar_2 == 1) & (bpp_restart.porb > 0)].groupby('bin_num', as_index=False).first()
In [55]: bpp_BH[['tphys', 'mass_1', 'mass_2', 'porb', 'ecc']]
Out[55]:
tphys mass_1 mass_2 porb ecc
0 4.558013 9.416665 0.931311 1102.835361 0.885739
1 4.558013 9.416665 0.931311 201.584381 0.252218
2 4.558013 9.416665 0.931310 423.460024 0.208799
3 4.558013 9.416665 0.931310 851.260735 0.327119
4 4.558013 9.416665 0.931310 619.292939 0.381620
.. ... ... ... ... ...
577 4.558013 9.416665 0.931311 104.281790 0.919105
578 4.558013 9.416665 0.931311 974.628894 0.657606
579 4.558013 9.416665 0.931311 150.255960 0.396784
580 4.558013 9.416665 0.931311 314.560934 0.528562
581 4.558013 9.416665 0.931310 448.920994 0.116250
[582 rows x 5 columns]