# Example use case: Hertzsprung-Russell diagrams

In this notebook we compute Hertzsprung-Russell diagrams (HRDs) of single and binary stars.

    [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_HRD")

## Setting up the Population object

First we set up a new population object. Our stars evolve to \(13.7\mathrm{Gyr}\), the age of the Universe, and we assume the metallicity \(Z=0.02\). These are rough approximations: a real population was born some finite time ago, so cannot possibly evolve to \(13.7\mathrm{Gyr}\), and stars are not really born with a metallicity of \(0.02\). These approximations only affect very low mass stars, so we assume all our stars have mass \(M\geq 1 \mathrm{M}_\odot\), and metallicity does not change evolution too much except in massive stars through the dependence of their winds on metallicity, so we limit our study to \(M\leq 10 \mathrm{M}_\odot\).

    [2]:

    # Create population object
    population = Population()

    # Setting values can be done via .set(<parameter_name>=<value>)
    # Values that are known to be binary_c_parameters are loaded into bse_options.
    # Those that are present in the default grid_options are set in grid_options
    # All other values that you set are put in a custom_options dict
    population.set(
        # binary_c physics 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
        tmp_dir=TMP_DIR,
        verbosity=1
    )

## Stellar Grid

We now construct a grid of stars, varying the mass from \(1\) to \(10\mathrm{M}_\odot\) in nine steps (so the masses are integers).

    [3]:

    import binarycpython.utils.distribution_functions
    # Set resolution and mass range that we simulate
    resolution = {"M_1": 10}
    massrange = (1, 11)

    population.add_grid_variable(
        name="M_1",
        longname="Primary mass", # == single-star mass
        valuerange=massrange,
        samplerfunc="const(1,10,9)", # space by unit masses
        probdist="1", # dprob/dm1 : we don't care, so just set it to 1
        dphasevol="dM_1",
        parameter_name="M_1",
        condition="",  # Impose a condition on this grid variable. Mostly for a check for yourself
        gridtype="edge"
    )

    Added grid variable: {
        "name": "M_1",
        "parameter_name": "M_1",
        "longname": "Primary mass",
        "valuerange": [
            1,
            11
        ],
        "samplerfunc": "const(1,10,9)",
        "precode": null,
        "postcode": null,
        "probdist": "1",
        "dphasevol": "dM_1",
        "condition": "",
        "gridtype": "edge",
        "branchpoint": 0,
        "branchcode": null,
        "topcode": null,
        "bottomcode": null,
        "grid_variable_number": 0
    }

## Setting logging and handling the output

We now construct the HRD output.

We choose stars prior to and including the thermally-pulsing asymptotic giant branch (TPAGB) phase that have \(>0.1\mathrm{M}_\odot\) of material in their outer hydrogen envelope (remember the core of an evolved star is made of helium or carbon/oxygen/neon). This prevents us showing the post-AGB phase which is a bit messy and we avoid the white-dwarf cooling track.

    [4]:

    custom_logging_statement = """
    Foreach_star(star)
    {
        if(star->stellar_type <= TPAGB &&
           star->mass - Outermost_core_mass(star) > 0.1)
        {
             double logTeff = log10(Teff_from_star_struct(star));
             double logL = log10(star->luminosity);
             double loggravity = log10(TINY+GRAVITATIONAL_CONSTANT*M_SUN*star->mass/Pow2(star->radius*R_SUN));
             Printf("HRD%d %30.12e %g %g %g %g\\n",
                    star->starnum, // 0
                    stardata->model.time, // 1
                    stardata->common.zero_age.mass[0], // 2 : note this is the primary mass
                    logTeff, // 3
                    logL, // 4
                    loggravity // 5
                    );

        }
    }
    """

    population.set(
        C_logging_code=custom_logging_statement
    )

    adding: C_logging_code=
    Foreach_star(star)
    {
        if(star->stellar_type <= TPAGB &&
           star->mass - Outermost_core_mass(star) > 0.1)
        {
             double logTeff = log10(Teff_from_star_struct(star));
             double logL = log10(star->luminosity);
             double loggravity = log10(TINY+GRAVITATIONAL_CONSTANT*M_SUN*star->mass/Pow2(star->radius*R_SUN));
             Printf("HRD%d %30.12e %g %g %g %g\n",
                    star->starnum, // 0
                    stardata->model.time, // 1
                    stardata->common.zero_age.mass[0], // 2 : note this is the primary mass
                    logTeff, // 3
                    logL, // 4
                    loggravity // 5
                    );

        }
    }
     to grid_options

The parse function must now catch lines that start with “HRDn”, where n is 0 (primary star) or 1 (secondary star, which doesn’t exist in single-star systems), and process the associated data.

    [5]:

    from binarycpython.utils.functions import datalinedict
    import re

    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", "zams_mass", "logTeff", "logL", "logg"]

        # Loop over the output.
        for line in output.splitlines():
            match = re.search('HRD(\d)',line)
            if match:
                nstar = match.group(1)

                # obtain the line of data in dictionary form
                linedata = datalinedict(line,parameters)

                # first time setup of the list of tuples
                if(len(self.grid_results['HRD'][nstar][linedata['zams_mass']])==0):
                    self.grid_results['HRD'][nstar][linedata['zams_mass']] = []

                # make the HRD be a list of tuples
                self.grid_results['HRD'][nstar][linedata['zams_mass']].append((linedata['logTeff'],
                                                                               linedata['logL']))

        # 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 0x14bf2250fdc0> to grid_options

## Evolving the grid

Now that we configured all the main parts of the population object, we can actually run the population! Doing this is straightforward: `population.evolve()`

This will start up the processing of all the systems. We can control how many cores are used by settings `num_cores`. By setting the `verbosity` of the population object to a higher value we can get a lot of verbose information about the run, but for now we will set it to 0.

There are many grid_options that can lead to different behaviour of the evolution of the grid. Please do have a look at the [grid options docs] (https://ri0005.pages.surrey.ac.uk/binary_c-python/grid_options_descriptions.html) for more details.

    [6]:

    # set number of threads
    population.set(
        # verbose output at level 1 is sufficient
        verbosity=1,
        # set number of threads (i.e. number of CPU cores we use)
        num_cores=1,
        )

    # Evolve the population - this is the slow, number-crunching step
    analytics = population.evolve()

    # Show the results (debugging)
    print (population.grid_results['HRD']['0'].keys())

    adding: verbosity=1 to grid_options
    adding: num_cores=1 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_HRD/binary_c_grid_8a7d52edc59346ce991063d9e05e9ba7.py [dry_run = True]
    Symlinked grid code to /tmp/binary_c_python-izzard/notebooks/notebook_HRD/binary_c_grid-latest0
    Loading grid code function from /tmp/binary_c_python-izzard/notebooks/notebook_HRD/binary_c_grid_8a7d52edc59346ce991063d9e05e9ba7.py
    Grid code loaded
    Dry run of the grid
    Grid has handled 9 stars with a total probability of 9
    **************************************
    *  Total starcount for this run is 9 *
    *       Total probability is 9       *
    **************************************

    [2021-11-01 09:49:52,836 DEBUG    MainProcess] --- setting up the system_queue_filler now

    Generating grid code
    Generating grid code

    [2021-11-01 09:49:52,838 DEBUG    Process-2] --- Setting up processor: process-0

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

    [2021-11-01 09:49:52,846 DEBUG    MainProcess] --- Signaling stop to processes
    [2021-11-01 09:49:53,800 DEBUG    Process-2] --- Process-0 is finishing.

    process 0 free memory and return
    ****************************************************
    *                Process 0 finished:               *
    *  generator started at 2021-11-01T09:49:52.838411 *
    * generator finished at 2021-11-01T09:49:53.802706 *
    *                   total: 0.96s                   *
    *           of which 0.94s with binary_c           *
    *                   Ran 9 systems                  *
    *           with a total probability of 9          *
    *         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:49:53,809 DEBUG    Process-2] --- Process-0 is finished.

    process 0 return
    **********************************************************
    *  Population-8a7d52edc59346ce991063d9e05e9ba7 finished! *
    *               The total probability is 9.              *
    *  It took a total of 1.21s to run 9 systems on 1 cores  *
    *                   = 1.21s of CPU time.                 *
    *              Maximum memory use 235.098 MB             *
    **********************************************************

    There were no errors found in this run.
    dict_keys([1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0])

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.

    [7]:

    print(analytics)

    {'population_name': '8a7d52edc59346ce991063d9e05e9ba7', 'evolution_type': 'grid', 'failed_count': 0, 'failed_prob': 0, 'failed_systems_error_codes': [], 'errors_exceeded': False, 'errors_found': False, 'total_probability': 9.0, 'total_count': 9, 'start_timestamp': 1635760192.8124797, 'end_timestamp': 1635760194.0176575, 'total_mass_run': 45.0, 'total_probability_weighted_mass_run': 45.0, 'zero_prob_stars_skipped': 0}

    [8]:

    # make an HRD using Seaborn and Pandas
    import seaborn as sns
    import pandas as pd
    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})

    hrd = population.grid_results['HRD']
    pd.set_option("display.max_rows", None, "display.max_columns", None)

    for nstar in sorted(hrd):
        print("star ",nstar)
        for zams_mass in sorted(hrd[nstar]):
            print("zams mass ",zams_mass)

            # get track data (list of tuples)
            track = hrd[nstar][zams_mass]

            # convert to Pandas dataframe
            data = pd.DataFrame(data=track,
                                columns = ['logTeff','logL'])

            # make seaborn plot
            p = sns.lineplot(data=data,
                             sort=False,
                             x='logTeff',
                             y='logL',
                             estimator=None)

            # set mass label at the zero-age main sequence (ZAMS) which is the first data point
            p.text(track[0][0],track[0][1],str(zams_mass))

    p.invert_xaxis()
    p.set_xlabel("$\log_{10} (T_\mathrm{eff} / \mathrm{K})$")
    p.set_ylabel("$\log_{10} (L/$L$_{☉})$")

    star  0
    zams mass  1.0
    zams mass  2.0
    zams mass  3.0
    zams mass  4.0
    zams mass  5.0
    zams mass  6.0
    zams mass  7.0
    zams mass  8.0
    zams mass  9.0

    [8]:

    Text(0, 0.5, '$\\log_{10} (L/$L$_{☉})$')

We now have an HRD. It took longer to make the plot than to run the stars with binary_c!

## Binary stars

Now we put a secondary star of mass \(0.5\mathrm{M}_\odot\) at a distance of \(10\mathrm{R}_\odot\) to see how this changes things. Then we rerun the population. At such short separations, we expect mass transfer to begin on or shortly after the main sequence.

    [9]:

    population.set(
        M_2 = 0.5, # Msun
        separation = 10, # Rsun
        multiplicity = 2, # binaries
    )
    population.clean()
    analytics = population.evolve()

    adding: M_2=0.5 to BSE_options
    adding: separation=10 to BSE_options
    adding: multiplicity=2 to BSE_options
    Custom logging library already loaded. Not setting them again.
    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_HRD/binary_c_grid_f4e5925effeb42e69b5be956f400ae41.py [dry_run = True]
    Symlinked grid code to /tmp/binary_c_python-izzard/notebooks/notebook_HRD/binary_c_grid-latest2
    Loading grid code function from /tmp/binary_c_python-izzard/notebooks/notebook_HRD/binary_c_grid_f4e5925effeb42e69b5be956f400ae41.py
    Grid code loaded
    Dry run of the grid
    Grid has handled 9 stars with a total probability of 9
    **************************************
    *  Total starcount for this run is 9 *
    *       Total probability is 9       *
    **************************************

    [2021-11-01 09:49:55,079 DEBUG    MainProcess] --- setting up the system_queue_filler now

    Generating grid code
    Generating grid code

    [2021-11-01 09:49:55,081 DEBUG    Process-4] --- Setting up processor: process-0

    Saving grid code to grid_options
    Writing grid code to /tmp/binary_c_python-izzard/notebooks/notebook_HRD/binary_c_grid_f4e5925effeb42e69b5be956f400ae41.py [dry_run = False]
    Symlinked grid code to /tmp/binary_c_python-izzard/notebooks/notebook_HRD/binary_c_grid-latest3
    Loading grid code function from /tmp/binary_c_python-izzard/notebooks/notebook_HRD/binary_c_grid_f4e5925effeb42e69b5be956f400ae41.py
    Grid code loaded

    [2021-11-01 09:49:55,093 DEBUG    MainProcess] --- Signaling stop to processes
    [2021-11-01 09:49:57,243 DEBUG    Process-4] --- Process-0 is finishing.

    process 0 free memory and return
    ****************************************************
    *                Process 0 finished:               *
    *  generator started at 2021-11-01T09:49:55.081024 *
    * generator finished at 2021-11-01T09:49:57.246671 *
    *                   total: 2.17s                   *
    *           of which 2.14s with binary_c           *
    *                   Ran 9 systems                  *
    *           with a total probability of 9          *
    *         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:49:57,256 DEBUG    Process-4] --- Process-0 is finished.

    process 0 return
    **********************************************************
    *  Population-f4e5925effeb42e69b5be956f400ae41 finished! *
    *               The total probability is 9.              *
    *  It took a total of 2.49s to run 9 systems on 1 cores  *
    *                   = 2.49s of CPU time.                 *
    *              Maximum memory use 294.551 MB             *
    **********************************************************

    There were no errors found in this run.

    [10]:

    hrd = population.grid_results['HRD']

    for nstar in sorted(hrd):
        print("star ",nstar)

        if nstar == '0': # choose only primaries

            for zams_mass in sorted(hrd[nstar]):
                print("zams mass ",zams_mass)

                # get track data (list of tuples)
                track = hrd[nstar][zams_mass]

                # convert to Pandas dataframe
                data = pd.DataFrame(data=track,
                                    columns = ['logTeff','logL'])

                # make seaborn plot
                p = sns.lineplot(data=data,
                                 sort=False,
                                 x='logTeff',
                                 y='logL',
                                 estimator=None)

                # set mass label at the zero-age main sequence (ZAMS) which is the first data point
                p.text(track[0][0],track[0][1],str(zams_mass))

    p.invert_xaxis()
    p.set_xlabel("$\log_{10} (T_\mathrm{eff} / \mathrm{K})$")
    p.set_ylabel("$\log_{10} (L/$L$_{☉})$")

    star  0
    zams mass  1.0
    zams mass  2.0
    zams mass  3.0
    zams mass  4.0
    zams mass  5.0
    zams mass  6.0
    zams mass  7.0
    zams mass  8.0
    zams mass  9.0
    star  1

    [10]:

    Text(0, 0.5, '$\\log_{10} (L/$L$_{☉})$')

We plot here the track for the primary star only. You can see immediately where stars merge on the main sequence: the tracks move very suddenly where usually evolution on the main sequence is smooth.

If we now set the separation to be longer, say \(100\mathrm{R}_\odot\), mass transfer should happen on the giant branch. We also set the secondary mass to be larger, \(1\mathrm{M}_\odot\), so that the interaction is stronger.

    [11]:

    population.set(
        M_2 = 1, # Msun
        separation = 100, # Rsun
        multiplicity = 2, # binaries
        alpha_ce = 1.0, # make common-envelope evolution quite efficient
    )
    population.clean()
    analytics = population.evolve()

    adding: M_2=1 to BSE_options
    adding: separation=100 to BSE_options
    adding: multiplicity=2 to BSE_options
    adding: alpha_ce=1.0 to BSE_options
    Custom logging library already loaded. Not setting them again.
    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_HRD/binary_c_grid_26171201413c43f5b91552a374879940.py [dry_run = True]
    Symlinked grid code to /tmp/binary_c_python-izzard/notebooks/notebook_HRD/binary_c_grid-latest4
    Loading grid code function from /tmp/binary_c_python-izzard/notebooks/notebook_HRD/binary_c_grid_26171201413c43f5b91552a374879940.py
    Grid code loaded
    Dry run of the grid
    Grid has handled 9 stars with a total probability of 9
    **************************************
    *  Total starcount for this run is 9 *
    *       Total probability is 9       *
    **************************************

    [2021-11-01 09:49:57,844 DEBUG    MainProcess] --- setting up the system_queue_filler now

    Generating grid code
    Generating grid code

    [2021-11-01 09:49:57,846 DEBUG    Process-6] --- Setting up processor: process-0

    Saving grid code to grid_options
    Writing grid code to /tmp/binary_c_python-izzard/notebooks/notebook_HRD/binary_c_grid_26171201413c43f5b91552a374879940.py [dry_run = False]
    Symlinked grid code to /tmp/binary_c_python-izzard/notebooks/notebook_HRD/binary_c_grid-latest5
    Loading grid code function from /tmp/binary_c_python-izzard/notebooks/notebook_HRD/binary_c_grid_26171201413c43f5b91552a374879940.py
    Grid code loaded

    [2021-11-01 09:49:57,858 DEBUG    MainProcess] --- Signaling stop to processes
    [2021-11-01 09:49:59,413 DEBUG    Process-6] --- Process-0 is finishing.

    process 0 free memory and return
    ****************************************************
    *                Process 0 finished:               *
    *  generator started at 2021-11-01T09:49:57.846045 *
    * generator finished at 2021-11-01T09:49:59.416141 *
    *                   total: 1.57s                   *
    *           of which 1.55s with binary_c           *
    *                   Ran 9 systems                  *
    *           with a total probability of 9          *
    *         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:49:59,424 DEBUG    Process-6] --- Process-0 is finished.

    process 0 return
    **********************************************************
    *  Population-26171201413c43f5b91552a374879940 finished! *
    *               The total probability is 9.              *
    *  It took a total of 1.89s to run 9 systems on 1 cores  *
    *                   = 1.89s of CPU time.                 *
    *              Maximum memory use 312.051 MB             *
    **********************************************************

    There were no errors found in this run.

    [12]:

    hrd = population.grid_results['HRD']

    for nstar in sorted(hrd):
        print("star ",nstar)

        if nstar == '0': # choose only primaries

            for zams_mass in sorted(hrd[nstar]):
                print("primary zams mass ",zams_mass)

                # get track data (list of tuples)
                track = hrd[nstar][zams_mass]

                # convert to Pandas dataframe
                data = pd.DataFrame(data=track,
                                    columns = ['logTeff','logL'])

                # make seaborn plot
                p = sns.lineplot(data=data,
                                 sort=False,
                                 x='logTeff',
                                 y='logL',
                                 estimator=None)

                # set mass label at the zero-age main sequence (ZAMS) which is the first data point
                p.text(track[0][0],track[0][1],str(zams_mass))

    p.invert_xaxis()
    p.set_xlabel("$\log_{10} (T_\mathrm{eff} / \mathrm{K})$")
    p.set_ylabel("$\log_{10} (L/$L$_{☉})$")

    star  0
    primary zams mass  1.0
    primary zams mass  2.0
    primary zams mass  3.0
    primary zams mass  4.0
    primary zams mass  5.0
    primary zams mass  6.0
    primary zams mass  7.0
    primary zams mass  8.0
    primary zams mass  9.0
    star  1

    [12]:

    Text(0, 0.5, '$\\log_{10} (L/$L$_{☉})$')

You now see the interaction in the jerky red-giant tracks where the stars interact. These probably, depending on the mass ratio at the moment of interaction, go through a common-envelope phase. The system can merge (most of the above do) but not all. The interaction is so strong on the RGB of the \(1\mathrm{M}_\odot\) star that the stellar evolution is terminated before it reaches the RGB tip, so it never ignites helium. This is how helium white dwarfs are probably made.

We can also plot the secondary stars’ HRD. Remember, the primary is star 0 in binary_c, while the secondary is star 1. That’s because all proper programming languages start counting at 0. We change the parsing function a little so we can separate the plots of the secondaries according to their primary mass.

    [13]:

    hrd = population.grid_results['HRD']

    for nstar in sorted(hrd):
        print("star ",nstar)

        if nstar == '1': # choose only secondaries

            for zams_mass in sorted(hrd[nstar]):
                print("primary zams mass ",zams_mass)

                # get track data (list of tuples)
                track = hrd[nstar][zams_mass]

                # convert to Pandas dataframe
                data = pd.DataFrame(data=track,
                                    columns = ['logTeff','logL'])

                # make seaborn plot
                p = sns.lineplot(data=data,
                                 sort=False,
                                 x='logTeff',
                                 y='logL',
                                 estimator=None)

    p.invert_xaxis()
    p.set_xlabel("$\log_{10} (T_\mathrm{eff} / \mathrm{K})$")
    p.set_ylabel("$\log_{10} (L/$L$_{☉})$")

    star  0
    star  1
    primary zams mass  1.0
    primary zams mass  2.0
    primary zams mass  3.0
    primary zams mass  4.0
    primary zams mass  5.0
    primary zams mass  6.0
    primary zams mass  7.0
    primary zams mass  8.0
    primary zams mass  9.0

    [13]:

    Text(0, 0.5, '$\\log_{10} (L/$L$_{☉})$')

Remember, all these stars start with a \(1\mathrm{M}_\odot\) binary, which begins at \(\log_{10}(T_\mathrm{eff}/\mathrm{K})\sim 3.750\), \(\log_{10}L/\mathrm{L}_\odot \sim 0\). The \(1\mathrm{M}_\odot\)-\(1\mathrm{M}_\odot\) binary evolves like two single stars until they interact up the giant branch at about \(\log_{10} (L/\mathrm{L}_\odot) \sim 2.5\), the others interact long before they evolve very far on the main sequence: you can just about see their tracks at the very start.

This is, of course, a very simple introduction to what happens in binaries. We haven’t talked about the remnants that are produced by interactions. When the stars do evolve on the giant branch, white dwarfs are made which can go on to suffer novae and (perhaps) thermonuclear explosions. The merging process itself leads to luminosus red novae and, in the case of neutron stars and black holes, kilonovae and gravitational wave events.

    [ ]: