Example use case: Common-envelope evolution
In this notebook we look at how common-envelope evolution (CEE) alters binary-star orbits. We construct a population of low- and intermediate-mass binaries and compare their orbital periods before and after CEE. Not all stars evolve into this phase, so we have to run a whole population to find those that do. We then have to construct the pre- and post-CEE distributions and plot them.
First, we import a few required Python modules.
[1]:
import os
import math
import matplotlib.pyplot as plt
from binarycpython.utils.functions import temp_dir
from binarycpython.utils.grid import Population
TMP_DIR = temp_dir("notebooks", "notebook_comenv")
Setting up the Population object
We set up a new population object. Our stars evolve to \(13.7\text{ }\mathrm{Gyr}\), the age of the Universe, and we assume the metallicity \(Z=0.02\). We also set the common-envelope ejection efficiency \(\alpha_\mathrm{CE}=1\) and the envelope structure parameter \(\lambda=0.5\). More complex options are available in binary_c, such as \(\lambda\) based on stellar mass, but this is just a demonstration example so let’s keep things simple.
[2]:
# Create population object
population = Population()
population.set(
# grid options
tmp_dir = TMP_DIR,
verbosity = 1,
log_dt = 10, # log every 10 seconds
# binary-star evolution options
max_evolution_time=13700, # maximum stellar evolution time in Myr (13700 Myr == 13.7 Gyr)
metallicity=0.02, # 0.02 is approximately Solar metallicity
alpha_ce = 1.0,
lambda_ce = 0.5,
)
adding: log_dt=10 to grid_options
adding: max_evolution_time=13700 to BSE_options
adding: metallicity=0.02 to BSE_options
adding: alpha_ce=1.0 to BSE_options
adding: lambda_ce=0.5 to BSE_options
Stellar Grid
We now construct a grid of stars, varying the mass from \(1\) to \(6\text{ }\mathrm{M}_\odot\). We avoid massive stars for now, and focus on the (more common) low- and intermediate-mass stars. We also limit the period range to \(10^4\text{ }\mathrm{d}\) because systems with longer orbital periods will probably not undergo Roche-lobe overflow and hence common-envelope evolution is impossible.
[3]:
import binarycpython.utils.distribution_functions
# Set resolution and mass range that we simulate
resolution = {"M_1": 10, "q" : 10, "per": 10}
massrange = [1, 6]
logperrange = [0.15, 4]
population.add_grid_variable(
name="lnm1",
longname="Primary mass",
valuerange=massrange,
resolution="{}".format(resolution["M_1"]),
spacingfunc="const(math.log({min}), math.log({max}), {res})".format(min=massrange[0],max=massrange[1],res=resolution["M_1"]),
precode="M_1=math.exp(lnm1)",
probdist="three_part_powerlaw(M_1, 0.1, 0.5, 1.0, 150, -1.3, -2.3, -2.3)*M_1",
dphasevol="dlnm1",
parameter_name="M_1",
condition="", # Impose a condition on this grid variable. Mostly for a check for yourself
)
# Mass ratio
population.add_grid_variable(
name="q",
longname="Mass ratio",
valuerange=["0.1/M_1", 1],
resolution="{}".format(resolution['q']),
spacingfunc="const({}/M_1, 1, {})".format(massrange[0],resolution['q']),
probdist="flatsections(q, [{{'min': {}/M_1, 'max': 1.0, 'height': 1}}])".format(massrange[0]),
dphasevol="dq",
precode="M_2 = q * M_1",
parameter_name="M_2",
condition="", # Impose a condition on this grid variable. Mostly for a check for yourself
)
# Orbital period
population.add_grid_variable(
name="log10per", # in days
longname="log10(Orbital_Period)",
valuerange=[0.15, 5.5],
resolution="{}".format(resolution["per"]),
spacingfunc="const({}, {}, {})".format(logperrange[0],logperrange[1],resolution["per"]),
precode="""orbital_period = 10.0 ** log10per
sep = calc_sep_from_period(M_1, M_2, orbital_period)
sep_min = calc_sep_from_period(M_1, M_2, 10**{})
sep_max = calc_sep_from_period(M_1, M_2, 10**{})""".format(logperrange[0],logperrange[1]),
probdist="sana12(M_1, M_2, sep, orbital_period, sep_min, sep_max, math.log10(10**{}), math.log10(10**{}), {})".format(logperrange[0],logperrange[1],-0.55),
parameter_name="orbital_period",
dphasevol="dlog10per",
)
Added grid variable: {
"name": "lnm1",
"longname": "Primary mass",
"valuerange": [
1,
6
],
"resolution": "10",
"spacingfunc": "const(math.log(1), math.log(6), 10)",
"precode": "M_1=math.exp(lnm1)",
"probdist": "three_part_powerlaw(M_1, 0.1, 0.5, 1.0, 150, -1.3, -2.3, -2.3)*M_1",
"dphasevol": "dlnm1",
"parameter_name": "M_1",
"condition": "",
"gridtype": "centred",
"branchpoint": 0,
"grid_variable_number": 0
}
Added grid variable: {
"name": "q",
"longname": "Mass ratio",
"valuerange": [
"0.1/M_1",
1
],
"resolution": "10",
"spacingfunc": "const(1/M_1, 1, 10)",
"precode": "M_2 = q * M_1",
"probdist": "flatsections(q, [{'min': 1/M_1, 'max': 1.0, 'height': 1}])",
"dphasevol": "dq",
"parameter_name": "M_2",
"condition": "",
"gridtype": "centred",
"branchpoint": 0,
"grid_variable_number": 1
}
Added grid variable: {
"name": "log10per",
"longname": "log10(Orbital_Period)",
"valuerange": [
0.15,
5.5
],
"resolution": "10",
"spacingfunc": "const(0.15, 4, 10)",
"precode": "orbital_period = 10.0 ** log10per\nsep = calc_sep_from_period(M_1, M_2, orbital_period)\nsep_min = calc_sep_from_period(M_1, M_2, 10**0.15)\nsep_max = calc_sep_from_period(M_1, M_2, 10**4)",
"probdist": "sana12(M_1, M_2, sep, orbital_period, sep_min, sep_max, math.log10(10**0.15), math.log10(10**4), -0.55)",
"dphasevol": "dlog10per",
"parameter_name": "orbital_period",
"condition": null,
"gridtype": "centred",
"branchpoint": 0,
"grid_variable_number": 2
}
Logging and handling the output
We now construct the pre- and post-common envelope evolution data for the first common envelope that forms in each binary. We look at the comenv_count variable, we can see that when it increases from 0 to 1 we have found our object. If this happens, we stop evolution of the system to save CPU time.
[4]:
custom_logging_statement = """
/*
* Detect when the comenv_count increased
*/
if(stardata->model.comenv_count == 1 &&
stardata->previous_stardata->model.comenv_count == 0)
{
/*
* We just had this system's first common envelope:
* output the time at which this happens,
* the system's probability (proportional to the number of stars),
* the previous timestep's (pre-comenv) orbital period (days) and
* the current timestep (post-comenv) orbital period (days)
*/
Printf("COMENV %g %g %g %g\\n",
stardata->model.time,
stardata->model.probability,
stardata->previous_stardata->common.orbit.period * YEAR_LENGTH_IN_DAYS,
stardata->common.orbit.period * YEAR_LENGTH_IN_DAYS);
/*
* We should waste no more CPU time on this system now we have the
* data we want.
*/
stardata->model.evolution_stop = TRUE;
}
"""
population.set(
C_logging_code=custom_logging_statement
)
adding: C_logging_code=
/*
* Detect when the comenv_count increased
*/
if(stardata->model.comenv_count == 1 &&
stardata->previous_stardata->model.comenv_count == 0)
{
/*
* We just had this system's first common envelope:
* output the time at which this happens,
* the system's probability (proportional to the number of stars),
* the previous timestep's (pre-comenv) orbital period (days) and
* the current timestep (post-comenv) orbital period (days)
*/
Printf("COMENV %g %g %g %g\n",
stardata->model.time,
stardata->model.probability,
stardata->previous_stardata->common.orbit.period * YEAR_LENGTH_IN_DAYS,
stardata->common.orbit.period * YEAR_LENGTH_IN_DAYS);
/*
* We should waste no more CPU time on this system now we have the
* data we want.
*/
stardata->model.evolution_stop = TRUE;
}
to grid_options
The parse function must now catch lines that start with “COMENV” and process the associated data. We set up the parse_data function to do just this.
[5]:
from binarycpython.utils.functions import bin_data,datalinedict
import re
# log-period distribution bin width (dex)
binwidth = 0.5
def parse_function(self, output):
"""
Parsing function to convert HRD data into something that Python can use
"""
# list of the data items
parameters = ["header", "time", "probability", "pre_comenv_period", "post_comenv_period"]
# Loop over the output.
for line in output.splitlines():
# obtain the line of data in dictionary form
linedata = datalinedict(line,parameters)
# choose COMENV lines of output
if linedata["header"] == "COMENV":
# bin the pre- and post-comenv log10-orbital-periods to nearest 0.5dex
binned_pre_period = bin_data(math.log10(linedata["pre_comenv_period"]), binwidth)
# but check if the post-comenv period is finite and positive: if
# not, the system has merged and we give it an aritifical period
# of 10^-100 days (which is very much unphysical)
if linedata["post_comenv_period"] > 0.0:
binned_post_period = bin_data(math.log10(linedata["post_comenv_period"]), binwidth)
else:
binned_post_period = bin_data(-100,binwidth) # merged!
# make the "histograms"
self.grid_results['pre'][binned_pre_period] += linedata["probability"]
self.grid_results['post'][binned_post_period] += linedata["probability"]
# verbose reporting
#print("parse out results_dictionary=",self.grid_results)
# Add the parsing function
population.set(
parse_function=parse_function,
)
adding: parse_function=<function parse_function at 0x14736bebc040> to grid_options
Evolving the grid
Now we actually run the population. This may take a little while. You can set num_cores higher if you have a powerful machine.
[6]:
# set number of threads
population.set(
# set number of threads (i.e. number of CPU cores we use)
num_cores=4,
)
# Evolve the population - this is the slow, number-crunching step
analytics = population.evolve()
# Show the results (debugging)
#print (population.grid_results)
adding: num_cores=4 to grid_options
Creating and loading custom logging functionality
Generating grid code
Generating grid code
Constructing/adding: lnm1
Constructing/adding: q
Constructing/adding: log10per
Saving grid code to grid_options
Writing grid code to /tmp/binary_c_python/notebooks/notebook_comenv/binary_c_grid_ad303100d719457c83256568f9a9887c.py
Loading grid code function from /tmp/binary_c_python/notebooks/notebook_comenv/binary_c_grid_ad303100d719457c83256568f9a9887c.py
Grid code loaded
Grid has handled 1000 stars
with a total probability of 0.0645905996773004
Total starcount for this run will be: 1000
[2021-09-12 18:07:39,950 DEBUG Process-2] --- Setting up processor: process-0
[2021-09-12 18:07:39,953 DEBUG Process-3] --- Setting up processor: process-1
[2021-09-12 18:07:39,959 DEBUG Process-4] --- Setting up processor: process-2
[2021-09-12 18:07:39,962 DEBUG MainProcess] --- setting up the system_queue_filler now
[2021-09-12 18:07:39,965 DEBUG Process-5] --- Setting up processor: process-3
Process 0 started at 2021-09-12T18:07:39.965721. Using store memaddr <capsule object "STORE" at 0x14736bee47e0>
Process 1 started at 2021-09-12T18:07:39.970949. Using store memaddr <capsule object "STORE" at 0x14736bee4870>
Process 2 started at 2021-09-12T18:07:39.978355. Using store memaddr <capsule object "STORE" at 0x14736bee4f30>
Process 3 started at 2021-09-12T18:07:39.983689. Using store memaddr <capsule object "STORE" at 0x14736bee4870>
[2021-09-12 18:07:40,066 DEBUG MainProcess] --- Signaling stop to processes
Generating grid code
Generating grid code
Constructing/adding: lnm1
Constructing/adding: q
Constructing/adding: log10per
Saving grid code to grid_options
Writing grid code to /tmp/binary_c_python/notebooks/notebook_comenv/binary_c_grid_ad303100d719457c83256568f9a9887c.py
Loading grid code function from /tmp/binary_c_python/notebooks/notebook_comenv/binary_c_grid_ad303100d719457c83256568f9a9887c.py
Grid code loaded
163/1000 16.3% complete 18:07:49 ETA= 51.5s tpr=6.16e-02 ETF=18:08:41 mem:594.9MB
322/1000 32.2% complete 18:07:59 ETA= 42.9s tpr=6.33e-02 ETF=18:08:42 mem:538.2MB
465/1000 46.5% complete 18:08:09 ETA= 38.1s tpr=7.12e-02 ETF=18:08:47 mem:538.2MB
586/1000 58.6% complete 18:08:19 ETA= 34.3s tpr=8.29e-02 ETF=18:08:54 mem:540.0MB
682/1000 68.2% complete 18:08:30 ETA= 34.0s tpr=1.07e-01 ETF=18:09:04 mem:540.1MB
784/1000 78.4% complete 18:08:40 ETA= 21.2s tpr=9.81e-02 ETF=18:09:01 mem:541.8MB
872/1000 87.2% complete 18:08:50 ETA= 15.0s tpr=1.17e-01 ETF=18:09:05 mem:546.1MB
963/1000 96.3% complete 18:09:00 ETA= 4.2s tpr=1.14e-01 ETF=18:09:04 mem:546.9MB
[2021-09-12 18:09:06,366 DEBUG Process-5] --- Process-3 is finishing.
Process 3 finished:
generator started at 2021-09-12T18:07:39.964604, done at 2021-09-12T18:09:06.370832 (total: 86.406228s of which 86.24177551269531s interfacing with binary_c).
Ran 222 systems with a total probability of 0.014137215791516371.
This thread had 0 failing systems with a total probability of 0.
Skipped a total of 0 systems because they had 0 probability
[2021-09-12 18:09:06,374 DEBUG Process-5] --- Process-3 is finished.
[2021-09-12 18:09:06,979 DEBUG Process-3] --- Process-1 is finishing.
Process 1 finished:
generator started at 2021-09-12T18:07:39.953039, done at 2021-09-12T18:09:06.982866 (total: 87.029827s of which 86.82909393310547s interfacing with binary_c).
Ran 273 systems with a total probability of 0.01877334232598154.
This thread had 0 failing systems with a total probability of 0.
Skipped a total of 0 systems because they had 0 probability
[2021-09-12 18:09:06,985 DEBUG Process-3] --- Process-1 is finished.
[2021-09-12 18:09:07,174 DEBUG Process-2] --- Process-0 is finishing.
Process 0 finished:
generator started at 2021-09-12T18:07:39.949775, done at 2021-09-12T18:09:07.176660 (total: 87.226885s of which 87.02672934532166s interfacing with binary_c).
Ran 268 systems with a total probability of 0.016469813170514686.
This thread had 0 failing systems with a total probability of 0.
Skipped a total of 0 systems because they had 0 probability
[2021-09-12 18:09:07,179 DEBUG Process-2] --- Process-0 is finished.
[2021-09-12 18:09:07,233 DEBUG Process-4] --- Process-2 is finishing.
Process 2 finished:
generator started at 2021-09-12T18:07:39.958802, done at 2021-09-12T18:09:07.236252 (total: 87.27745s of which 87.0905077457428s interfacing with binary_c).
Ran 237 systems with a total probability of 0.015210228389288167.
This thread had 0 failing systems with a total probability of 0.
Skipped a total of 0 systems because they had 0 probability
[2021-09-12 18:09:07,238 DEBUG Process-4] --- Process-2 is finished.
Population-ad303100d719457c83256568f9a9887c finished! The total probability was: 0.06459059967730076. It took a total of 87.54819011688232s to run 1000 systems on 4 cores
There were no errors found in this run.
After the run is complete, some technical report on the run is returned. I stored that in analytics
. As we can see below, this dictionary is like a status report of the evolution. Useful for e.g. debugging. We check this, and then set about making the plot of the orbital period distributions using Seaborn.
[7]:
print(analytics)
{'population_name': 'ad303100d719457c83256568f9a9887c', 'evolution_type': 'grid', 'failed_count': 0, 'failed_prob': 0, 'failed_systems_error_codes': [], 'errors_exceeded': False, 'errors_found': False, 'total_probability': 0.06459059967730076, 'total_count': 1000, 'start_timestamp': 1631462859.9342952, 'end_timestamp': 1631462947.4824853, 'total_mass_run': 4680.235689312421, 'total_probability_weighted_mass_run': 0.22611318083528567, 'zero_prob_stars_skipped': 0}
[8]:
# make a plot of the distributions
import seaborn as sns
import pandas as pd
import copy
pd.set_option("display.max_rows", None, "display.max_columns", None)
from binarycpython.utils.functions import pad_output_distribution
# set up seaborn for use in the notebook
sns.set(rc={'figure.figsize':(20,10)})
sns.set_context("notebook",
font_scale=1.5,
rc={"lines.linewidth":2.5})
pd.set_option("display.max_rows", None, "display.max_columns", None)
# remove the merged objects
probability = { "merged" : 0.0, "unmerged" : 0.0}
# copy the results so we can change the copy
results = copy.deepcopy(population.grid_results)
for distribution in ['post']:
for logper in population.grid_results[distribution]:
dprob = results[distribution][logper]
if logper < -90:
# merged system
probability["merged"] += dprob
del results[distribution][logper]
else:
# unmerged system
probability["unmerged"] += dprob
print(probability)
# pad the final distribution with zero
for distribution in population.grid_results:
pad_output_distribution(results[distribution],
binwidth)
# make pandas dataframe
plot_data = pd.DataFrame.from_dict(results, orient='columns')
# make the plot
p = sns.lineplot(data=plot_data)
p.set_xlabel("$\log_{10} (P_\mathrm{orb} / \mathrm{day})$")
p.set_ylabel("Number of stars")
#p.set(xlim=(-5,5)) # might be necessary?
{'merged': 0.035263029200000025, 'unmerged': 0.019388724199999995}
[8]:
Text(0, 0.5, 'Number of stars')

You can see that common-envelope evolution shrinks stellar orbits, just as we expect. Pre-CEE, most orbits are in the range \(10\) to \(1000\text{ }\mathrm{d}\), while after CEE the distribution peaks at about \(1\text{ }\mathrm{d}\). Some of these orbits are very short: \(\log_{10}(-2) = 0.01\text{ }\mathrm{d}\sim10\text{ }\mathrm{minutes}\). Such systems are prime candidates for exciting astrophysics: novae, type Ia supernovae and gravitational wave sources.
Things to try: * Extend the logging to output more data than just the orbital period. * What are the stellar types of the post-common envelope systems? Are they likely to undergo novae or a type-Ia supernova? * What are the lifetimes of the systems in close (\(<1\text{ }\mathrm{d}\)) binaries? Are they likely to merge in the life of the Universe? * How much mass is lost in common-envelope interactions? * Extend the grid to massive stars. Do you see many NS and BH compact binaries? * Try different \(\alpha_\mathrm{CE}\) and \(\lambda_\mathrm{CE}\) options… * … and perhaps increased resolution to obtain smoother curves. * Why do long-period systems not reach common envelope evolution?