Varying stellar evolution parameters#

cogsworth offers all of the same flexibility as COSMIC for changing stellar evolution parameters. In this tutorials we’ll demonstrate how to change them.

Learning Goals#

By the end of this tutorial you should know how to:

Vary stellar evolution parameters using the BSE_settings

[1]:

import cogsworth
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd

[2]:

# this all just makes plots look nice
%config InlineBackend.figure_format = 'retina'

plt.rc('font', family='serif')
plt.rcParams['text.usetex'] = False
fs = 24

# update various fontsizes to match
params = {'figure.figsize': (12, 8),
          'legend.fontsize': 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}
plt.rcParams.update(params)
pd.options.display.max_columns = 999

The `BSE_settings`#

All population synthesis settings are changed using the BSE_settings input to the Population (BSE here stands for binary stellar evolution). This dictionary contains the default COSMIC settings but any of these settings can be overwritten by the user.

Supernova kick example#

Let’s try an example where we change the strength of supernova natal kicks for two populations (which I’ve skewed to target systems that are likely to go supernova)

[3]:

p = cogsworth.pop.Population(n_binaries=100,
                             final_kstar1=[13, 14], final_kstar2=[13, 14],
                             BSE_settings={"sigma": 30},
                             use_default_BSE_settings=True)
p.sample_initial_binaries()
p.perform_stellar_evolution()

First let’s see what fraction of the population were disrupted when drawing kicks from a Maxwellian centred on $30 \, {\rm km / s}$

[4]:

# count the number of binaries that got kicks while they hadn't already merged
n_kicked = len(p.bpp[((p.bpp["evol_type"] == 15) | (p.bpp["evol_type"] == 16)) & (p.bpp["sep"] != 0.0)]["bin_num"].unique())
f'{p.disrupted.sum() / n_kicked:1.2f}'

[4]:

'0.92'

Now let’s change the Maxwellian to be centred on $1000 \, {\rm km / s}$ and do the same for electron-capture supernovae.

[5]:

# we can drop the BSE_settings and just change the initial_binaries settings
p.BSE_settings = {}
p.initial_binaries["sigma"] = 1000
p.initial_binaries["sigmadiv"] = -1000
p.perform_stellar_evolution()

[6]:

# count the number of binaries that got kicks while they hadn't already merged
n_kicked = len(p.bpp[((p.bpp["evol_type"] == 15) | (p.bpp["evol_type"] == 16)) & (p.bpp["sep"] != 0.0)]["bin_num"].unique())
f'{p.disrupted.sum() / n_kicked:1.2f}'

[6]:

'0.92'

As you’d expect, the fraction is now much higher (and any that aren’t disrupting are likely getting fallback adjusted kicks, got lucky directions, or had extremely high orbital velocities).

Tip

Notice that I changed the setting without needing to redefine a new population in this last example (I just updated p.BSE_settings). This also means that I evolved the same initial population with different physics. This can be a useful way to ensure that a change in your results is directly due to different physics rather than some sort of stochastic sampling effect.

Common-envelope example#

Let’s do another example in which we change the efficiency of the common-envelope phase (alpha1), which is defined as the fraction of orbital energy that goes into unbinding the envelope.

[7]:

# create a population of 100 binaries
p = cogsworth.pop.Population(n_binaries=1000, use_default_BSE_settings=True)
p.sample_initial_binaries()
p.perform_stellar_evolution()

# for a series of different alphas
for alpha in [10, 2, 1, 0.5, 0.1]:
    # change the setting and perform evolution
    p.BSE_settings = {}
    p.initial_binaries["alpha1"] = alpha
    p.perform_stellar_evolution()

    # find all binaries that start a CE and how many end as a mergers
    ce_nums = p.bpp[p.bpp["evol_type"] == 7]["bin_num"].unique()
    n_merged = (p.final_bpp.loc[ce_nums]["sep"] == 0.0).sum()

    # get separations after the first CE for each binary
    seps = p.bpp[p.bpp["evol_type"] == 8].drop_duplicates(subset="bin_num", keep="first")["sep"].values

    print(f"Run for alpha = {alpha}")
    print(f"  {len(ce_nums)} binaries experience at least one common-envelope event")
    print(f"  {n_merged / len(ce_nums) * 100:1.1f}% of these resulted in a stellar merger")
    print(f"  For surviving systems, mean separation immediately after the CE is {np.mean(seps):1.2f} Rsun\n")

Run for alpha = 10
  117 binaries experience at least one common-envelope event
  50.4% of these resulted in a stellar merger
  For surviving systems, mean separation immediately after the CE is 24.59 Rsun

Run for alpha = 2
  117 binaries experience at least one common-envelope event
  74.4% of these resulted in a stellar merger
  For surviving systems, mean separation immediately after the CE is 9.92 Rsun

Run for alpha = 1
  117 binaries experience at least one common-envelope event
  79.5% of these resulted in a stellar merger
  For surviving systems, mean separation immediately after the CE is 5.68 Rsun

Run for alpha = 0.5
  117 binaries experience at least one common-envelope event
  85.5% of these resulted in a stellar merger
  For surviving systems, mean separation immediately after the CE is 2.88 Rsun

Run for alpha = 0.1
  117 binaries experience at least one common-envelope event
  92.3% of these resulted in a stellar merger
  For surviving systems, mean separation immediately after the CE is 0.28 Rsun

So what have we learned? Well first of all, the number of binaries that experiences a common-envelope is unaffected by the efficiency - that’s definitely what we expect (and we should be worried if it changed!).

We see two trends which are that the percentage of systems that merge during the common-envelope increases with lower efficiencies, whilst the mean separation of the surviving systems decreases. Why do you think this is?

And the reason for these trends is…

$\alpha_{\rm CE}$ defines the fraction of the orbital energy that goes into unbinding the envelope. So a lower $\alpha_{\rm CE}$ means that that a binary must inspiral significantly closer before it can achieve the same unbinding energy (and end the common-envelope phase). In many cases, this results in the binary inspiraling too far and merging. Fun!

Definitions of all settings#

All of the many varied potential population synthesis settings are listed in great detail on the COSMIC documentation. You can check them out here and find out how to change your favourite bit of physics!

Wrap-up#

And that brings us to the end of this short but sweet tutorial. Be sure to look at the COSMIC docs linked above to find out about all of the parameters (and importantly, their default values!). Read on to learn about changing other settings like the galactic potential or galactic star formation model.

Note

This tutorial was generated from a Jupyter notebook that can be found here.