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]: SSEDict = {'stellar_engine': 'sse'}
In [5]: 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, "LBV_flag": 1, "alpha1": [1.0, 1.0],
...: "lambdaf": 0.0, "ce2stageflag": 0, "ceflag": 0, "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,
...: "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,
...: "fryer_mass_limit": 0, "mxns": 3.0, "fryer_fmix": 1.0,
...: "fryer_mcrit_nsbh": 5.75, "rembar_massloss": 0.5, "wd_mass_lim": 1,
...: "maltsev_mode": 0, "maltsev_fallback": 0.5, "maltsev_pf_prob": 0.1,
...: "pisn": -2, "ppi_co_shift": 0.0, "ppi_extra_ml": 0.0, "bhspinflag": 0,
...: "bhspinmag": 0.0, "grflag": 1, "eddfac": 10, "gamma": -2, "don_lim": -1,
...: "acc_lim": [-1, -1], "smt_periastron_check": 0, "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 [6]: bpp, bcm, initC, kick_info = Evolve.evolve(initialbinarytable=single_binary, BSEDict=BSEDict, SSEDict=SSEDict)
In [7]: print("From beginning")
From beginning
In [8]: 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.209233 20.702742 ... 0.0 0.0 0
0 8.698975 18.945631 20.795256 ... 0.0 0.0 0
0 8.698975 10.604460 20.795256 ... 0.0 0.0 0
0 9.948208 10.604460 20.728562 ... 0.0 0.0 0
0 9.964123 10.604460 20.725065 ... 0.0 0.0 0
0 10.918726 10.604460 17.776137 ... 0.0 0.0 0
0 10.936001 10.604460 17.592566 ... 0.0 0.0 0
0 10.936001 10.604460 7.408961 ... 0.0 0.0 0
0 13700.000000 10.604460 7.408961 ... 0.0 0.0 0
[13 rows x 46 columns]
# assign row indices and find SNe
In [9]: bpp["row_iloc"] = list(range(len(bpp)))
In [10]: first_SN_iloc = bpp[bpp["evol_type"] == 15].iloc[0]["row_iloc"]
In [11]: second_SN_iloc = bpp[bpp["evol_type"] == 16].iloc[0]["row_iloc"]
In [12]: print("First SN Index: ", first_SN_iloc)
First SN Index: 5.0
In [13]: print("Second SN Index: ", second_SN_iloc)
Second SN Index: 10.0
# choose some inds to restart from
In [14]: before_SN_1 = int(first_SN_iloc // 2)
In [15]: between_SNe = int((first_SN_iloc + second_SN_iloc) // 2)
In [16]: after_SN_2 = int((second_SN_iloc + bpp["row_iloc"].max()) // 2)
In [17]: print("Restart Indices: ", before_SN_1, between_SNe, after_SN_2)
Restart Indices: 2 7 11
# restart from different points
In [18]: 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.209233 20.702742 ... 0.0 0.0 0
0 8.698975 18.945631 20.795256 ... 0.0 0.0 0
0 8.698975 10.604460 20.795256 ... 0.0 0.0 0
0 9.948208 10.604460 20.728562 ... 0.0 0.0 0
0 9.964123 10.604460 20.725065 ... 0.0 0.0 0
0 10.918726 10.604460 17.776137 ... 0.0 0.0 0
0 10.936001 10.604460 17.592566 ... 0.0 0.0 0
0 10.936001 10.604460 7.408961 ... 0.0 0.0 0
0 13707.886566 10.604460 7.408961 ... 0.0 0.0 0
[10 rows x 46 columns]
Started in middle at Index 7
tphys mass_1 mass_2 ... bhspin_1 bhspin_2 bin_num
0 9.964123 10.60446 20.725065 ... 0.0 0.0 0
0 10.918726 10.60446 17.776137 ... 0.0 0.0 0
0 10.936001 10.60446 17.592566 ... 0.0 0.0 0
0 10.936001 10.60446 7.408961 ... 0.0 0.0 0
0 13709.948208 10.60446 7.408961 ... 0.0 0.0 0
[5 rows x 46 columns]
Started in middle at Index 11
tphys mass_1 mass_2 ... bhspin_1 bhspin_2 bin_num
0 13710.936001 10.60446 7.408961 ... 0.0 0.0 0
0 13710.936001 10.60446 7.408961 ... 0.0 0.0 0
[2 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 [19]: from cosmic import utils
In [20]: import pandas as pd
In [21]: 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 [22]: bpp, bcm, initC, kick_info = Evolve.evolve(
....: initialbinarytable=single_binary, BSEDict=BSEDict, SSEDict=SSEDict
....: )
....:
In [23]: for column in bpp.columns:
....: initC = initC.assign(**{column:bpp.iloc[6][column]})
....:
In [24]: initC = pd.concat([initC]*1000)
In [25]: initC['natal_kick_1'] = np.random.uniform(0, 100, 1000)
In [26]: initC['phi_1'] = np.random.uniform(-90, 90, 1000)
In [27]: initC['theta_1'] = np.random.uniform(0, 360, 1000)
In [28]: initC['mean_anomaly_1'] = np.random.uniform(0, 360, 1000)
In [29]: initC['porb'] = np.random.uniform(50, 190, 1000)
In [30]: initC['sep'] = utils.a_from_p(p=initC.porb.values, m1=initC.mass_1.values, m2=initC.mass_2.values)
In [31]: initC['bin_num'] = np.linspace(0, 1000, 1000)
In [32]: bpp_restart, bcm_restart, initC_restart, kick_info_restart = Evolve.evolve(
....: initialbinarytable=initC, BSEDict=BSEDict, SSEDict=SSEDict
....: )
....:
In [33]: 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 [34]: bpp_BH[['tphys', 'mass_1', 'mass_2', 'porb', 'ecc']]
Out[34]:
Empty DataFrame
Columns: [tphys, mass_1, mass_2, porb, ecc]
Index: []