Creating your first population#

Welcome to the cogsworth tutorials! In this tutorial we will create a simple population of binaries.

Learning Goals#

By the end of this tutorial you should know:

[15]:

import cogsworth
import gala.potential as gp
import astropy.units as u
import matplotlib.pyplot as plt
import numpy as np

[3]:

# 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)

The Population Class#

The core functionality of cogsworth is to simulate populations of binaries in the context of the galaxy. This is implemented primarily through the Population class, which will give you access to the vast majority of cogsworth’s features.

Let’s initialise a simple population of binaries using the default input settings:

[5]:

p = cogsworth.pop.Population(n_binaries=100, processes=4,
                             final_kstar1=list(range(16)), final_kstar2=list(range(16)),
                             sfh_model=cogsworth.sfh.Wagg2022,
                             galactic_potential=gp.MilkyWayPotential(version='v2'),
                             v_dispersion=5 * u.km / u.s,
                             max_ev_time=12.0 * u.Gyr, timestep_size=1 * u.Myr,
                             BSE_settings={}, sampling_params={}, store_entire_orbits=True,
                             use_default_BSE_settings=True)

Basic input settings#

The Population class has a number of input settings that can be used to control the way in which the population is created. Many of these directly affect the way in which a population is evolved or the initial conditions are sampled and these are covered in later tutorials.

For now, we will focus on the basic settings:

  • n_binaries: By default, cogsworth samples a population by targetting a number of binaries, this is that target value

  • processes: This is the number of CPUs to use for the population, which uses a multiprocessing pool for a variety of tasks. In general, the speed of the simulation scales with the number of processors.

  • max_ev_time: How long to evolve a simulation for

The rest of the settings will be addressed in those other tutorials. Try running the next part with different values for these parameters and seeing what happens

Creating a population#

The manual way#

First let’s look a how you create a population manually, step-by-step (which can be useful for custom populations later on!).

[6]:

p = cogsworth.pop.Population(n_binaries=1000, processes=4, max_ev_time=12.0 * u.Gyr,
                             use_default_BSE_settings=True)

Initial stellar population#

Let’s start by sampling initial stellar population. This will create the _initial_binaries table as well as store values for the mass_singles, mass_binaries, n_singles_req, n_bin_req, which track the mass and number of single stars and binary stars necessary to sample this population.

Check out the sampling tutorial for more info and changing the sampling settings.

[7]:

# now let's get the initital stellar population
p.sample_initial_binaries()
p._initial_binaries

[7]:
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 0 0.0 0.0 0.157566 0.132221 197194.064203 0.622412 0.013086 5682.087153 0.157566 ... 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0
1 1 0.0 0.0 0.660025 0.220672 47.162822 0.027374 0.010722 6579.674348 0.660025 ... 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0
2 2 1.0 1.0 11.815256 8.041803 8.478898 0.019202 0.019260 6643.735592 11.815256 ... 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0
3 3 0.0 0.0 0.472437 0.202251 2110.567622 0.143982 0.013821 6660.107450 0.472437 ... 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0
4 4 1.0 1.0 1.405933 1.368043 36.503566 0.351799 0.011648 3456.999292 1.405933 ... 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0
... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ...
1303 1303 0.0 0.0 0.321245 0.315190 1.667328 0.097813 0.007657 10601.159078 0.321245 ... 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0
1304 1304 0.0 0.0 0.282806 0.087128 18.389173 0.522408 0.012857 8281.116937 0.282806 ... 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0
1305 1305 0.0 0.0 0.241308 0.154405 2.889133 0.587927 0.016882 8470.409857 0.241308 ... 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0
1306 1306 1.0 0.0 1.099550 0.618352 12160.670995 0.027828 0.008546 9869.178612 1.099550 ... 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0
1307 1307 1.0 0.0 0.815015 0.493761 69.382317 0.654700 0.009505 7792.279021 0.815015 ... 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0

1308 rows × 39 columns

[8]:

fig, ax = plt.subplots()
ax.scatter(p._initial_binaries["mass_1"], p._initial_binaries["mass_2"])
ax.set(xscale="log", yscale="log",
       xlabel=r"Initial primary mass $[\rm M_\odot]$", ylabel=r"Initial secondary mass $[\rm M_\odot]$")
plt.show()
../../_images/tutorials_basics_pops_10_0.png

Initial galaxy sampling#

Now let’s try sampling the initial galactic positions, kinematics and birth times. This will create the initial_galaxy attribute, which is a Galaxy that was instantiated based on your input choice of sfh_model when creating the Population.

Check out the star formation history tutorial for more info and changing initial galaxy sampling settings.

[9]:

# first let's sample the initial galaxy
p.sample_initial_galaxy()

# and check distributions look normal
plt.hist(p.initial_galaxy.rho.value, bins="fd");
plt.xlabel(r"Galactocentric $\rho \, [\rm kpc]$");
../../_images/tutorials_basics_pops_12_0.png

Stellar evolution#

Now we can use COSMIC to evolve the stellar population from birth until present day. This will create a series of output tables: initC (the initial conditions for each binary), bpp (the evolution table each binary), bcm (the conditioned timestep evolution table for each binary) and kick_info (the kicks experienced by each binary).

Check out this tutorial for more info and changing the population synthesis settings.

[10]:

p.perform_stellar_evolution()
p.bpp

[10]:
tphys mass_1 mass_2 kstar_1 kstar_2 sep porb ecc RRLO_1 RRLO_2 ... B_2 bacc_1 bacc_2 tacc_1 tacc_2 epoch_1 epoch_2 bhspin_1 bhspin_2 bin_num
0 0.000000 0.157566 0.132221 0 0 9431.345261 197194.064203 0.622412 0.000138 0.000130 ... 0.0 0.0 0.0 0.0 0.0 0.000000 0.000000 0.0 0.0 0
0 5682.087153 0.157566 0.132221 0 0 9431.345261 197194.064203 0.622412 0.000138 0.000130 ... 0.0 0.0 0.0 0.0 0.0 0.000000 0.000000 0.0 0.0 0
1 0.000000 0.660025 0.220672 0 0 52.637351 47.162822 0.027374 0.025274 0.016495 ... 0.0 0.0 0.0 0.0 0.0 0.000000 0.000000 0.0 0.0 1
1 6579.674348 0.660025 0.220672 0 0 52.637334 47.162798 0.027374 0.025789 0.016525 ... 0.0 0.0 0.0 0.0 0.0 0.000000 0.000000 0.0 0.0 1
2 0.000000 11.815256 8.041803 1 1 47.367793 8.478898 0.019202 0.238627 0.226264 ... 0.0 0.0 0.0 0.0 0.0 0.000000 0.000000 0.0 0.0 2
... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ...
1306 7792.922709 0.849999 0.619220 6 0 3106.837676 16557.958413 0.027791 0.152780 0.000553 ... 0.0 0.0 0.0 0.0 0.0 -6207.630424 37.576069 0.0 0.0 1306
1306 7793.902178 0.460000 0.624964 11 0 3794.605436 26008.446229 0.027525 0.000011 0.000396 ... 0.0 0.0 0.0 0.0 0.0 7793.902178 284.499601 0.0 0.0 1306
1306 9869.178612 0.460000 0.624964 11 0 5102.149166 40550.311918 0.027525 0.000008 0.000296 ... 0.0 0.0 0.0 0.0 0.0 7793.902178 284.499601 0.0 0.0 1306
1307 0.000000 0.815015 0.493761 1 0 77.697338 69.382317 0.654700 0.065259 0.050086 ... 0.0 0.0 0.0 0.0 0.0 0.000000 0.000000 0.0 0.0 1307
1307 7792.279021 0.815015 0.493761 1 0 77.695127 69.379355 0.654690 0.070742 0.050747 ... 0.0 0.0 0.0 0.0 0.0 0.000000 0.000000 0.0 0.0 1307

4139 rows × 46 columns

[11]:

# print out some evolution information for the most massive binary
# we'll cover more on how this masking works in the next few tutorials
p.bpp.loc[p.initial_binaries[p.initial_binaries["mass_1"] == p.initial_binaries["mass_1"].max()]["bin_num"].iloc[0]]

[11]:
tphys mass_1 mass_2 kstar_1 kstar_2 sep porb ecc RRLO_1 RRLO_2 ... B_2 bacc_1 bacc_2 tacc_1 tacc_2 epoch_1 epoch_2 bhspin_1 bhspin_2 bin_num
966 0.000000 21.290014 9.558130 1 1 6363.572493 10592.760098 0.330790 0.002544 0.002257 ... 0.000000e+00 0.0 0.0 0.0 0.0 0.000000 0.000000 0.0 0.0 966
966 9.339993 20.982598 9.556822 2 1 6427.900606 10808.003948 0.330790 0.006110 0.002608 ... 0.000000e+00 0.0 0.0 0.0 0.0 -0.094812 -0.001113 0.0 0.0 966
966 9.355676 20.979662 9.556820 4 1 6428.518958 10810.083668 0.330790 0.023743 0.002608 ... 0.000000e+00 0.0 0.0 0.0 0.0 -0.094812 -0.001117 0.0 0.0 966
966 10.295175 18.631822 9.643188 5 1 5828.143893 9697.642571 0.069340 0.660799 0.002052 ... 0.000000e+00 0.0 0.0 0.0 0.0 -0.094812 0.148776 0.0 0.0 966
966 10.312262 18.476582 9.652698 5 1 5751.590690 9531.797508 0.049206 0.660799 0.002052 ... 0.000000e+00 0.0 0.0 0.0 0.0 -0.094812 0.165256 0.0 0.0 966
966 10.312262 8.032066 9.652698 14 1 -1.000000 -1.000000 -1.000000 0.000000 -2.000000 ... 0.000000e+00 0.0 0.0 0.0 0.0 10.312262 0.165256 0.0 0.0 966
966 27.860784 8.032066 9.634363 14 2 -1.000000 -1.000000 -1.000000 0.000100 0.000100 ... 0.000000e+00 0.0 0.0 0.0 0.0 10.312262 0.089712 0.0 0.0 966
966 27.942244 8.032066 9.632528 14 3 -1.000000 -1.000000 -1.000000 0.000100 0.000100 ... 0.000000e+00 0.0 0.0 0.0 0.0 10.312262 0.080960 0.0 0.0 966
966 27.947231 8.032066 9.631883 14 4 -1.000000 -1.000000 -1.000000 0.000100 0.000100 ... 0.000000e+00 0.0 0.0 0.0 0.0 10.312262 0.080960 0.0 0.0 966
966 30.845812 8.032066 9.406922 14 5 -1.000000 -1.000000 -1.000000 0.000100 0.000100 ... 0.000000e+00 0.0 0.0 0.0 0.0 10.312262 0.080960 0.0 0.0 966
966 30.951651 8.032066 9.298136 14 5 -1.000000 -1.000000 -1.000000 0.000100 0.000100 ... 0.000000e+00 0.0 0.0 0.0 0.0 10.312262 0.080960 0.0 0.0 966
966 30.951651 8.032066 1.277584 14 13 -1.000000 -1.000000 -1.000000 0.000100 0.000100 ... 0.000000e+00 0.0 0.0 0.0 0.0 10.312262 30.951651 0.0 0.0 966
966 11947.433362 8.032066 1.277584 14 13 -1.000000 -1.000000 -1.000000 0.000100 0.000100 ... 4.737887e+10 0.0 0.0 0.0 0.0 10.312262 30.951651 0.0 0.0 966

13 rows × 46 columns

Galactic evolution#

Finally, let’s bring Gala into the mix and evolve this population of binaries through a galactic potential. This will create orbits for each binary in orbits that can be accessed easily as primary_orbits and secondary_orbits.

Check out the potentials tutorial for more info and changing the galactic evolution settings (in particular the potential).

[12]:

p.perform_galactic_evolution()
p.orbits

Integrating orbits: 100%|██████████| 1314/1314 [00:02<00:00, 556.00it/s]

[12]:

array([<Orbit cartesian, dim=3, shape=(7771,)>,
       <Orbit cartesian, dim=3, shape=(762,)>,
       <Orbit cartesian, dim=3, shape=(6339,)>, ...,
       <Orbit cartesian, dim=3, shape=(10750,)>,
       <Orbit cartesian, dim=3, shape=(11112,)>,
       <Orbit cartesian, dim=3, shape=(10325,)>],
      shape=(1314,), dtype=object)

[16]:

fig, ax = plt.subplots()
p.orbits[np.random.randint(0, len(p.orbits))].cylindrical.plot(["rho", "z"], axes=[ax]);
../../_images/tutorials_basics_pops_18_0.png

The easy way#

Now that was useful to understand what’s happening behind the scenes, but for many cases you just want to efficiently create a population. This can be done all at once with the create_population().

[19]:

p.create_population()

Run for 1000 binaries
Sampled 1348 binaries
[1e-02s] Sample initial binaries
[0.6s] Evolve binaries (run COSMIC)

Integrating orbits: 100%|██████████| 1350/1350 [00:02<00:00, 584.73it/s]

[4.4s] Integrate galactic orbits (run gala)
Overall: 5.0s

Check out the next tutorial to learn more about interpretting the outputs of these simulations.

Wrap-up#

And that’s all for this first tutorial of the basic series, nice job! We covered a lot here and much of the functionality I used will get explained in more detail in the later tutorials. Press your right arrow key for the next installment or head over to the main tutorials page for more options.

Note

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