Restarting a binary¶
COSMIC allows you to restart a binary from any point in its evolution from a COSMIC generated bpp array.
In [1]: from cosmic.sample.initialbinarytable import InitialBinaryTable
In [2]: from cosmic.evolve import Evolve
Below we provide an example of the same evolutionary track started from the beginning and three different points in the evolution:
sometime between the beginning and the first object going supernova
between the first and second supernova
after both supernova
In [3]: single_binary = InitialBinaryTable.InitialBinaries(
...: m1=25, m2=20, porb=6000, ecc=0.0,
...: tphysf=13700.0, kstar1=1, kstar2=1, metallicity=0.002
...: )
...:
In [4]: BSEDict = {
...: "pts1": 0.001, "pts2": 0.01, "pts3": 0.02, "zsun": 0.014, "windflag": 3,
...: "eddlimflag": 0, "neta": 0.5, "bwind": 0.0, "hewind": 0.5, "beta": 0.125,
...: "xi": 0.5, "acc2": 1.5, "alpha1": 1.0, "lambdaf": 0.0, "ceflag": 1,
...: "cekickflag": 2, "cemergeflag": 1, "cehestarflag": 0, "qcflag": 5,
...: "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],
...: "kickflag": 5, "sigma": 265.0, "bhflag": 1, "bhsigmafrac": 1.0,
...: "sigmadiv": -20.0, "ecsn": 2.25, "ecsn_mlow": 1.6, "aic": 1, "ussn": 1,
...: "pisn": -2, "polar_kick_angle": 90.0,
...: "natal_kick_array": [[-100.0, -100.0, -100.0, -100.0, 0.0], [-100.0, -100.0, -100.0, -100.0, 0.0]],
...: "mm_mu_ns": 400.0, "mm_mu_bh": 200.0, "remnantflag": 4, "mxns": 3.0,
...: "rembar_massloss": 0.5, "wd_mass_lim": 1, "maltsev_mode": 0,
...: "maltsev_fallback": 0.5, "maltsev_pf_prob": 0.1, "bhspinflag": 0,
...: "bhspinmag": 0.0, "grflag": 1, "eddfac": 10, "gamma": -2, "don_lim": -1,
...: "acc_lim": -1, "tflag": 1, "ST_tide": 1,
...: "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],
...: "ifflag": 1, "wdflag": 1, "epsnov": 0.001, "bdecayfac": 1,
...: "bconst": 3000, "ck": 1000, "rejuv_fac": 1.0, "rejuvflag": 0,
...: "bhms_coll_flag": 0, "htpmb": 1, "ST_cr": 1, "rtmsflag": 0
...: }
...:
# evolve the binary
In [5]: bpp, bcm, initC, kick_info = Evolve.evolve(initialbinarytable=single_binary, BSEDict=BSEDict)
In [6]: print("From beginning")
From beginning
In [7]: print(bpp)
tphys mass_1 mass_2 ... bhspin_1 bhspin_2 bin_num
0 0.000000 25.000000 20.000000 ... 0.0 0.0 0
0 7.874629 24.411688 19.823083 ... 0.0 0.0 0
0 7.886566 24.407128 19.822598 ... 0.0 0.0 0
0 8.419117 22.865321 19.931095 ... 0.0 0.0 0
0 8.684362 19.208942 20.702707 ... 0.0 0.0 0
0 8.698975 18.945363 20.795179 ... 0.0 0.0 0
0 8.698975 10.604316 20.795179 ... 0.0 0.0 0
0 8.698975 10.604316 20.795179 ... 0.0 0.0 0
0 9.948142 10.604324 20.728490 ... 0.0 0.0 0
0 9.964057 10.604327 20.724993 ... 0.0 0.0 0
0 10.271866 10.604525 20.637428 ... 0.0 0.0 0
0 10.436765 10.991291 19.245419 ... 0.0 0.0 0
0 10.436765 10.991291 19.245419 ... 0.0 0.0 0
0 10.436765 10.991291 6.844259 ... 0.0 0.0 0
0 10.436765 10.991291 6.844259 ... 0.0 0.0 0
0 10.879530 11.036427 6.619769 ... 0.0 0.0 0
0 10.926937 11.046294 6.576965 ... 0.0 0.0 0
0 10.926937 11.046294 2.363262 ... 0.0 0.0 0
0 503.719302 11.046294 2.363262 ... 0.0 0.0 0
0 503.719302 0.000000 13.409556 ... 0.0 0.0 0
0 503.719302 0.000000 13.409556 ... 0.0 0.0 0
0 13700.000000 0.000000 13.409556 ... 0.0 0.0 0
[22 rows x 46 columns]
# assign row indices and find SNe
In [8]: bpp["row_iloc"] = list(range(len(bpp)))
In [9]: first_SN_iloc = bpp[bpp["evol_type"] == 15].iloc[0]["row_iloc"]
In [10]: second_SN_iloc = bpp[bpp["evol_type"] == 16].iloc[0]["row_iloc"]
In [11]: print("First SN Index: ", first_SN_iloc)
First SN Index: 5.0
In [12]: print("Second SN Index: ", second_SN_iloc)
Second SN Index: 16.0
# choose some inds to restart from
In [13]: before_SN_1 = int(first_SN_iloc // 2)
In [14]: between_SNe = int((first_SN_iloc + second_SN_iloc) // 2)
In [15]: after_SN_2 = int((second_SN_iloc + bpp["row_iloc"].max()) // 2)
In [16]: print("Restart Indices: ", before_SN_1, between_SNe, after_SN_2)
Restart Indices: 2 10 18
# restart from different points
In [17]: for i in [before_SN_1, between_SNe, after_SN_2]:
....: new_initC = initC.copy()
....: for column in bpp.columns:
....: new_initC = new_initC.assign(**{column:bpp.iloc[i][column]})
....: bpp_mid, bcm_mid, initC_mid, kick_info = Evolve.evolve(initialbinarytable=new_initC, BSEDict={})
....: print("Started in middle at Index {0}".format(i))
....: print(bpp_mid)
....:
Started in middle at Index 2
tphys mass_1 mass_2 ... bhspin_1 bhspin_2 bin_num
0 8.419117 22.865321 19.931095 ... 0.0 0.0 0
0 8.684362 19.208942 20.702707 ... 0.0 0.0 0
0 8.698975 18.945363 20.795179 ... 0.0 0.0 0
0 8.698975 10.604316 20.795179 ... 0.0 0.0 0
0 8.698975 10.604316 20.795179 ... 0.0 0.0 0
0 9.948142 10.604324 20.728490 ... 0.0 0.0 0
0 9.964057 10.604327 20.724993 ... 0.0 0.0 0
0 10.271866 10.604525 20.637428 ... 0.0 0.0 0
0 10.436765 10.991291 19.245419 ... 0.0 0.0 0
0 10.436765 10.991291 19.245419 ... 0.0 0.0 0
0 10.436765 10.991291 6.844259 ... 0.0 0.0 0
0 10.436765 10.991291 6.844259 ... 0.0 0.0 0
0 10.879530 11.036427 6.619769 ... 0.0 0.0 0
0 10.926937 11.046294 6.576965 ... 0.0 0.0 0
0 10.926937 11.046294 2.363262 ... 0.0 0.0 0
0 503.727510 11.046294 2.363262 ... 0.0 0.0 0
0 503.727510 0.000000 13.409556 ... 0.0 0.0 0
0 503.727510 0.000000 13.409556 ... 0.0 0.0 0
0 13707.886566 0.000000 13.409556 ... 0.0 0.0 0
[19 rows x 46 columns]
Started in middle at Index 10
tphys mass_1 mass_2 ... bhspin_1 bhspin_2 bin_num
0 10.271866 10.604525 20.637428 ... 0.0 0.0 0
0 10.435873 10.985748 19.289632 ... 0.0 0.0 0
0 10.435873 10.985748 19.289632 ... 0.0 0.0 0
0 10.435873 10.985748 6.842286 ... 0.0 0.0 0
0 10.435873 10.985748 6.842286 ... 0.0 0.0 0
0 10.879633 11.031398 6.617598 ... 0.0 0.0 0
0 10.927067 11.041367 6.574806 ... 0.0 0.0 0
0 10.927067 11.041367 2.362235 ... 0.0 0.0 0
0 484.437881 11.041367 2.362235 ... 0.0 0.0 0
0 484.437881 0.000000 13.403602 ... 0.0 0.0 0
0 484.437881 0.000000 13.403602 ... 0.0 0.0 0
0 13710.271866 0.000000 13.403602 ... 0.0 0.0 0
[12 rows x 46 columns]
Started in middle at Index 18
tphys mass_1 mass_2 ... bhspin_1 bhspin_2 bin_num
0 503.719302 11.046294 2.363262 ... 0.0 0.0 0
0 503.719302 0.000000 13.409556 ... 0.0 0.0 0
0 503.719302 0.000000 13.409556 ... 0.0 0.0 0
0 14114.819302 0.000000 13.409556 ... 0.0 0.0 0
[4 rows x 46 columns]
Natal kick example¶
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 and Gaia BH2 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 [18]: from cosmic import utils
In [19]: import pandas as pd
In [20]: 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 [21]: bpp, bcm, initC, kick_info = Evolve.evolve(
....: initialbinarytable=single_binary, BSEDict=BSEDict
....: )
....:
In [22]: for column in bpp.columns:
....: initC = initC.assign(**{column:bpp.iloc[6][column]})
....:
In [23]: initC = pd.concat([initC]*1000)
In [24]: initC['natal_kick_1'] = np.random.uniform(0, 100, 1000)
In [25]: initC['phi_1'] = np.random.uniform(-90, 90, 1000)
In [26]: initC['theta_1'] = np.random.uniform(0, 360, 1000)
In [27]: initC['mean_anomaly_1'] = np.random.uniform(0, 360, 1000)
In [28]: initC['porb'] = np.random.uniform(50, 190, 1000)
In [29]: initC['sep'] = utils.a_from_p(p=initC.porb.values, m1=initC.mass_1.values, m2=initC.mass_2.values)
In [30]: initC['bin_num'] = np.linspace(0, 1000, 1000)
In [31]: bpp_restart, bcm_restart, initC_restart, kick_info_restart = Evolve.evolve(
....: initialbinarytable=initC, BSEDict={}
....: )
....:
In [32]: 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 [33]: bpp_BH[['tphys', 'mass_1', 'mass_2', 'porb', 'ecc']]
Out[33]:
tphys mass_1 mass_2 porb ecc
0 4.557585 9.701219 0.931296 136.567702 0.827225
1 4.557585 9.701219 0.931294 200.886971 0.950187
2 4.557585 9.701219 0.931303 350.883911 0.400289
3 4.557585 9.701219 0.931302 1273.689122 0.738289
4 4.557585 9.701219 0.931295 1161.514175 0.554208
.. ... ... ... ... ...
586 4.557585 9.701219 0.931301 335.800870 0.323143
587 4.557585 9.701219 0.931295 669.273282 0.508196
588 4.557585 9.701219 0.931296 412.078441 0.368801
589 4.557585 9.701219 0.931295 536.634944 0.176515
590 4.557585 9.701219 0.931298 356.751625 0.230096
[591 rows x 5 columns]