Modifying internal timesteps¶
In addition to modifying the resolution of the output data, users can also modify the internal timesteps used by COSMIC.
Changing the default timesteps (not recommended)¶
By default, COSMIC timesteps are set by the user-specified settings pts1, pts2, and pts3 (see Configuration files for more details). These specify timesteps for different evolutionary phases of the binary.
One could edit these directly in the BSEDict as we’ve seen in other tutorials.
In [1]: from cosmic.sample.initialbinarytable import InitialBinaryTable
In [2]: from cosmic.evolve import Evolve
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]: 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, "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,
...: "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
...: }
...:
In [6]: bpp, bcm, initC, kick_info = Evolve.evolve(
...: initialbinarytable=single_binary,
...: BSEDict=BSEDict,
...: SSEDict=SSEDict
...: )
...:
# take twice as large timesteps for all phases
In [7]: BSEDict['pts1'] *= 2
In [8]: BSEDict['pts2'] *= 2
In [9]: BSEDict['pts3'] *= 2
In [10]: bpp_mod, bcm_mod, initC_mod, kick_info_mod = Evolve.evolve(
....: initialbinarytable=single_binary,
....: BSEDict=BSEDict,
....: SSEDict=SSEDict
....: )
....:
From inspection of the resulting tables, you may start to see strange behaviour if the timesteps are too large.
Modifying timesteps as a function of mass (recommended)¶
Users can also specify additional modifiers to these timesteps that will be applied during the evolution. The original BSE timesteps can be modified for higher masses, where extrapolation beyond the original stellar grids can lead to unexpected behaviour.
These modifiers are specified as dt_mass_modifiers in the evolve() function when performing evolution. They are written as a list of tuples of the form (min, max, modifier), where the fractional modifier is applied to the default timesteps for stars with masses between min and max.
For example, to take half as large timesteps for stars with masses between 50 and 100 solar masses, you can do:
In [11]: dt_mass_modifiers = [(50, 100, 0.5)]
In [12]: bpp_mod, bcm_mod, initC_mod, kick_info_mod = Evolve.evolve(
....: initialbinarytable=single_binary,
....: BSEDict=BSEDict,
....: SSEDict=SSEDict,
....: dt_mass_modifiers=dt_mass_modifiers
....: )
....:
The default choice in COSMIC is to set dt_mass_modifiers = [(40, 70, 0.3), (70, np.inf, 0.1)], which we find leads to more numerically stable evolution for very massive stars.
Example: BH masses¶
Let’s go through an example to see how these modifiers can affect the results of our simulations. Specifically, let’s explore the relation between initial mass and BH mass for a high-mass grid.
First, we can create a grid of 500 single stars with initial masses between 50 and 150 solar masses.
In [13]: import matplotlib.pyplot as plt
In [14]: import numpy as np
In [15]: N_BINARIES = 500
# create an initial binary table of single stars
In [16]: IBT_bh = InitialBinaryTable.InitialBinaries(
....: m1=np.linspace(50, 150, N_BINARIES),
....: m2=np.zeros(N_BINARIES),
....: porb=np.zeros(N_BINARIES),
....: ecc=np.zeros(N_BINARIES),
....: kstar1=[1] * N_BINARIES,
....: kstar2=[0] * N_BINARIES,
....: metallicity=[0.01] * N_BINARIES,
....: tphysf=[200] * N_BINARIES,
....: )
....:
Now let’s copy our BSEDict but turn of pair-instability supernovae (PISN) to simplify the top end of the plot and make sure any differences are just coming from timesteps.
# switch off PISN to simplify the plot
In [17]: no_PISN_BSEDict = BSEDict.copy()
In [18]: no_PISN_BSEDict["pisn"] = 0
We can set up 4 different choices of timesteps to compare:
The default choice in
COSMIC:dt_mass_modifiers = [(40, 70, 0.3), (70, np.inf, 0.1)]Half the default BSE timesteps:
dt_mass_modifiers = [(0, np.inf, 0.5)]The default BSE timesteps:
dt_mass_modifiers = []Double the default BSE timesteps:
dt_mass_modifiers = [(0, np.inf, 2)]
In [19]: bpp_timesteps = []
In [20]: labels = ['Default dt_mass_modifiers', 'Half BSE',
....: 'BSE defaults', 'Double BSE']
....:
In [21]: dt_mass_modifiers_vals = [
....: [(40, 70, 0.3), (70, np.inf, 0.1)],
....: [(0, np.inf, 0.5)],
....: [],
....: [(0, np.inf, 2)]
....: ]
....:
Now let’s use those to evolve our grid of single stars 4 different times and store the resulting BPPs.
In [22]: for dtmm in dt_mass_modifiers_vals:
....: bpp, _, _, _ = Evolve.evolve(
....: initialbinarytable=IBT_bh,
....: SSEDict=SSEDict,
....: BSEDict=BSEDict,
....: dt_mass_modifiers=dtmm,
....: )
....: bpp_timesteps.append(bpp)
....:
With our evolution done, let’s find the final BH mass and the corresponding initial mass in each grid for plotting.
In [23]: m_inits = []
In [24]: m_bhs = []
In [25]: for bpp in bpp_timesteps:
....: bh = bpp.loc[bpp['kstar_1'] == 14].groupby(level=0).last()
....: init = IBT_bh.loc[bh.index, 'mass_1']
....: m_init, m_bh = init.values, bh['mass_1'].values
....: m_inits.append(m_init)
....: m_bhs.append(m_bh)
....:
And finally, let’s plot the results!
In [26]: plt.rc('font', family='serif')
In [27]: plt.rcParams['text.usetex'] = False
In [28]: fs = 24
In [29]: params = {'figure.figsize': (12, 8),
....: 'legend.fontsize': 0.7*fs,
....: 'legend.title_fontsize': 0.8*fs,
....: 'axes.labelsize': fs,
....: 'xtick.labelsize': 0.9 * fs,
....: 'ytick.labelsize': 0.9 * fs,
....: 'axes.linewidth': 1.1,
....: 'xtick.major.size': 7,
....: 'xtick.minor.size': 4,
....: 'ytick.major.size': 7,
....: 'ytick.minor.size': 4}
....:
In [30]: plt.rcParams.update(params)
In [31]: fig, axes = plt.subplots(2, 1, figsize=(14, 8), gridspec_kw={'height_ratios': [3, 1]})
In [32]: for m_init, m_bh, label, zorder in zip(m_inits, m_bhs, labels, [100, 50, 25, 1]):
....: axes[0].plot(m_init, m_bh, zorder=zorder, label=label)
....:
In [33]: axes[0].legend();
In [34]: for m_init, m_bh in zip(m_inits, m_bhs):
....: axes[1].plot(m_init, m_bh - m_bhs[0])
....:
In [35]: axes[1].set_xlabel('Initial primary mass $M_1^\\mathrm{init}$ [$M_\\odot$]');
In [36]: for ax in axes:
....: ax.grid(True, lw=0.4, alpha=0.4);
....:
In [37]: axes[0].set_ylabel('BH mass [$M_\\odot$]');
In [38]: axes[1].set_ylabel('Difference [$M_\\odot$]');
In [39]: plt.tight_layout();
In [40]: plt.show()
So we can see here that there are clear numerical artifacts in the relation between initial mass and BH mass when using the default BSE timesteps (and especially when they are larger!). This motivated our choice of the default dt_mass_modifiers in COSMIC to ensure more stable evolution for very massive stars.