Evolve individual binaries 🐞

It can be extremely useful to evolve a single binary from specific initial conditions or to start at a specific evolutionary state. This is useful for debugging binaries and checking if your custom steps or flow are correctly working.

This tutorial will cover:

• How to initialise a zero-age main sequence binary (ZAMS).

• How to re-evolve a binary in an existing population.

• How to evolve a binary from specific initial states.

Evolve a ZAMS binary

To evolve a binary from ZAMS, we will:

1. Load the standard simulation properties.

2. Load the steps.

3. Initialise the binary.

4. Evolve it.

import os
import shutil
from posydon.config import PATH_TO_POSYDON

path_to_params = os.path.join(PATH_TO_POSYDON, "posydon/popsyn/population_params_default.ini")
shutil.copyfile(path_to_params, './population_params.ini')

# load the function to load the simulation properties from the ini file
from posydon.popsyn.io import simprop_kwargs_from_ini
from posydon.binary_evol.simulationproperties import SimulationProperties

# Load the simulation properties from the default ini file.
sim_kwargs = simprop_kwargs_from_ini('population_params.ini')
# manually add the metallicity to each step that requires it
metallicity = {'metallicity':1}

sim_kwargs['step_HMS_HMS'][1].update(metallicity)
sim_kwargs['step_CO_HeMS'][1].update(metallicity)
sim_kwargs['step_CO_HMS_RLO'][1].update(metallicity)
sim_kwargs['step_CO_HeMS_RLO'][1].update(metallicity)
sim_kwargs['step_detached'][1].update(metallicity)
sim_kwargs['step_disrupted'][1].update(metallicity)
sim_kwargs['step_merged'][1].update(metallicity)
sim_kwargs['step_initially_single'][1].update(metallicity)

sim_prop = SimulationProperties(**sim_kwargs)
# Load the steps and required data
sim_prop.load_steps(verbose=True)

# load the binary and single star classes
from posydon.binary_evol.singlestar import SingleStar
from posydon.binary_evol.binarystar import BinaryStar

STAR1 = SingleStar(**{'mass': 30.782576,
                      'state': 'H-rich_Core_H_burning'})
STAR2 = SingleStar(**{'mass':20.273864,
                      'state': 'H-rich_Core_H_burning'})

BINARY = BinaryStar(STAR1, STAR2,
                    **{'time': 0.0, 'state': 'detached', 'event': 'ZAMS', 'orbital_period':3513.150157, 'eccentricity': 0.0},
                    properties = sim_prop)

# Note: depending on the evolution of the binary, you might get some warnings about the Roche lobe calculation.
BINARY.evolve()

You’re now able to access the binary and evolutionary information. Using the function, BinaryStar.to_df() we can inspect the evolution of the systm.

# we will add the step names to the dataframe
col = ['time', 'step_names', 'state', 'event', 'orbital_period', 'eccentricity', 'S1_state', 'S2_state', 'S1_mass', 'S2_mass']

BINARY.to_df(extra_columns={'step_names':'string'})[col]

Re-evolving the Binary

The above binary might have been disrupted in the first or second SN due to the strength of the supernova kick. You can restore the binary completely or to a specific state and re-evolve it.

# restore to the original state
BINARY.restore()
BINARY.to_df()

BINARY.evolve()

Even thought we reset the binary, the natal kick properties stored in the natal kick array are not reset. This is because we want to be able to re-evolve the binary to the same final state. In case you want to reset the natal kick array values, then set the list to None values.

# kick magnitude km/s, azimuthal angle rad, polar angle rad and mean anomaly rad (TODO check documentation)
print(BINARY.star_1.natal_kick_array)
print(BINARY.star_2.natal_kick_array)

# This will be exactly the same evolution as the previous one
BINARY.to_df(extra_columns={'step_names':'string'})[col]

Changing the kick

We will restore the BINARY to its initial state and increase the first kick velocity to disrupt the binary after the first supernova.

BINARY.restore()
BINARY.to_df()

# [velocity, azimuthal angle, polar angle, phase]
BINARY.star_1.natal_kick_array = [1000., 0.858129274334538, 1.9157148786534735, 1.8675467897282945]
BINARY.star_2.natal_kick_array = [None, None, None, None] # default
print(BINARY.star_1.natal_kick_array)
print(BINARY.star_2.natal_kick_array)

BINARY.evolve()
BINARY.to_df(extra_columns={'step_names':'string'})[col]

Evolve from specific step

You can provide the row number (starting at 0), to reset the binary to that specific evolutionary phase. This can be useful for debugging a specific evolutionary state.

BINARY.restore(3)
BINARY.to_df()

Loading a binary from an existing population

Assume we run the population somewhere else and you have the h5 file associated with it, and you want to load the population in memory and re-evolve a binary. We will have to load the population_params.ini file associated with the Population file to re-evolve the binary in the same way.

We will use the BBH_contact.h5 file in the example dataset and use the standard population_params.ini SimulationProperties from earlier in this notebook.

from posydon.popsyn.synthetic_population import Population


data_path = f'PATH/TO/TUTORIAL/POPULATIONS/'

pop = Population(f'{data_path}/BBH_contact.h5')

# Because the original population contains two extre columns. These have to be defined manually here.
# Otherwise, the reseting of the binaries fails.
BINARY = BinaryStar.from_df(pop.history[0],
                            extra_columns={'step_names':'string','step_times':'float'})

# you have to set the simulation properties for the binary
BINARY.properties = sim_prop

# Set the natal kick arrays from the oneline to get the exact same evolution as the original binary
BINARY.star_1.natal_kick_array = pop.oneline[0][['S1_natal_kick_array_0', 'S1_natal_kick_array_1', 'S1_natal_kick_array_2', 'S1_natal_kick_array_3']].values.flatten()
BINARY.star_2.natal_kick_array = pop.oneline[0][['S2_natal_kick_array_0', 'S2_natal_kick_array_1', 'S2_natal_kick_array_2', 'S2_natal_kick_array_3']].values.flatten()

BINARY.to_df(extra_columns={'step_names':'string'})[col]

BINARY.restore()

BINARY.evolve()
BINARY.to_df(extra_columns={'step_names':'string'})[col]

Evolving a Binary Star starting from an arbitrary state

If you want to evolve a custom binary star, you can do so by crafting yourself the BinaryStar object and providing the SimulationProperties object. For example, let’s evolve a neutron star with a low mass helium star in Roche lobe overflow. And start a binary in the detached step, which requires a bit of additional input data.

We will use the standard SimulationProperties, loaded earlier in this notebook.

Caution!

While you can evolve a binary from an arbitrary state, you will need to provide data on the internal structure of the star, if you’re starting in the detached step. Otherwise, the matching to a single star model will not work!

# we initialise the binary to be tight and Roche lobe overflowing in a circular orbit.
binary = BinaryStar(star_1 = SingleStar(**{'state' : 'NS',
                                           'mass' : 1.1,
                                           'spin' : 0.,}),
                    star_2 = SingleStar(**{'state' : 'H-rich_Core_H_burning',
                                           'mass' : 2.5,
                                           'natal_kick_array' : [10., 0., 0., 0.]}),
                    **{'time' : 0.,
                       'state' : 'RLO2',
                       'event' : 'oRLO2',
                       'orbital_period' : 1.,
                       'eccentricity' : 0.},
                    properties = sim_prop,
                    )
binary.evolve()
binary.to_df(extra_columns={'step_names':'string'})[col]

# we perform the same evolution, but with a stripped helium star to go to the CO_HeMS_RLO grid
binary = BinaryStar(star_1 = SingleStar(**{'state' : 'NS',
                                           'mass' : 1.1,
                                           'spin' : 0.,}),
                    star_2 = SingleStar(**{'state' : 'stripped_He_Core_He_burning',
                                           'mass' : 2.5,
                                           'natal_kick_array' : [10., 0., 0., 0.]}),
                    **{'time' : 0.,
                       'state' : 'RLO2',
                       'event' : 'oRLO2',
                       'orbital_period' : 1.,
                       'eccentricity' : 0.},
                    properties = sim_prop,
                    )
binary.evolve()
binary.to_df(extra_columns={'step_names':'string'})[col]

We will change the binary to a much wider system and have it start in the detached state. This requires additional starting information to allow the star to be matched to the single star grids.

from posydon.utils.constants import Zsun
from posydon.utils.common_functions import orbital_separation_from_period

import numpy as np
Z = 1.0

# Setup the central abundances
zams_table = {2.: 2.915e-01,
                      1.: 2.703e-01,
                      0.45: 2.586e-01,
                      0.2: 2.533e-01,
                      0.1: 2.511e-01,
                      0.01: 2.492e-01,
                      0.001: 2.49e-01,
                      0.0001: 2.49e-01}

Y = zams_table[Z]
Z = Z*Zsun
X = 1 - Y - Z

STAR1 = SingleStar(**{'mass':1.2,
                      'state': 'NS'})

STAR2 = SingleStar(**{'mass': 17.782576,
                      'state': 'H-rich_Core_H_burning',
                      # add the metallicity and central abundances
                      'metallicity':Z,
                      'center_h1':X,
                      'center_he4':Y,
                      # add a numerical value for the radius
                      'log_R': np.nan,
                      # add the helium core mass
                      'he_core_mass': 0.0,
                      })

binary = BinaryStar(STAR1, STAR2,
                    **{'time' : 0.,
                       'state' : 'detached',
                       'event' : None,
                       'orbital_period' : 5000.,
                        # calculate the separation; current bug that the separation is not automatically calculated if orbital period is given
                       'separation': orbital_separation_from_period(5000., 17.782576, 1.2),
                       'eccentricity' : 0.},
                    properties = sim_prop,
                    )
binary.evolve()
binary.to_df(extra_columns={'step_names':'string'})[col]