Changing the way in which stars are sampled#

With cogsworth you can change all sorts of aspects of how you sample your initial stellar population, let’s talk about how this works.

Learning Goals#

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

  • Vary initial stellar sampling parameters using the sampling_params

  • Target a sample mass instead of a number of binaries

  • Include single stars in your population instead of just binaries

[1]:
import cogsworth
import numpy as np
import astropy.units as u
import matplotlib.pyplot as plt
[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)

Initial parameter distributions#

You can change the distribution of initial parameters using the sampling_params of a Population. The potential arguments are listed in the documentation for COSMIC.

Some to highlight are that you can change the distribution of primary masses, eccentricities and orbital periods with primary_model, ecc_model and porb_model respectively.

Let’s try sampling two different populations with different initial mass functions:

[3]:
p_kroupa = cogsworth.pop.Population(100000, sampling_params={"primary_model": "kroupa01"},
                                    use_default_BSE_settings=True)
p_kroupa.sample_initial_binaries()
[4]:
p_salpeter = cogsworth.pop.Population(100000, sampling_params={"primary_model": "salpeter55"},
                                      use_default_BSE_settings=True)
p_salpeter.sample_initial_binaries()
[5]:
fig, ax = plt.subplots()
for i, p in enumerate([p_kroupa, p_salpeter]):
    ax.hist(p._initial_binaries['mass_1'], bins=np.geomspace(0.08, 150, 100),
            alpha=1 - 0.5 * i, label=p.sampling_params["primary_model"]);
ax.set(xscale="log", yscale="log",
       xlabel=r"Initial mass $[\rm M_\odot]$", ylabel="Number of binaries")
ax.legend()
plt.show()
../../_images/tutorials_pop_settings_sampling_6_0.png

Sampling targets#

By default, cogsworth simulations target a number of binaries. However, it is possible to instead target a total sampled mass. This can be useful in a situation where, for example, you might want to sample a population for a cluster of a known mass from a hydrodynamical zoom-in simulation.

To do this you need to:

  • Pass sampling_target = "total_mass" in sampling params

  • Specify a total_mass in your sampling params

Let’s try an example:

[6]:
# n_binaries has to be passed so just pass None
p = cogsworth.pop.Population(n_binaries=None,
                             sampling_params={"sampling_target": "total_mass",
                                              "total_mass": 10000},
                             use_default_BSE_settings=True)
p.sample_initial_binaries()

Great, now we’ve got a population that resulted from sampling at least 10,000 solar masses. You can see this from checking mass_binaries and mass_singles.

[7]:
p.mass_binaries + p.mass_singles
[7]:
np.float64(17737.902908041917)

Note

Note that this is the total mass sampled, which is not the same as the total mass in the resulting sampled binaries (which don’t include single stars or the lowest mass stars). You could confirm this by summing up the mass_1 and mass_2 in the initC table and seeing that it’s smaller than the sum of mass_binaries + mass_singles

Targetting exact values#

You may have noticed that when specifying a target total mass that COSMIC overshoots the value you ask for. It is possible to ask cogsworth to throw away any extra samples to ensure you get as close as possible to the desired value (it’s not possible to get the exact value since each binary has some finite mass).

You do this by passing trim_extra_samples = True to your sampling params.

[8]:
p = cogsworth.pop.Population(n_binaries=None,
                             sampling_params={"sampling_target": "total_mass",
                                              "total_mass": 10000,
                                              "trim_extra_samples": True},  # <--- The new bit
                             use_default_BSE_settings=True)
p.sample_initial_binaries()

And now the total sampled mass should be nearly exactly 10,000

[9]:
p.mass_binaries + p.mass_singles
[9]:
np.float64(9997.57621266758)

Retaining single stars#

You’ll likely have seen that I’m rather biased towards binaries, but you need not be restricted to them! By default COSMIC throws away every sampled single star but you can keep by simply passing keep_singles = True in your sampling params

[10]:
p_with_singles = cogsworth.pop.Population(n_binaries=None,
                                          sampling_params={"sampling_target": "total_mass",
                                                           "total_mass": 10000,
                                                           "trim_extra_samples": True,
                                                           "keep_singles": True},  # <--- the new bit
                                          use_default_BSE_settings=True)
p_with_singles.sample_initial_binaries()

Now this population has single stars in its initial_binaries table. These still are written as binaries for consistency but the secondary stars are marked as massless remnants and the orbital period is exactly -1.0.

[11]:
p_with_singles._initial_binaries[p_with_singles._initial_binaries["porb"] == -1.0]
[11]:
index kstar_1 kstar_2 mass_1 mass_2 porb ecc metallicity tphysf mass0_1 ... tacc_1 tacc_2 epoch_1 epoch_2 tms_1 tms_2 bhspin_1 bhspin_2 tphys binfrac

0 rows × 39 columns

We can compare this to the earlier population which doesn’t have any single stars:

[12]:
(p._initial_binaries["porb"] > 0.0).all()
[12]:
np.True_
[13]:
(p_with_singles._initial_binaries["porb"] > 0.0).all()
[13]:
np.True_

Wrap-up#

And that’s all for this one - hopefully you’ve got a better idea of how to change the way in which you sample your initial stellar population. Be sure to read on to learn about how to change other settings in your simulations.

Note

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