<jupyter_start><jupyter_text>**How Bright are Accreting Black Holes in Binary Systems?** ***A Computational Investigation***Leader: <NAME>  **Background**Binary star systems consist of two stars that orbit one another. Depending on the mass of these stars, one or more of them may explode in a supernova and form a black hole late in their lifetime. This is generally the case for stars with initial masses greater than 20 times that of the sun, or 20 solar masses ($M_\odot$). In the majority of these binary systems, black hole formation occurs asynchronously. This implies that the binary undergoes a period time where one object is a black hole and one object is still a star. This stage of evolution is called the x-ray binary phase, and the system is referred to as an **x-ray binary**.Why do we call these systems x-ray binaries? During this evolutionary stage, the black hole accretes matter from its companion star due to its strong gravitational pull. The accretion process heats transferred material through gravitational and frictional forces, which in turn triggers the emission of x-rays. Below is a helpful illustration of this process:![photo1](https://drive.google.com/uc?export=view&id=1ya2mtVhdkbPInwV-kysQbDaqj-DYc--v) Because the matter surrounding these black holes emits brightly in the x-ray band of the electromagnetic spectrum, x-ray observatories (such as the Chandra Observatory) search for this emission to identify binary systems that house black holes. The brightness, or luminosity, that astronomers observe in these binaries is called the **accretion luminosity**, as it is the luminosity that results from the accretion process of stellar material onto the black hole. These observations help astronomers learn about the properties of these systems and predict their possible courses of evolution. For example, some astronomers are interested in finding systems like these that will eventually form colliding black holes, which can be observed with gravitational wave detections.In this project, we will take a closer look at how we can predict the luminosity of these black hole systems using computational methods and data from simulations.  Theoretical calculations like these help astronomers target their observations towards certain luminosities so that they can better identify these systems in the sky. These calculations also help astronomers compare their theoretical predictions with observations so that they can better understand if there are flaws in their theoretical or observational methods.  **Building a Luminosity Function**First, we need to write a function that calculates and returns the accretion luminosity of a black hole in a binary system. Once we have this function, we can apply it to simulation data in multiple different ways. Previous theoretical work tells us that the black hole accretion luminosity is calculated by: $L_{acc} = \eta \dot{M}_{acc}c^2$where $L_{acc}$ is the accretion luminosity, $\eta$ is a scale factor derived from relativistic properties of the black hole, $\dot{M}_{acc}$ is the accretion rate of material onto the black hole, and $c$ is the speed of light. The parameters $\eta$ and $\dot{M}_{acc}$ are calculated by: $\eta = 1 - \sqrt{1 - \bigg(\frac{M_{BH}}{3M_{BH,0}}\bigg)^2}$ $\dot{M}_{acc} = \frac{\dot{M}_{BH}}{1 - \eta}$where $M_{BH}$ is the current black hole mass, $M_{BH,0}$ is the initial black hole mass, and $\dot{M}_{BH}$ is the rate of change of the black hole mass. The parameters provided in our simulation output are $M_{BH}$, $M_{BH,0}$, and $\dot{M}_{BH}$. Using this information and the equations above, define a function in python that takes the appropriate input parameters, calculates $L_{acc}$, $\eta$, and $\dot{M}_{acc}$ in the proper order, and returns the accretion luminosity. Don't worry about defining numeric quantities for constants yet.<jupyter_code><jupyter_output><empty_output><jupyter_text>There are a couple refinements we need to make to our function before we can apply it to our simulation data. The first is that there is an upper limit for $\eta$. Specifically, if $M_{BH}$ exceeds $\sqrt{6}M_{BH,0}$, then $\eta$ is taken to be $0.42$. Copy and paste your old function below, and then add this condition into your calculation for $\eta$.<jupyter_code><jupyter_output><empty_output><jupyter_text>Next, we need to apply a **bolometric correction** to our calculated luminosity. A bolometric correction is a multiplicative factor that accounts for differences between the visible luminosity that we observe and the true brightness of the observed source. Bolometric corrections usually vary depending on the observing apparatus and specific properties of the observed object. However, since we are simulating binaries and not observing them, we take the bolometric correction factor to be constant for all sources. For binaries with accreting black holes, this correction is $0.8$. Copy and paste your old function below, and then multiply your final accretion luminosity by this factor before returning it.<jupyter_code><jupyter_output><empty_output><jupyter_text>Lastly, we need to adjust the units of our calculation so that our final accretion luminosity has the traditional units of x-ray observations. To do this, we need to know the units of our simulation variables. Below is a list of the units for our simulation output data:$M_{BH}$ ---> [$M_\odot$]$M_{BH,0}$ ---> [$M_\odot$]$\dot{M}_{BH}$ ---> [$M_\odot yr^{-1}$]We want our final luminosity to be in **ergs per second**, or $erg$  $s^{-1}$, which is a common unit of measure for luminosities in astronomy. A helpful unit conversion for one $erg$ is:[$erg$] = [$g$ $cm^2$ $s^{-2}$] Unforunately, it looks like we have a lot of converting to do! Alas, if you continue with astronomy or astrophysics, you'll find this is a common practice, as astronomers like to use a plethora of unique and often inconvenient units.  Using dimensional analysis, work out any additional multiplicative constants you may need to add to your calculations so that your returned accretion luminosity has the correct units. You can choose the units of $c$ to be whatever is most convenient for your calculations. Copy and paste your old function below, and then write in your unit adjustments. Once you are done, define numeric quantities for your constants at the beginning of your function.<jupyter_code><jupyter_output><empty_output><jupyter_text>Your function is now ready to go! Let's put it to work with some simulation data.  **Applying Your Luminosity Function to Data: ```for``` Loop Method** First, read-in the simulation data file called ```"cosmic_data_1.csv"```, located in this project folder, and store it as a ```pandas``` dataframe. Make sure to import ```pandas``` and any other packages you might need first!<jupyter_code>import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from google.colab import files
uploaded = files.upload()

df1 = pd.read_csv("cosmic_data_1.csv")<jupyter_output><empty_output><jupyter_text>Now, open the file and examine its format, including column titles, the number of rows, etc.<jupyter_code><jupyter_output><empty_output><jupyter_text>This file contains the evolutionary data for a single simulated binary system during its x-ray binary phase. Each row represents one timestep of the simulation, which you can see by looking at the ```tphys``` column. ```tphys``` is the time, in mega years (Myr), since the binary began its evolution. As a warm-up to practice working with the data, use ```pandas``` operations to calculate how long this system was in the x-ray binary phase, and print your result.<jupyter_code>time_xrb = df1['tphys'].iloc[-1] - df1['tphys'].iloc[0]
print (time_xrb, "Myr")<jupyter_output>1.309999999999972 Myr<jupyter_text>The parameters we will need to calculate the accretion luminosity are ```mass_2```, which represents the black hole mass, and ```deltam_2```, which represents the rate of change in the black hole mass. One way we can acquire luminosity information is to loop through the dataframe, row by row, and calculate the luminosity at each time step using our function.First, write a ```for``` loop that goes through each row of the data frame and applies your luminosity function at each iteration. Save your calculated luminosities in an array. Print your final array of luminosities to check that they are reasonable. Reasonable luminosities should fall between $10^{30}$ and $10^{45}$ $erg$ $s^{-1}$. *Hint: you will need to use a special iterator for your dataframe. Try the pandas function iterrows().*<jupyter_code>def calc_luminosity_loop(current_BH_mass, initial_BH_mass, mdot_BH):

  ## phsyical constants                                                                                                                                   
  c = 2.99e10         ## speed of light in cm/s                                                                                                           
  secyr = 3.154e7    ## seconds per year                                                                                                                  
  Myr = 1e6        ## years per Myr                                                                                                                       
  Msun = 1.989e33    ## grams per solar mass    

  bolometric_correction = 0.8

  if current_BH_mass > np.sqrt(6)*initial_BH_mass:
    eta = 0.42
  else:
    eta = 1 - np.sqrt(1-(current_BH_mass/(3*initial_BH_mass))**2)

  acc_rate = mdot_BH/(1-eta)
  luminosity = bolometric_correction*eta*acc_rate*c**2*Msun/secyr 
  
  return luminosity

acc_lum_loop = []

M_bh_0 = df1['mass_2'].iloc[0]

for index, row in df1.iterrows():

  M_bh = row['mass_2']
  Mdot_bh = row['deltam_2']

  luminosity = calc_luminosity_loop(M_bh, M_bh_0, Mdot_bh)
  acc_lum_loop.append(luminosity)

print (acc_lum_loop)<jupyter_output>[1.1253823492823668e+35, 1.1408901655388753e+35, 1.1566749933938256e+35, 1.1727424293171707e+35, 1.189098170691318e+35, 1.205748015670417e+35, 1.2226978628376603e+35, 1.2399537106417342e+35, 1.257521656592063e+35, 1.2754078961913838e+35, 1.293618721582208e+35, 1.3121605198824428e+35, 1.3310397711835939e+35, 1.350263046183022e+35, 1.369837003420074e+35, 1.3897683860834096e+35, 1.4100640183552084e+35, 1.4307308012550807e+35, 1.4517757079446097e+35, 1.4732057784503908e+35, 1.4950281137611946e+35, 1.5172498692516261e+35, 1.539878247382079e+35, 1.5629204896211678e+35, 1.5863838675340488e+35, 1.610275672976188e+35, 1.6346032073286615e+35, 1.6593737697074791e+35, 1.6845946440752246e+35, 1.7102730851793114e+35, 1.7364163032371813e+35, 1.763031447283878e+35, 1.7901255870933556e+35, 1.81770569358e+35, 1.8457786175819872e+35, 1.874351066923311e+35, 1.9034295816462324e+35, 1.9330205073007557e+35, 1.9631299661726803e+35, 1.9937638263264203e+35, 2.024927668333835e+35, 2.0566267495548505e+35, 2.08886[...]<jupyter_text>Now that we have the accretion luminosity data for our binary, let's plot it to see how the system's emission changes over time. Using ```matplotlib```, make a time series plot using the ```tphys``` column from the simulation output along with your array of luminosities. Make sure to properly label your plot! *Hint: Use a log scale when plotting your luminosities so that you can better visualize the emission behavior.*<jupyter_code>fig, ax = plt.subplots(figsize=(8,6))
plt.plot(df1['tphys'], np.log10(acc_lum_loop))
plt.xlabel("Time (Myr)", size=13)
plt.ylabel("Accretion Luminosity (log$_{10}$ erg/s)", size=13)
plt.title("Accretion Luminosity vs. Time", size=14)
plt.show()<jupyter_output><empty_output><jupyter_text>One question researchers might ask is if a given binary system is observable. For a system to be observed by an x-ray observatory, it must reach a maximum luminosity above an observable threshold. Let's take our observable threshold to be $10^{37} erg$ $s^{-1}$, which is a common threshold used for the Chandra Observatory. Using your calculations, find the maximum luminosity of your binary system. Is it observable?<jupyter_code>print (np.max(acc_lum_loop), "erg/s")<jupyter_output>1.0352989720162052e+39 erg/s<jupyter_text>**Applying Your Luminosity Function to Data: Array Method** It turns out that ```for``` loops are too computationally slow to analyze very large files or many files in one run. Luckily, we can speed-up our code by passing arrays to our luminosity function instead of passing values at each timestep of the data. This will allow us to calculate all luminosities for the binary in one step! However, this method will require some alteration to your luminosity function.Primarily, now that your function takes arrays as parameters instead of values, you can no longer use conditional logic to determine what $\eta$ should be. Instead, you should use boolean array logic, which is a standard application of ```numpy``` arrays.  For example, let's say I have two arrays of the same size called ```x``` and ```y```. If want to know which values of ```x``` are less than their corresponding values in ```y```, I can write:```x_lower_values = x < y``````x_lower_values``` is an array of the same size as ```x``` and ```y```. It it contains ```True/False``` entries that tell me if the corresponding entry in ```x``` is less than that in ```y```. Another useful tool that might help you rework your luminosity function is ```np.where()```. This function takes a boolean array (such as ```x_lower_values``` in the example above) and returns an array that yields entires with two possible outcomes. The outcomes are supplied by you and depend on if a given boolean array entry is ```True``` or ```False```. Copy and paste your old luminosity function below. Then, using the tools described above and/or any other methods you prefer, alter you luminosity function so that it takes array parameters and then calculates and returns an array of accretion luminosities using array operations.<jupyter_code>def calc_luminosity_arrays(current_BH_mass, initial_BH_mass, mdot_BH):

    ## phsyical constants                                                                                                                                   
    c = 2.99e10         ## speed of light in cm/s                                                                                                           
    secyr = 3.154e7    ## seconds per year                                                                                                                  
    Myr = 1e6        ## years per Myr                                                                                                                       
    Msun = 1.989e33    ## grams per solar mass 

    bolometric_correction = 0.8
    where_lower_masses = current_BH_mass < np.sqrt(6)*initial_BH_mass

    eta_lower_masses = 1 - np.sqrt(1-(current_BH_mass/(3*initial_BH_mass))**2)
    eta = np.where(where_lower_masses, eta_lower_masses, 0.42)

    acc_rate = mdot_BH/(1-eta)
    luminosity = bolometric_correction*eta*acc_rate*c**2*Msun/secyr   ## accretion luminosity in erg/sec                                                                                       

    return luminosity<jupyter_output><empty_output><jupyter_text>Now, apply your new and improved function to the data set above. You should now only need one line of code to calculate your array of accretion luminosities. Print your final array of luminosities to check that they are reasonable. *Hint: you may need to create an array for $M_{BH,0}$ that is the same length as your arrays for $M_{BH}$ and $\dot M_{BH}$ in order to pass it to your function.*<jupyter_code>M_bh_0 = np.ones(len(df1['mass_2']))*df1['mass_2'].iloc[0]

acc_lum_array = calc_luminosity_arrays(df1['mass_2'], M_bh_0, df1['deltam_2'])

print (acc_lum_array)<jupyter_output>0      1.125382e+35
1      1.140890e+35
2      1.156675e+35
3      1.172742e+35
4      1.189098e+35
           ...     
127    1.032864e+39
128    1.029047e+39
129    1.023241e+39
130    1.015235e+39
131    1.004801e+39
Name: deltam_2, Length: 132, dtype: float64<jupyter_text>To double check that your improved function is working correctly, make another time series plot using ```matplotlib```. This time, plot the accretion luminosities you obtained with the ```for``` loop method as well as those you obtained with the array method. They should appear to be the same line on your plot! Make sure to properly label your plot.<jupyter_code>fig, ax = plt.subplots(figsize=(8,6))
plt.plot(df1['tphys'], np.log10(acc_lum_loop))
plt.plot(df1['tphys'], np.log10(acc_lum_array))
plt.xlabel("Time (Myr)", size=13)
plt.ylabel("Accretion Luminosity (log$_{10}$ erg/s)", size=13)
plt.title("Accretion Luminosity vs. Time", size=14)
plt.show()<jupyter_output><empty_output><jupyter_text>One fun way to examine the computational difference between your two luminosity functions is to time how long each of them takes to run. Copy and paste your pieces of code that run your different luminosity functions below. Import the python ```time``` module and use the ```time()``` function to determine the runtime of each code segment. Print both runtimes to compare them.<jupyter_code>import time

start = time.time()
M_bh_0 = np.ones(len(df1['mass_2']))*df1['mass_2'].iloc[0]
acc_lum_array = calc_luminosity_arrays(df1['mass_2'], M_bh_0, df1['deltam_2'])
end = time.time()
array_time = end-start
print ("Arrays: ", array_time, "sec")


start = time.time()
acc_lum = []
M_bh_0 = df1['mass_2'].iloc[0]

for index, row in df1.iterrows():

  M_bh = row['mass_2']
  Mdot_bh = row['deltam_2']

  luminosity = calc_luminosity_loop(M_bh, M_bh_0, Mdot_bh)
  acc_lum.append(luminosity)
end = time.time()
loop_time = end-start
print ("For Loop: ", loop_time, "sec")<jupyter_output>Arrays:  0.002843141555786133 sec
For Loop:  0.014995574951171875 sec<jupyter_text>Even though the difference in runtimes may seem small, it is actually quite noticeable when applied to many systems. For example, one standard simulation run contains a population of ~50,000 binaries. Calculate how long each method would take to analyze the luminosities of a full population of binaries, and convert the resulting times (in seconds) to more useful units.<jupyter_code>print ("Arrays: ", 50000*array_time/60, "min")
print ("For Loop: ", 50000*loop_time/60, "min")<jupyter_output>Arrays:  2.3692846298217773 min
For Loop:  12.496312459309896 min<jupyter_text>Now, suppose we run a total of 10 different simulations, and need to analyze the luminosities of all systems in each simulation. How long does each method take to complete these calculations? Were the improvements to your luminosity function worth it?<jupyter_code>print ("Arrays: ", 10*50000*array_time/60, "min")
print ("For Loop: ", 10*50000*loop_time/60/60, "hrs")<jupyter_output>Arrays:  23.692846298217773 min
For Loop:  2.082718743218316 hrs<jupyter_text>**Bonus: Extended Data File** The file ```cosmic_data_1.csv``` is actually a very small portion of a full simulation output file. Normally, researchers want information about the entire evolution of a binary system so that characteristics of its behavior, such as its merge status, emission properties, etc. are placed in context with one other. This allows correlations to be drawn between evolutionary behaviors and facilitates a more complete understanding of how binary evolution works.Read-in the simulation data file ```cosmic_data_2.csv```, also located in this project folder, and store it as a ```pandas``` dataframe (it may take a couple of minutes to upload). Open the file and examine its format, including column titles, the number of rows, etc.<jupyter_code>from google.colab import files
uploaded = files.upload()

import pandas as pd
import numpy as np

df2 = pd.read_csv("cosmic_data_2.csv")<jupyter_output><empty_output><jupyter_text>This file is clearly far more massive than the previous one - and it is still a reduced version of the full simulation output! To work with the complete data set, we would need a supercomputer. However, this reduction will be suitable for the remainder of our investigation. The bonus portion of this project will involve identifying the x-ray binary phase of the system ourselves and applying our luminosity function to the appropriate subset of the extended data file. First, we need to identify when the x-ray binary phase begins. We can do this by finding the timestep when the first black hole is formed, which is recorded in the ```kstar_1``` and ```kstar_2``` columns of the data file. If ```kstar_1``` or ```kstar_2``` equals $14$, then that respective object is a black hole. The x-ray binary phase begins when the first black hole's accretion rate becomes greater than zero.Write a piece of code that identifies the timestep or index of the data file at which each object in the binary becomes a black hole. Store each index as a variable.<jupyter_code>BH1_index = np.where(df2['kstar_1'] == 14)[0][0]
BH2_index = np.where(df2['kstar_2'] == 14)[0][0]<jupyter_output><empty_output><jupyter_text>Using conditional logic, compare the indicies or timesteps you found above, and determine which one occurred first. Then, using this information, identify the data columns for $M_{BH}$ and $\dot{M}_{BH}$ that correspond to the first black hole object. Also store the first and second black hole formation indicies as separate variables.<jupyter_code>if BH2_index > BH1_index: 
  first_BH = BH1_index
  second_BH = BH2_index
  BH_obj = "kstar_1"
  BH_mass = "mass_1"
  BH_mdot = "deltam_1"

else: 
  first_BH = BH2_index
  second_BH = BH1_index
  BH_obj = "kstar_2"
  BH_mass = "mass_2"
  BH_mdot = "deltam_2"<jupyter_output><empty_output><jupyter_text>Now, using the appropriate data columns, identify the index when the black hole accretion rate becomes positive. This index signifies the beginning of the x-ray binary phase. Store this index as a separate variable.<jupyter_code>XRB_start = np.where(df2[BH_mdot] > 0)[0][0]<jupyter_output><empty_output><jupyter_text>The x-ray binary phase terminates when the black hole accretion rate is no longer positive. Using similar methods as above, identify the index when the x-ray binary phase ends. Store this index as a separate variable. To check that you indicies make sense, print the indicies of black hole formation as well as those identifying the x-ray binary phase.<jupyter_code>XRB_end = np.where(df2[BH_mdot] > 0)[0][-1]

print ("XRB phase begins: ", XRB_start) 
print ("XRB phase ends: ", XRB_end) 
print ("First BH forms: ", first_BH)
print ("Second BH forms: ", second_BH)<jupyter_output>XRB phase begins:  12
XRB phase ends:  127
First BH forms:  12
Second BH forms:  128<jupyter_text>Cut the data you need to calculate the accretion luminosity using the indices you found above so that the data begins and ends with the x-ray binary phase. Store your reduced data in separate arrays.<jupyter_code>M_bh_0_value = df2[BH_mass][first_BH]
Mbh = df2[BH_mass][XRB_start:XRB_end]
Mdot_bh = df2[BH_mdot][XRB_start:XRB_end]<jupyter_output><empty_output><jupyter_text>You should now be able to pass your reduced data arrays directly to your function! Do this to calculate the accretion luminosity of the binary during the x-ray binary phase, and print the resulting array of luminosities. Determine if this system is observable.<jupyter_code>M_bh_0 = np.ones(len(Mbh))*M_bh_0_value
acc_lum_extended = calc_luminosity_arrays(Mbh, M_bh_0, Mdot_bh)

print (acc_lum_extended, "\n")
print ("Maximum luminosity: ", np.max(acc_lum_extended))<jupyter_output>12     1.191027e+30
13     1.209642e+30
14     1.227897e+30
15     1.246501e+30
16     1.265461e+30
           ...     
122    9.023618e+38
123    9.013496e+38
124    8.988866e+38
125    8.948127e+38
126    8.889543e+38
Name: deltam_2, Length: 115, dtype: float64 

Maximum luminosity:  9.023618244430517e+38<jupyter_text>As a last step, compare your new luminosities to those you found above with the ```cosmic_data_1.csv``` data set. Use ```matplotlib``` to plot a time series of the binary emission over time, and label your plot accordingly. *Hint: you will need to apply the same cuts to the ```tphys``` column of ```cosmic_data_2.csv``` that you applied to your luminosity input parameters in order to plot your luminosities.*What are some similarities and differences between the emission of the first and second binary system?<jupyter_code>time = df2['tphys'][XRB_start:XRB_end]

fig, ax = plt.subplots(figsize=(8,6))
plt.plot(df1['tphys'], np.log10(acc_lum_array))
plt.plot(time, np.log10(acc_lum_extended))
plt.xlabel("Time (Myr)", size=13)
plt.ylabel("Accretion Luminosity (log$_{10}$ erg/s)", size=13)
plt.title("Accretion Luminosity vs. Time", size=14)
plt.show()<jupyter_output><empty_output>
StarcoderdataJupyterStructuredCleanDedup