Example use case: Binary black hole systems

In this notebook we will set up binary systems to find binary black hole systems. We will make use of the log_files to find out whether a system has become a BBH and to store the summary of their evolution.

We use the run_system function to individual binary systems so we can play around with the input settings and see how that affects the formation of the BHBH system. There are other, more appropriate, ways to find populations of binary black hole systems. A good start for that is setting up a custom logging code and a parse function, and then run a grid of systems.

[1]:
import os
from binarycpython.utils.functions import temp_dir
from binarycpython.utils.run_system_wrapper import run_system

TMP_DIR = temp_dir("notebooks", "notebook_BHBH")
VERBOSITY = 0

To run a system we can use the run_system function, which can parse function arguments for binary_c and will run a system.

[2]:
# set filename
log_filename = "log_file.txt"

# Run via run_system
output = run_system(M_1=60,
                    BH_prescription='BH_BELCZYNSKI',
                    log_filename=log_filename,
                    api_log_filename_prefix=TMP_DIR)

# Readout the contents
with open(log_filename, 'r') as f:
    print(f.read())
      TIME      M1       M2   K1  K2           SEP        PER   ECC  R1/ROL1 R2/ROL2  TYPE RANDOM_SEED=13868 RANDOM_COUNT=0
     0.0000   60.000    0.000  1  15            -1       -1   -1.00   0.000   0.000  "INITIAL "
     6.6492   26.756    0.000  2  15            -1       -1   -1.00   0.000   0.000  "OFF_MS"
     6.6492   26.756    0.000  2  15            -1       -1   -1.00   0.000   0.000  "TYPE_CHNGE"
     6.6588   26.687    0.000  4  15            -1       -1   -1.00   0.000   0.000  "TYPE_CHNGE"
     7.2135    9.972    0.000  7  15            -1       -1   -1.00   0.000   0.000  "TYPE_CHNGE"
     7.5298    7.325    0.000  8  15            -1       -1   -1.00   0.000   0.000  "TYPE_CHNGE"
     7.5700    2.903    0.000 14  15            -1       -1   -1.00   0.000   0.000  Randbuf=14456 - Mers(0)=0.415747 - Mers(1)=0.403489 - Mers(2)=0.273236 - Mers(3)=0.627902
     7.5700    2.903    0.000 14  15            -1       -1   -1.00   0.000   0.000  SN kick Ib/c (SN type 11 11, pre-explosion M=6.92015 Mc"CO"=5.31213 type=8) -> kick 1(190) vk=223.39 vr=0 omega=3.94523 phi=-0.470721 -> vn=223.39 ; final sep 0 ecc -1 (random count 0) - Runaway v=(-155.054,-143.325,72.9341) |v|=223.39 : companion v=(0,0,0), |v|=0 ;
     7.5700    2.903    0.000 14  15            -1       -1   -1.00   0.000   0.000  "TYPE_CHNGE"
     7.5700    2.903    0.000 14  15            -1       -1   -1.00   0.000   0.000  "SN"
 15000.0000    2.903    0.000 14  15            -1       -1   -1.00   0.000   0.000  "MAX_TIME"

From the output of this file we can see that the stellar types (K1, K2) are not entirely what we are looking for. One of them is 15, which means the system must have merged. To prevent this, we can star the system with a wider initial orbit, so that the stars do not interact (except by a little wind accretion).

[3]:
# run via run_system
output = run_system(M_1=60,
                    M_2=40,
                    orbital_period=2000, # days
                    BH_prescription='BH_BELCZYNSKI',
                    log_filename=log_filename,
                    wind_mass_loss='WIND_ALGORITHM_BINARY_C_2020',
                    api_log_filename_prefix=TMP_DIR)

# Readout contents
with open(log_filename, 'r') as f:
    print(f.read())
      TIME      M1       M2   K1  K2           SEP        PER   ECC  R1/ROL1 R2/ROL2  TYPE RANDOM_SEED=71985 RANDOM_COUNT=0
     0.0000   60.000   40.000  1   1        3101.2     5.48 y  0.00   0.009   0.009  "INITIAL "
     0.0000   60.000   40.000  1   1        3101.2     5.48 y  0.00   0.009   0.009  "BEG_SYMB"
     4.3583   42.081   35.930  1   1          3971     8.98 y  0.00   0.023   0.014  "Start tidal lock 1"
     4.3928   41.632   35.894  1   1        3995.8      9.1 y  0.00   0.023   0.014  "End tidal lock 1"
     6.4219   27.615   27.603  1   1        5603.4     17.9 y  0.00   0.013   0.013  "Start tidal lock 2"
     6.5403   27.157   27.151  1   1        5696.9     18.5 y  0.00   0.012   0.012  "End tidal lock 2"
     6.6475   26.763   26.731  2   1        5791.5     19.1 y  0.00   0.010   0.011  "TYPE_CHNGE"
     6.6548   26.730   26.729  2   1        5794.3     19.1 y  0.00   0.203   0.010  "Start tidal lock 2"
     6.6550   26.728   26.729  2   1        5794.5     19.1 y  0.00   0.225   0.010  "q-inv"
     6.6552   26.725   26.729  2   2        5795.5     19.1 y  0.00   0.246   0.010  "OFF_MS"
     6.6552   26.725   26.729  2   2        5795.5     19.1 y  0.00   0.246   0.010  "TYPE_CHNGE"
     6.6552   26.725   26.729  2   2        5795.5     19.1 y  0.00   0.246   0.010  "End tidal lock 2"
     6.6556   26.721   26.728  2   2          5796     19.1 y  0.00   0.282   0.011  "Start tidal lock 2"
     6.6559   26.716   26.727  2   2        5796.6     19.1 y  0.00   0.321   0.013  "End tidal lock 2"
     6.6571   26.691   26.724  4   2        5798.4     19.2 y  0.00   0.546   0.021  "TYPE_CHNGE"
     6.6649   26.502   26.683  4   4          5812     19.3 y  0.00   0.548   0.543  "TYPE_CHNGE"
     7.2663   10.178   10.480  4   4        8828.2     57.9 y  0.00   0.000   0.460  "END_SYMB"
     7.2670   10.171   10.458  7   4        8837.8       58 y  0.00   0.000   0.455  "TYPE_CHNGE"
     7.2733   10.111   10.272  7   4        8921.3     59.2 y  0.00   0.000   0.265  "BEG_SYMB"
     7.2801   10.033   10.166  7   7          8987     60.1 y  0.00   0.000   0.000  "TYPE_CHNGE"
     7.4925    7.910    7.995  8   7         11434     97.2 y  0.00   0.000   0.000  "TYPE_CHNGE"
     7.4989    7.850    7.942  8   8         11516     98.6 y  0.00   0.000   0.000  "TYPE_CHNGE"
     7.5282    3.966    7.616 14   8       -16.539       -1   -1.00   0.000   0.000  Randbuf=54990 - Mers(0)=0.431329 - Mers(1)=0.844304 - Mers(2)=0.737837 - Mers(3)=0.957392 - Mers(4)=0.745668
     7.5282    3.966    7.616 14   8       -16.539       -1   -1.00   0.000   0.000  SN kick Ib/c (SN type 11 11, pre-explosion M=7.49451 Mc"CO"=5.76448 type=8) -> kick 1(190) vk=366.674 vr=15.5058 omega=4.68517 phi=1.15497 -> vn=366.077 ; final sep -16.5393 ecc -1 (random count 0) - Runaway v=(-2.16623,-147.729,-335.088) |v|=366.214 : companion v=(-7.68852,-0.172227,-0.112321), |v|=7.69127 ;
     7.5282    3.966    7.616 14   8       -16.539       -1   -1.00   0.000   0.000  "TYPE_CHNGE"
     7.5282    3.966    7.616 14   8       -16.539       -1   -1.00   0.000   0.000  "DISRUPT "
     7.5282    3.966    7.616 14   8       -16.539       -1   -1.00   0.000   0.000  "END_SYMB"
     7.5282    3.966    7.616 14   8       -16.539       -1   -1.00   0.000   0.000  "SN"
     7.5344    3.966    4.029 14  14       -16.539       -1   -1.00   0.000   0.000  Mers(5)=0.787956 - Mers(6)=0.497544 - Mers(7)=0.721903 - Mers(8)=0.10755
     7.5344    3.966    4.029 14  14       -16.539       -1   -1.00   0.000   0.000  SN kick Ib/c (SN type 11 11, pre-explosion M=7.526 Mc"CO"=5.7893 type=8) -> kick 1(190) vk=410.158 vr=0 omega=0.675756 phi=0.459842 -> vn=410.158 ; final sep -16.5393 ecc -1 (random count 5) - Runaway v=(312.331,229.727,113.746) |v|=410.158 : companion v=(-6.54865e-51,-8.54769e-104,5.39903e-104), |v|=6.54865e-51 ;
     7.5344    3.966    4.029 14  14       -16.539       -1   -1.00   0.000   0.000  "TYPE_CHNGE"
     7.5344    3.966    4.029 14  14       -16.539       -1   -1.00   0.000   0.000  "SN"
 15000.0000    3.966    4.029 14  14       -16.539       -1   -1.00   0.000   0.000  "MAX_TIME"

So let’s write a function to detect BHBH pairs.

[4]:
def detect_BHBH(filename):
    """
    Function to detect whether the system is a binary black hole system, based on the logfile
    """

    with open(filename, 'r') as f:
        log_file_lines = f.readlines()

    # Check if the content has more than 1 line
    if not len(log_file_lines) > 1:
        print("not long enough")
        return None

    # loop over lines
    for line in log_file_lines[1:]:
        data = line.split()

        # Capture stellar types
        stellar_type1 = int(data[3])
        stellar_type2 = int(data[4])

        # remember: stellar type 14 == BH
        if stellar_type1 == 14 and stellar_type2 == 14:

            # BHBH system
            names = ['time','M1','M2','K1','K2','separation']
            d = {}
            for name in names:
                d[name] = float(data.pop(0))

            return_dict = {
                'system_properties': d,
                'log_file_contents': ''.join(log_file_lines)
            }

            return return_dict

data = detect_BHBH(log_filename)

if data is None or len(data) == 0:
    print("Oops: no data found")
else:
    if data['system_properties']['separation'] < 0.0:
        print("BHBH system is unbound")
    else:
        print("BHBH system is bound")
BHBH system is unbound

Now that we have a function to detect a BHBH system, we want to set up a search for these systems.

[5]:
def search_for_BHBH(maxcount, verbosity=0, **opts):
    """
    Function search for a BHBH system
    """

    found = False
    count = 0

    while found == False and count < maxcount:
        count = count + 1

        if verbosity: print("system {} / {}".format(count,maxcount))

        output = run_system(**opts,
                            log_filename=log_filename,
                            api_log_filename_prefix=TMP_DIR,
                            )
        data = detect_BHBH(log_filename)

        if data is None or len(data) == 0:
            # no BHBH found
            found = False

        else:
            if data['system_properties']['separation'] > 0.0:
                # found bound BHBH system
                print('Found bound BHBH system')
                found = True

                return data


    if found == False:
        print("No BHBH systems found")
        return {}
[6]:
args = dict(
            M_1=60,
            M_2=40,
            orbital_period=2000, # days
            BH_prescription='BH_FRYER12_RAPID',
            wind_mass_loss='WIND_ALGORITHM_BINARY_C_2020',
           )

search_result = search_for_BHBH(100, verbosity=VERBOSITY, **args)
No BHBH systems found

How can we help the system become a BHBH merger? We can try turning off SN kicks.

[7]:
args = dict(
            M_1=60,
            M_2=40,
            orbital_period=2000, # days
            BH_prescription='BH_BELCZYNSKI',
            wind_mass_loss='WIND_ALGORITHM_BINARY_C_2020',
            sn_kick_dispersion_II=0,
            sn_kick_dispersion_IBC=0,
            sn_kick_dispersion_GRB_COLLAPSAR=0,
           )

search_results = search_for_BHBH(100, verbosity=VERBOSITY, **args)
if search_results:
    print(search_results['log_file_contents'])
Found bound BHBH system
      TIME      M1       M2   K1  K2           SEP        PER   ECC  R1/ROL1 R2/ROL2  TYPE RANDOM_SEED=16379 RANDOM_COUNT=0
     0.0000   60.000   40.000  1   1        3101.2     5.48 y  0.00   0.009   0.009  "INITIAL "
     0.0000   60.000   40.000  1   1        3101.2     5.48 y  0.00   0.009   0.009  "BEG_SYMB"
     4.3583   42.081   35.930  1   1          3971     8.98 y  0.00   0.023   0.014  "Start tidal lock 1"
     4.3928   41.632   35.894  1   1        3995.8      9.1 y  0.00   0.023   0.014  "End tidal lock 1"
     6.4219   27.615   27.603  1   1        5603.4     17.9 y  0.00   0.013   0.013  "Start tidal lock 2"
     6.5403   27.157   27.151  1   1        5696.9     18.5 y  0.00   0.012   0.012  "End tidal lock 2"
     6.6475   26.763   26.731  2   1        5791.5     19.1 y  0.00   0.010   0.011  "TYPE_CHNGE"
     6.6548   26.730   26.729  2   1        5794.3     19.1 y  0.00   0.203   0.010  "Start tidal lock 2"
     6.6550   26.728   26.729  2   1        5794.5     19.1 y  0.00   0.225   0.010  "q-inv"
     6.6552   26.725   26.729  2   2        5795.5     19.1 y  0.00   0.246   0.010  "OFF_MS"
     6.6552   26.725   26.729  2   2        5795.5     19.1 y  0.00   0.246   0.010  "TYPE_CHNGE"
     6.6552   26.725   26.729  2   2        5795.5     19.1 y  0.00   0.246   0.010  "End tidal lock 2"
     6.6556   26.721   26.728  2   2          5796     19.1 y  0.00   0.282   0.011  "Start tidal lock 2"
     6.6559   26.716   26.727  2   2        5796.6     19.1 y  0.00   0.321   0.013  "End tidal lock 2"
     6.6571   26.691   26.724  4   2        5798.4     19.2 y  0.00   0.546   0.021  "TYPE_CHNGE"
     6.6649   26.502   26.683  4   4          5812     19.3 y  0.00   0.548   0.543  "TYPE_CHNGE"
     7.2663   10.178   10.480  4   4        8828.2     57.9 y  0.00   0.000   0.460  "END_SYMB"
     7.2670   10.171   10.458  7   4        8837.8       58 y  0.00   0.000   0.455  "TYPE_CHNGE"
     7.2733   10.111   10.272  7   4        8921.3     59.2 y  0.00   0.000   0.265  "BEG_SYMB"
     7.2801   10.033   10.166  7   7          8987     60.1 y  0.00   0.000   0.000  "TYPE_CHNGE"
     7.4925    7.910    7.995  8   7         11434     97.2 y  0.00   0.000   0.000  "TYPE_CHNGE"
     7.4989    7.850    7.942  8   8         11516     98.6 y  0.00   0.000   0.000  "TYPE_CHNGE"
     7.5282    3.966    7.616 14   8         17247      211 y  0.30   0.000   0.000  Randbuf=81114 - Mers(0)=0.159326 - Mers(1)=0.970141 - Mers(2)=0.340351 - Mers(3)=0.132638 - Mers(4)=0.502053
     7.5282    3.966    7.616 14   8         17247      211 y  0.30   0.000   0.000  SN kick Ib/c (SN type 11 11, pre-explosion M=7.49451 Mc"CO"=5.76448 type=8) -> kick 1(0) vk=0 vr=15.5058 omega=3.15449 phi=-0.825261 -> vn=15.5058 ; final sep 17247.1 ecc 0.304595 (random count 0) - Runaway v=(0,0,0) |v|=0 : companion v=(0,0,0), |v|=0 ;
     7.5282    3.966    7.616 14   8         17247      211 y  0.30   0.000   0.000  "TYPE_CHNGE"
     7.5282    3.966    7.616 14   8         17247      211 y  0.30   0.000   0.000  "SN"
     7.5300    3.966    7.592 14   8         17209      210 y  0.30   0.000   0.000  "END_SYMB"
     7.5344    3.966    4.029 14  14         22992      391 y  0.13   0.000   0.000  Mers(5)=0.584197 - Mers(6)=0.452365 - Mers(7)=0.0846357 - Mers(8)=0.848993 - Mers(9)=0.111944
     7.5344    3.966    4.029 14  14         22992      391 y  0.13   0.000   0.000  SN kick Ib/c (SN type 11 11, pre-explosion M=7.526 Mc"CO"=5.7893 type=8) -> kick 1(0) vk=0 vr=8.44077 omega=0.703362 phi=0.772581 -> vn=8.44077 ; final sep 22992.4 ecc 0.130746 (random count 5) - Runaway v=(0,0,0) |v|=0 : companion v=(0,0,0), |v|=0 ;
     7.5344    3.966    4.029 14  14         22992      391 y  0.13   0.000   0.000  "TYPE_CHNGE"
     7.5344    3.966    4.029 14  14         22992      391 y  0.13   0.000   0.000  "SN"
     7.5444    3.966    4.029 14  14         22992      391 y  0.13   0.000   0.000  "END_SYMB"
 15000.0000    3.966    4.029 14  14         22992      391 y  0.13   0.000   0.000  "MAX_TIME"

You should now have found a BHBH system but usually they have very wide orbits. This is caused by mass loss before, and during, the supernova.

We can reduce the former by setting the wind mass loss to zero. This can be done unphysically by setting wind_mass_loss=0, or you can reduce the massive-star winds by setting the metallicity to be small, e.g. \(10^{-3}\).

[8]:
args = dict(
            M_1=60,
            M_2=40,
            metallicity=0.001,
            orbital_period=2000, # days
            BH_prescription='BH_BELCZYNSKI',
            wind_mass_loss='WIND_ALGORITHM_BINARY_C_2020',
            sn_kick_dispersion_II=0,
            sn_kick_dispersion_IBC=0,
            sn_kick_dispersion_GRB_COLLAPSAR=0,
           )

search_results = search_for_BHBH(100, verbosity=VERBOSITY, **args)
if search_results:
    print(search_results['log_file_contents'])
Found bound BHBH system
      TIME      M1       M2   K1  K2           SEP        PER   ECC  R1/ROL1 R2/ROL2  TYPE RANDOM_SEED=3486 RANDOM_COUNT=0
     0.0000   60.000   40.000  1   1        3101.2     5.48 y  0.00   0.006   0.006  "INITIAL "
     0.0000   60.000   40.000  1   1        3101.2     5.48 y  0.00   0.006   0.006  "BEG_SYMB"
     3.8674   57.840   39.505  2   1        3185.7     5.78 y  0.00   0.017   0.010  "TYPE_CHNGE"
     3.8704   57.600   39.509  2   1        3192.2      5.8 y  0.00   0.393   0.010  "Start tidal lock 1"
     3.8704   57.592   39.509  2   1        3192.4      5.8 y  0.00   0.418   0.010  "End tidal lock 1"
     3.8706   57.568   39.510  4   1          3193     5.81 y  0.00   0.496   0.010  "TYPE_CHNGE"
     3.8910   54.510   39.758  4   1        3246.3     6.04 y  0.00   1.000   0.010  "BEG_RCHE 1>2"
     3.9650   42.495   42.571  4   1        3223.5     6.29 y  0.00   1.913   0.009  "q-inv"
     4.0585   24.961   50.453  4   1        3211.5     6.65 y  0.00   0.999   0.007  "END_RCHE 1!>2"
     4.0649   24.747   50.479  7   1        3216.1     6.67 y  0.00   0.002   0.007  "TYPE_CHNGE"
     4.2699   22.235   50.402 14   1        3331.3     7.15 y  0.00   0.000   0.007  Randbuf=32131 - Mers(0)=0.175268 - Mers(1)=0.383437 - Mers(2)=0.361343 - Mers(3)=0.457958 - Mers(4)=0.403947
     4.2699   22.235   50.402 14   1        3331.3     7.15 y  0.00   0.000   0.007  SN kick Ib/c (SN type 11 11, pre-explosion M=22.2347 Mc"CO"=16.8374 type=8) -> kick 1(0) vk=0 vr=64.5055 omega=2.53807 phi=-0.084183 -> vn=64.5055 ; final sep 3331.3 ecc 2.10734e-08 (random count 0) - Runaway v=(0,0,0) |v|=0 : companion v=(0,0,0), |v|=0 ;
     4.2699   22.235   50.402 14   1        3331.3     7.15 y  0.00   0.000   0.007  "TYPE_CHNGE"
     4.2699   22.235   50.402 14   1        3331.3     7.15 y  0.00   0.000   0.007  "SN"
     4.2799   22.235   50.398 14   1        3331.1     7.15 y  0.00   0.000   0.007  "END_SYMB"
     4.2946   22.235   50.392 14   1        3331.7     7.16 y  0.00   0.000   0.007  "BEG_BSS "
     5.9074   22.235   49.493 14   2        3373.9     7.34 y  0.00   0.000   0.013  "OFF_MS"
     5.9074   22.235   49.493 14   2        3373.9     7.34 y  0.00   0.000   0.013  "TYPE_CHNGE"
     5.9074   22.235   49.493 14   2        3373.9     7.34 y  0.00   0.000   0.013  "END_BSS"
     5.9113   22.235   49.243 14   4        3393.2     7.41 y  0.00   0.000   0.238  "TYPE_CHNGE"
     5.9192   22.236   48.059 14   4        3447.2     7.65 y  0.00   0.000   0.236  "Start tidal lock 2"
     5.9203   22.236   47.904 14   4        3454.5     7.69 y  0.00   0.000   0.236  "End tidal lock 2"
     5.9786   22.242   39.145 14   4        3893.9     9.83 y  0.00   0.000   1.001  "BEG_RCHE 2>1"
     6.0907   22.255   22.249 14   4          4523     14.5 y  0.00   0.000   1.296  "q-inv"
     6.1018   22.256   20.583 14   4        4540.5     14.8 y  0.00   0.000   0.997  "END_RCHE 2!>1"
     6.1137   22.256   20.232 14   7        4551.1     14.9 y  0.00   0.000   0.001  "TYPE_CHNGE"
     6.3373   22.256   19.925 14  14        4596.3     15.2 y  0.00   0.000   0.000  Mers(5)=0.101099 - Mers(6)=0.660012 - Mers(7)=0.331365 - Mers(8)=0.36634 - Mers(9)=0.576176
     6.3373   22.256   19.925 14  14        4596.3     15.2 y  0.00   0.000   0.000  SN kick Ib/c (SN type 11 11, pre-explosion M=19.9248 Mc"CO"=15.0519 type=8) -> kick 1(0) vk=0 vr=41.8488 omega=3.62022 phi=-0.27061 -> vn=41.8488 ; final sep 4596.27 ecc 2.58096e-08 (random count 5) - Runaway v=(0,0,0) |v|=0 : companion v=(0,0,0), |v|=0 ;
     6.3373   22.256   19.925 14  14        4596.3     15.2 y  0.00   0.000   0.000  "TYPE_CHNGE"
     6.3373   22.256   19.925 14  14        4596.3     15.2 y  0.00   0.000   0.000  "SN"
     6.3473   22.256   19.925 14  14        4596.3     15.2 y  0.00   0.000   0.000  "END_SYMB"
 15000.0000   22.256   19.925 14  14        4596.3     15.2 y  0.00   0.000   0.000  "MAX_TIME"

Oh dear! While the mass loss is reduced, there is now Roche-lobe overflow (RLOF)! This increases the secondary’s mass, but that means it also has RLOF. We can try making the system wider to prevent this, but we want a shorter period! Instead, let’s shorten the initial period to force common-envelope evolution. This will shrink the system. We will also turn off the mass loss to give us the best change of acquiring the system we want, and turn the SN kicks back on (they have less effect in closer, more grvitationally-bound systems).

[9]:
args = dict(
    M_1=60,
    M_2=40,
    metallicity=0.0001,
    orbital_period=2000, # days
    BH_prescription='BH_BELCZYNSKI',
    wind_mass_loss='WIND_ALGORITHM_NONE',
    alpha_ce = 0.1,
    lambda_ce = 0.5,
)


search_results = search_for_BHBH(100, verbosity=VERBOSITY, **args)
if search_results:
    print(search_results['log_file_contents'])
Found bound BHBH system
      TIME      M1       M2   K1  K2           SEP        PER   ECC  R1/ROL1 R2/ROL2  TYPE RANDOM_SEED=39800 RANDOM_COUNT=0
     0.0000   60.000   40.000  1   1        3101.2     5.48 y  0.00   0.006   0.005  "INITIAL "
     3.8571   60.000   40.000  2   1        3101.2     5.48 y  0.00   0.014   0.009  "TYPE_CHNGE"
     3.8600   60.000   40.000  4   1        3101.2     5.48 y  0.00   0.096   0.009  "TYPE_CHNGE"
     4.0451   60.000   40.000  4   1        3101.2     5.48 y  0.00   0.503   0.010  "Start tidal lock 1"
     4.0491   60.000   40.000  4   1        3101.2     5.48 y  0.00   0.561   0.010  "End tidal lock 1"
     4.0678   60.000   40.000  4   1        3099.3     5.47 y  0.00   1.000   0.010  "BEG_RCHE 1>2"
     4.0681   60.000   40.000  4   1        3099.1     5.47 y  0.00   1.009   0.010  "BEG_SYMB"
     4.2094   56.978   43.022  5   1          2884     4.91 y  0.00   1.672   0.010  "TYPE_CHNGE"
     4.2094   56.978   43.022  5   1          2884     4.91 y  0.00   1.672   0.010  Unstable RLOF : q=Mdonor/Maccretor=1.32438 > qc=0.733417 (prescription=float, donor [st=5 M=56.9778 Mc=25.8207 R=1945.85 Rc=1.67126], accretor [st=1 M=43.0222 Mc=0 R=10.1192 Rc=0] sep=2884.02 )
     4.2094   56.978   43.022  5   1          2884     4.91 y  0.00   1.672   0.010  COMENV n=0 presc=0 56.9778(Mc"He"=25.8207,"CO"=25.8207;EAGB)+43.0222(Mc;MS) - LAMBDA(fixed,st=5,m01=60,mc1=25.8207,l1=1.14773e+06,r1=1945.85,rzams=7.21541,convfrac=0.560301,lambda_ion=0.5)=0.5 - alpha=0.1 lambda=0.5 - a_in=2884.02 P_in=1793.89 Jtot_in=8.09752e+55 Jorb_in=7.95297e+55 EorbI=7.31037e+47 EbindI=-6.92616e+48 - sepF=30.1221 1:Rc=1.67126~RL=10.1116 2:Rc=10.1193~RL=12.7646 - : Mf1=25.8207(19.6094,9) Mf2=43.0222(0,1) af=30.1221 Jf=4.44174e+54 Jej=7.82682e+55 Eorbf=-6.99926e+49 dE=6.92616e+48 Eej=1.89371e+48 alpha_ej=0.0273415
     4.2095   25.821   43.022  9   1        30.122     2.31 d  0.00 647.092   0.793  post-COMENV M=25.8207(Mc"CO"=19.6094;HeGB;Menv=19.6094,Mstart=25.8207,Teff=0.407726,age=3102.24) + M=43.0222(Mc;MS;Menv=43.0222,Mstart=43.0222,Teff=3.62709,age=46857.8) a_out=30.1221 P_out=0.00631839 (2.30747 d) Jorb_out=4.44174e+54 R1=6543.21 RL1=10.1116 (R1/RL1 = 647.099) R2=10.1193 RL2=12.7646 (R2/RL2 = 0.792766) SN? 0 0 Single? No : merge in 1.60828 Gyr at t = 1.61249 Gyr, events pending? 1
     4.2095   25.821   43.022  9   1        30.122     2.31 d  0.00 647.092   0.793  "TYPE_CHNGE"
     4.2095   25.821   43.022  9   1        30.122     2.31 d  0.00 647.092   0.793  "q-inv"
     4.2095   24.803   44.040 14   1        38.566     3.34 d  0.32   0.000   0.457  Randbuf=79375 - Mers(0)=0.301545 - Mers(1)=0.0266946 - Mers(2)=0.404793 - Mers(3)=0.189909 - Mers(4)=0.0438563
     4.2095   24.803   44.040 14   1        38.566     3.34 d  0.32   0.000   0.457  SN kick Ib/c (SN type 11 11, pre-explosion M=24.803 Mc"CO"=19.6094 type=9) -> kick 1(190) vk=126.625 vr=509.248 omega=0.275557 phi=-0.668975 -> vn=421.918 ; final sep 38.5655 ecc 0.319412 (random count 0) - Runaway v=(0,0,0) |v|=0 : companion v=(0,0,0), |v|=0 ;
     4.2095   24.803   44.040 14   1        38.566     3.34 d  0.32   0.000   0.457  "TYPE_CHNGE"
     4.2095   24.803   44.040 14   1        38.566     3.34 d  0.32   0.000   0.457  "END_RCHE 1!>2"
     4.2095   24.803   44.040 14   1        38.566     3.34 d  0.32   0.000   0.457  "SN"
     4.2095   24.803   44.040 14   1        44.506     4.14 d  0.32   0.000   0.521  "END_SYMB"
     4.8539   24.803   44.040 14   1        46.203     4.38 d  0.15   0.000   0.629  "BEG_BSS "
     4.9016   24.803   44.040 14   1        45.848     4.33 d  0.13   0.000   0.648  "Start tidal lock 2"
     5.4816   24.803   44.040 14   1        44.633     4.16 d  0.00   0.000   0.954  "Circularized"
     5.5343   24.803   44.040 14   1        44.533     4.15 d  0.00   0.000   1.000  "BEG_RCHE 2>1"
     5.5348   24.803   44.040 14   1        44.533     4.15 d  0.00   0.000   1.001  "BEG_SYMB"
     5.5599   24.803   44.040 14   1        44.534     4.15 d  0.00   0.000   1.000  "END_RCHE 2!>1"
     5.5604   24.803   44.040 14   1        44.533     4.15 d  0.00   0.000   0.999  "END_SYMB"
     5.5867   24.803   44.040 14   2        44.899      4.2 d  0.00   0.000   0.781  "OFF_MS"
     5.5867   24.803   44.040 14   2        44.899      4.2 d  0.00   0.000   0.781  "TYPE_CHNGE"
     5.5867   24.803   44.040 14   2        44.899      4.2 d  0.00   0.000   0.781  "END_BSS"
     5.5867   24.803   44.040 14   2        44.899      4.2 d  0.00   0.000   0.781  "End tidal lock 2"
     5.5875   24.803   44.040 14   2        44.903      4.2 d  0.00   0.000   1.001  "BEG_RCHE 2>1"
     5.5878   24.803   44.040 14   2        44.904      4.2 d  0.00   0.000   1.106  "BEG_SYMB"
     5.5914   24.803   43.901 14   4        45.258     4.25 d  0.00   0.000   3.428  "TYPE_CHNGE"
     5.5916   24.804   43.873 14   4        45.461     4.28 d  0.00   0.000   3.413  "Start tidal lock 2"
     5.6245   24.808   39.603 14   4        48.625     4.89 d  0.00   0.000   3.245  "End tidal lock 2"
     5.7494   24.823   24.806 14   4        63.569     8.33 d  0.00   0.000   3.454  "q-inv"
     5.8032   24.830   17.333 14   4        75.087     11.6 d  0.00   0.000   0.998  "END_RCHE 2!>1"
     5.8033   24.830   17.333 14   4        75.087     11.6 d  0.00   0.000   0.990  "END_SYMB"
     5.8290   24.830   17.333 14   7        75.085     11.6 d  0.00   0.000   0.053  "TYPE_CHNGE"
     6.0403   24.830   17.333 14   8        75.085     11.6 d  0.00   0.000   0.050  "TYPE_CHNGE"
     6.0427   24.830   17.333 14  14        37.893     4.16 d  0.99   0.000   0.000  Mers(5)=0.552236 - Mers(6)=0.103282 - Mers(7)=0.393936 - Mers(8)=0.539534 - Mers(9)=0.864167
     6.0427   24.830   17.333 14  14        37.893     4.16 d  0.99   0.000   0.000  SN kick Ib/c (SN type 11 11, pre-explosion M=17.3329 Mc"CO"=13.0483 type=8) -> kick 1(190) vk=179.074 vr=327.351 omega=5.42972 phi=0.0791504 -> vn=44.4989 ; final sep 37.8925 ecc 0.991187 (random count 5) - Runaway v=(66.5305,-36.3905,-4.82105) |v|=75.9857 : companion v=(35.7333,-68.5019,-4.08275), |v|=77.3696 ;
     6.0427   24.830   17.333 14  14        37.893     4.16 d  0.99   0.000   0.000  "TYPE_CHNGE"
     6.0427   24.830   17.333 14  14        37.893     4.16 d  0.99   0.000   0.000  "SN"
     6.0742   24.830   17.333 14  14        0.1051   0.0146 h  0.06   0.002   0.002  Contact reached R/RL = 0.00319559 0.00262864, st 14 14
     6.0742   42.163    0.000 14  15            -1       -1   -1.00   0.000   0.000  Mers(10)=0.57683 - Mers(11)=0.0798097
     6.0742   42.163    0.000 14  15            -1       -1   -1.00   0.000   0.000  SN kick BH_BH (SN type 17 17, pre-explosion M=42.1626 Mc"CO"=42.1626,"BH"=42.1626 type=14) -> kick 0(0) vk=0 vr=0 omega=0.501459 phi=0.154271 -> vn=0 ; final sep 0 ecc -1 (random count 10) - Runaway v=(35.7333,-68.5019,-4.08275) |v|=77.3696 : companion v=(66.5305,-36.3905,-4.82105), |v|=75.9857 ;
     6.0742   42.163    0.000 14  15            -1       -1   -1.00   0.000   0.000  "TYPE_CHNGE"
     6.0742   42.163    0.000 14  15            -1       -1   -1.00   0.000   0.000  "COALESCE"
     6.0742   42.163    0.000 14  15            -1       -1   -1.00   0.000   0.000  "SN"
 15000.0000   42.163    0.000 14  15            -1       -1   -1.00   0.000   0.000  "MAX_TIME"

Note that this system has a far shorter period. Let’s run a number of these to see what we find.

Lets now try to find a system which at the formation of the final black hole has a separation of 10 \(R_{\odot}\) or less

[10]:
args = dict(
            M_1=60,
            M_2=40,
            metallicity=0.001,
            orbital_period=2000, # days
            BH_prescription='BH_BELCZYNSKI',
            wind_mass_loss='WIND_ALGORITHM_NONE',
            alpha_ce=0.1,
           )

found = False
while found == False:
    search_result = search_for_BHBH(100,**args)

    if search_result and search_result['system_properties']['separation'] < 100:
        found = True

        print(search_result['log_file_contents'])
Found bound BHBH system
Found bound BHBH system
      TIME      M1       M2   K1  K2           SEP        PER   ECC  R1/ROL1 R2/ROL2  TYPE RANDOM_SEED=13312 RANDOM_COUNT=0
     0.0000   60.000   40.000  1   1        3101.2     5.48 y  0.00   0.006   0.006  "INITIAL "
     3.7617   60.000   40.000  2   1        3101.2     5.48 y  0.00   0.018   0.010  "TYPE_CHNGE"
     3.7647   60.000   40.000  2   1        3101.2     5.48 y  0.00   0.572   0.010  "Start tidal lock 1"
     3.7647   60.000   40.000  4   1        3101.2     5.48 y  0.00   0.583   0.010  "TYPE_CHNGE"
     3.7721   60.000   40.000  4   1        3101.2     5.48 y  0.00   0.633   0.010  "End tidal lock 1"
     3.7785   60.000   40.000  4   1        3100.8     5.48 y  0.00   1.001   0.010  "BEG_RCHE 1>2"
     3.7785   60.000   40.000  4   1        3100.8     5.48 y  0.00   1.001   0.010  "BEG_SYMB"
     4.1326   52.305   47.695  5   1        2791.6     4.68 y  0.00   2.159   0.011  "TYPE_CHNGE"
     4.1326   52.305   47.695  5   1        2791.6     4.68 y  0.00   2.159   0.011  Unstable RLOF : q=Mdonor/Maccretor=1.09665 > qc=0.714712 (prescription=float, donor [st=5 M=52.3048 Mc=27.0276 R=2332.11 Rc=1.71818], accretor [st=1 M=47.6952 Mc=0 R=10.9757 Rc=0] sep=2791.57 )
     4.1326   52.305   47.695  5   1        2791.6     4.68 y  0.00   2.159   0.011  COMENV n=0 presc=0 52.3048(Mc"He"=27.0276,"CO"=27.0276;EAGB)+47.6952(Mc;MS) - LAMBDA(fixed,st=5,m01=60,mc1=27.0276,l1=1.44223e+06,r1=2332.11,rzams=8.33552,convfrac=0.467851,lambda_ion=0.5)=0.5 - alpha=0.1 lambda=0.5 - a_in=2791.57 P_in=1708.33 Jtot_in=8.09706e+55 Jorb_in=7.96557e+55 EorbI=8.76416e+47 EbindI=-4.30388e+48 - sepF=55.7114 1:Rc=1.71818~RL=18.4385 2:Rc=10.9758~RL=23.8926 - : Mf1=27.0276(20.5423,9) Mf2=47.6952(0,1) af=55.7114 Jf=6.72832e+54 Jej=7.63566e+55 Eorbf=-4.39152e+49 dE=4.30388e+48 Eej=1.03996e+48 alpha_ej=0.0241634
     4.1327   27.028   47.695  9   1        55.711     5.57 d  0.00 367.232   0.459  post-COMENV M=27.0276(Mc"CO"=20.5423;HeGB;Menv=20.5423,Mstart=27.0276,Teff=0.397754,age=3084.58) + M=47.6952(Mc;MS;Menv=47.6952,Mstart=47.6952,Teff=2.9071,age=46089.5) a_out=55.7114 P_out=0.0152545 (5.57093 d) Jorb_out=6.72832e+54 R1=6771.28 RL1=18.4385 (R1/RL1 = 367.236) R2=10.9758 RL2=23.8926 (R2/RL2 = 0.459382) SN? 0 0 Single? No : merge in 14.9411 Gyr at t = 14.9452 Gyr, events pending? 1
     4.1327   27.028   47.695  9   1        55.711     5.57 d  0.00 367.232   0.459  "TYPE_CHNGE"
     4.1327   27.028   47.695  9   1        55.711     5.57 d  0.00 367.232   0.459  "q-inv"
     4.1327   26.215   48.508 14   1        68.001     7.51 d  0.76   0.000   0.284  Randbuf=67260 - Mers(0)=0.643477 - Mers(1)=0.441321 - Mers(2)=0.365945 - Mers(3)=0.352254 - Mers(4)=0.0493316
     4.1327   26.215   48.508 14   1        68.001     7.51 d  0.76   0.000   0.284  SN kick Ib/c (SN type 11 11, pre-explosion M=26.2147 Mc"CO"=20.5423 type=9) -> kick 1(190) vk=265.461 vr=401.272 omega=0.309959 phi=-0.299971 -> vn=335.175 ; final sep 68.0005 ecc 0.759863 (random count 0) - Runaway v=(236.986,186.272,-13.8186) |v|=301.746 : companion v=(8.55202,-58.8609,-5.46244), |v|=59.7292 ;
     4.1327   26.215   48.508 14   1        68.001     7.51 d  0.76   0.000   0.284  "TYPE_CHNGE"
     4.1327   26.215   48.508 14   1        68.001     7.51 d  0.76   0.000   0.284  "END_RCHE 1!>2"
     4.1327   26.215   48.508 14   1        68.001     7.51 d  0.76   0.000   0.284  "SN"
     4.1327   26.215   48.508 14   1         78.19     9.26 d  0.76   0.000   0.322  "END_SYMB"
     4.3844   26.215   48.508 14   1        38.009     3.14 d  0.12   0.000   0.711  "Start tidal lock 2"
     4.4261   26.215   48.508 14   1        37.629     3.09 d  0.05   0.000   0.727  "BEG_BSS "
     4.8255   26.215   48.508 14   1        37.339     3.06 d  0.00   0.000   0.848  "Circularized"
     5.1375   26.215   48.508 14   1        37.026     3.02 d  0.00   0.000   1.000  "BEG_RCHE 2>1"
     5.1380   26.215   48.508 14   1        37.026     3.02 d  0.00   0.000   1.001  "BEG_SYMB"
     5.7437   26.276   47.488 14   2        37.702     3.12 d  0.00   0.000   1.185  "OFF_MS"
     5.7437   26.276   47.488 14   2        37.702     3.12 d  0.00   0.000   1.185  "TYPE_CHNGE"
     5.7437   26.276   47.488 14   2        37.702     3.12 d  0.00   0.000   1.185  "END_BSS"
     5.7437   26.276   47.488 14   2        37.702     3.12 d  0.00   0.000   1.185  "End tidal lock 2"
     5.7464   26.276   47.110 14   2        38.756     3.26 d  0.00   0.000   6.815  "Start tidal lock 2"
     5.7466   26.276   46.978 14   2        38.938     3.29 d  0.00   0.000   7.881  "End tidal lock 2"
     5.7479   26.276   45.466 14   4        40.188     3.48 d  0.00   0.000  18.114  "TYPE_CHNGE"
     5.7598   26.278   26.248 14   4        54.754     6.47 d  0.00   0.000  14.114  "q-inv"
     5.7691   26.279   18.231 14   4        64.697     9.03 d  0.00   0.000   0.995  "END_RCHE 2!>1"
     5.7692   26.279   18.231 14   4        64.697     9.03 d  0.00   0.000   0.987  "END_SYMB"
     5.7854   26.279   18.231 14   7        64.696     9.03 d  0.00   0.000   0.061  "TYPE_CHNGE"
     5.9061   26.279   18.231 14   7        64.695     9.03 d  0.00   0.000   0.063  "Start tidal lock 2"
     6.1521   26.279   18.231 14   7        64.695     9.03 d  0.00   0.000   0.063  "End tidal lock 2"
     6.1942   26.279   18.231 14   8        64.695     9.03 d  0.00   0.000   0.060  "TYPE_CHNGE"
     6.1949   26.279   18.231 14  14        49.936     6.13 d  0.30   0.000   0.000  Mers(5)=0.534629 - Mers(6)=0.163092 - Mers(7)=0.907792 - Mers(8)=0.0909357 - Mers(9)=0.00121347
     6.1949   26.279   18.231 14  14        49.936     6.13 d  0.30   0.000   0.000  SN kick Ib/c (SN type 11 11, pre-explosion M=18.2315 Mc"CO"=13.7429 type=8) -> kick 1(190) vk=140.491 vr=362.342 omega=0.00762447 phi=-0.958149 -> vn=304.115 ; final sep 49.9359 ecc 0.295576 (random count 5) - Runaway v=(8.55202,-58.8609,-5.46244) |v|=59.7292 : companion v=(236.986,186.272,-13.8186), |v|=301.746 ;
     6.1949   26.279   18.231 14  14        49.936     6.13 d  0.30   0.000   0.000  "TYPE_CHNGE"
     6.1949   26.279   18.231 14  14        49.936     6.13 d  0.30   0.000   0.000  "SN"
 15000.0000   26.279   18.231 14  14         41.31     4.61 d  0.23   0.000   0.000  "MAX_TIME"

Playing around with the input for an individual system gives us some understanding of which parameters affect the formation of black hole systems, but, as mentioned in the introduction, this is nto the most appropriate way to find many such systems. We want to make use of the Population utilities, and set up custom logging and a parse function.

Things to try next: * Set up a population object with a grid in M1, q and orbital period (10 x 10 x 10) * Set up a custom logging function to catch BHBH systems, and include the ZAMS properties in the output. * Set up a parse function to process the output, try writing the output to a file. * Run this population at several different metallicities. Do you see a trend?