# 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,
        samplerfunc="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],
         samplerfunc="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],
        samplerfunc="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",
        "parameter_name": "M_1",
        "longname": "Primary mass",
        "valuerange": [
            1,
            6
        ],
        "samplerfunc": "const(math.log(1), math.log(6), 10)",
        "precode": "M_1=math.exp(lnm1)",
        "postcode": null,
        "probdist": "three_part_powerlaw(M_1, 0.1, 0.5, 1.0, 150, -1.3, -2.3, -2.3)*M_1",
        "dphasevol": "dlnm1",
        "condition": "",
        "gridtype": "centred",
        "branchpoint": 0,
        "branchcode": null,
        "topcode": null,
        "bottomcode": null,
        "grid_variable_number": 0
    }
    Added grid variable: {
        "name": "q",
        "parameter_name": "M_2",
        "longname": "Mass ratio",
        "valuerange": [
            "0.1/M_1",
            1
        ],
        "samplerfunc": "const(1/M_1, 1, 10)",
        "precode": "M_2 = q * M_1",
        "postcode": null,
        "probdist": "flatsections(q, [{'min': 1/M_1, 'max': 1.0, 'height': 1}])",
        "dphasevol": "dq",
        "condition": "",
        "gridtype": "centred",
        "branchpoint": 0,
        "branchcode": null,
        "topcode": null,
        "bottomcode": null,
        "grid_variable_number": 1
    }
    Added grid variable: {
        "name": "log10per",
        "parameter_name": "orbital_period",
        "longname": "log10(Orbital_Period)",
        "valuerange": [
            0.15,
            5.5
        ],
        "samplerfunc": "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)",
        "postcode": null,
        "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",
        "condition": null,
        "gridtype": "centred",
        "branchpoint": 0,
        "branchcode": null,
        "topcode": null,
        "bottomcode": null,
        "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 0x149c95c56c10> 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
    Doing dry run to calculate total starcount and probability
    Generating grid code
    Generating grid code
    Saving grid code to grid_options
    Writing grid code to /tmp/binary_c_python-izzard/notebooks/notebook_comenv/binary_c_grid_2b66f805db424c48a1d29c45092b6e3c.py [dry_run = True]
    Symlinked grid code to /tmp/binary_c_python-izzard/notebooks/notebook_comenv/binary_c_grid-latest0
    Loading grid code function from /tmp/binary_c_python-izzard/notebooks/notebook_comenv/binary_c_grid_2b66f805db424c48a1d29c45092b6e3c.py
    Grid code loaded
    Dry run of the grid
    Grid has handled 1000 stars with a total probability of 0.0645906
    ****************************************
    * Total starcount for this run is 1000 *
    *    Total probability is 0.0645906    *
    ****************************************

    [2021-11-01 09:50:56,156 DEBUG    Process-2] --- Setting up processor: process-0
    [2021-11-01 09:50:56,159 DEBUG    Process-3] --- Setting up processor: process-1
    [2021-11-01 09:50:56,163 DEBUG    Process-4] --- Setting up processor: process-2
    [2021-11-01 09:50:56,166 DEBUG    MainProcess] --- setting up the system_queue_filler now

    Generating grid code
    Generating grid code
    Saving grid code to grid_options

    [2021-11-01 09:50:56,168 DEBUG    Process-5] --- Setting up processor: process-3

    Writing grid code to /tmp/binary_c_python-izzard/notebooks/notebook_comenv/binary_c_grid_2b66f805db424c48a1d29c45092b6e3c.py [dry_run = False]
    Symlinked grid code to /tmp/binary_c_python-izzard/notebooks/notebook_comenv/binary_c_grid-latest1
    Loading grid code function from /tmp/binary_c_python-izzard/notebooks/notebook_comenv/binary_c_grid_2b66f805db424c48a1d29c45092b6e3c.py
    Grid code loaded

    [2021-11-01 09:50:56,258 DEBUG    MainProcess] --- Signaling stop to processes

    177/1000  17.7% complete 09:51:05 ETA=   46.9s tpr=2.28e-01 ETF=09:51:52 mem:941.1MB M1=1.3 M2=1.2 P=1.1e+3
    348/1000  34.8% complete 09:51:15 ETA=   38.1s tpr=2.34e-01 ETF=09:51:53 mem:679.5MB M1=1.9 M2=1.4 P=2.6e+3
    451/1000  45.1% complete 09:51:25 ETA=   32.1s tpr=2.34e-01 ETF=09:51:57 mem:682.6MB M1=2.2 M2=1.7 P=5.3
    555/1000  55.5% complete 09:51:35 ETA=   32.8s tpr=2.95e-01 ETF=09:52:08 mem:684.4MB M1=2.7 M2=1.9 P=1.9e+2
    631/1000  63.1% complete 09:51:46 ETA=   29.7s tpr=3.22e-01 ETF=09:52:15 mem:685.3MB M1=3.2 M2=1.8 P=5.3
    682/1000  68.2% complete 09:51:56 ETA=   28.4s tpr=3.58e-01 ETF=09:52:24 mem:687.5MB M1=3.2 M2=2.9 P=13
    757/1000  75.7% complete 09:52:06 ETA=   24.4s tpr=4.02e-01 ETF=09:52:30 mem:687.9MB M1=3.8 M2=2.6 P=1.1e+3
    810/1000  81.0% complete 09:52:16 ETA=   19.9s tpr=4.20e-01 ETF=09:52:36 mem:688.0MB M1=4.6 M2=1.5 P=2.2
    863/1000  86.3% complete 09:52:27 ETA=   15.3s tpr=4.48e-01 ETF=09:52:42 mem:691.9MB M1=4.6 M2=3.3 P=31
    915/1000  91.5% complete 09:52:37 ETA=   10.1s tpr=4.76e-01 ETF=09:52:47 mem:693.2MB M1=5.5 M2=1.7 P=1.9e+2
    963/1000  96.3% complete 09:52:47 ETA=    4.6s tpr=4.97e-01 ETF=09:52:52 mem:694.2MB M1=5.5 M2=3.9 P=31

    [2021-11-01 09:52:56,316 DEBUG    Process-5] --- Process-3 is finishing.

    process 3 free memory and return
    ****************************************************
    *                Process 3 finished:               *
    *  generator started at 2021-11-01T09:50:56.168379 *
    * generator finished at 2021-11-01T09:52:56.323359 *
    *                  total: 2m 0.15s                 *
    *          of which 2m 0.03s with binary_c         *
    *                  Ran 251 systems                 *
    *       with a total probability of 0.0163656      *
    *         This thread had 0 failing systems        *
    *       with a total failed probability of 0       *
    *   Skipped a total of 0 zero-probability systems  *
    *                                                  *
    ****************************************************

    process 3 queue put output_dict

    [2021-11-01 09:52:56,329 DEBUG    Process-5] --- Process-3 is finished.

    process 3 return

    [2021-11-01 09:52:56,986 DEBUG    Process-2] --- Process-0 is finishing.

    process 0 free memory and return
    ****************************************************
    *                Process 0 finished:               *
    *  generator started at 2021-11-01T09:50:56.155678 *
    * generator finished at 2021-11-01T09:52:56.991657 *
    *                  total: 2m 0.84s                 *
    *          of which 2m 0.70s with binary_c         *
    *                  Ran 267 systems                 *
    *       with a total probability of 0.0175264      *
    *         This thread had 0 failing systems        *
    *       with a total failed probability of 0       *
    *   Skipped a total of 0 zero-probability systems  *
    *                                                  *
    ****************************************************

    process 0 queue put output_dict

    [2021-11-01 09:52:56,996 DEBUG    Process-2] --- Process-0 is finished.

    process 0 return

    [2021-11-01 09:52:57,094 DEBUG    Process-3] --- Process-1 is finishing.

    process 1 free memory and return
    ****************************************************
    *                Process 1 finished:               *
    *  generator started at 2021-11-01T09:50:56.158640 *
    * generator finished at 2021-11-01T09:52:57.099417 *
    *                  total: 2m 0.94s                 *
    *          of which 2m 0.83s with binary_c         *
    *                  Ran 234 systems                 *
    *       with a total probability of 0.0143896      *
    *         This thread had 0 failing systems        *
    *       with a total failed probability of 0       *
    *   Skipped a total of 0 zero-probability systems  *
    *                                                  *
    ****************************************************

    process 1 queue put output_dict

    [2021-11-01 09:52:57,104 DEBUG    Process-3] --- Process-1 is finished.

    process 1 return

    [2021-11-01 09:52:57,728 DEBUG    Process-4] --- Process-2 is finishing.

    process 2 free memory and return
    ****************************************************
    *                Process 2 finished:               *
    *  generator started at 2021-11-01T09:50:56.163481 *
    * generator finished at 2021-11-01T09:52:57.732244 *
    *                  total: 2m 1.57s                 *
    *          of which 2m 1.45s with binary_c         *
    *                  Ran 248 systems                 *
    *       with a total probability of 0.016309       *
    *         This thread had 0 failing systems        *
    *       with a total failed probability of 0       *
    *   Skipped a total of 0 zero-probability systems  *
    *                                                  *
    ****************************************************

    process 2 queue put output_dict

    [2021-11-01 09:52:57,736 DEBUG    Process-4] --- Process-2 is finished.

    process 2 return
    **************************************************************
    *    Population-2b66f805db424c48a1d29c45092b6e3c finished!   *
    *             The total probability is 0.0645906.            *
    * It took a total of 2m 1.85s to run 1000 systems on 4 cores *
    *                   = 8m 7.40s of CPU time.                  *
    *                Maximum memory use 941.109 MB               *
    **************************************************************

    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': '2b66f805db424c48a1d29c45092b6e3c', 'evolution_type': 'grid', 'failed_count': 0, 'failed_prob': 0, 'failed_systems_error_codes': [], 'errors_exceeded': False, 'errors_found': False, 'total_probability': 0.06459059967730083, 'total_count': 1000, 'start_timestamp': 1635760256.1239555, 'end_timestamp': 1635760377.9739752, 'total_mass_run': 4680.235689312423, 'total_probability_weighted_mass_run': 0.22611318083528548, '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.03502960360000004, 'unmerged': 0.019715467199999996}

    [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?