tangos_bhs_intro.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "db454c01-d072-4109-b4e7-5497a7718eac",
   "metadata": {},
   "source": [
    "<h1> Black Holes in the Database <\\h1>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dd90c26e-e875-47a5-b31b-e298c17b17b4",
   "metadata": {},
   "source": [
    "Many galaxies/halos have supermassive black holes associated with them in the simulation. Here we will go over some of the basic ways we can access information about black holes and relate them back to their host galaxies/halos (and vice versa)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a76760ba-1c85-47ef-a4ed-4492f95fe1aa",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "#start by importing some needed modules\n",
    "#You may see some warnings... unless it spits out an error, you can ignore them!\n",
    "import tangos as db\n",
    "import numpy as np\n",
    "import pynbody\n",
    "import matplotlib.pyplot as plt\n",
    "plt.ion()\n",
    "\n",
    "#Set up some properties to make the plots look nice\n",
    "plt.rcParams['figure.figsize'] = (12, 9) #set figure size\n",
    "plt.rcParams['font.size'] = 20 #set font size so things are readible and not too small or bit\n",
    "plt.rcParams['font.family'] = 'sans-serif' # just a font preference here\n",
    "plt.rcParams['xtick.major.size'] = 10 #size and widths of major/minor tick marks on plots\n",
    "plt.rcParams['xtick.major.width'] = 2\n",
    "plt.rcParams['ytick.major.size'] = 10\n",
    "plt.rcParams['ytick.major.width'] = 2\n",
    "plt.rcParams['xtick.minor.size'] = 5\n",
    "plt.rcParams['xtick.minor.width'] = 1\n",
    "plt.rcParams['ytick.minor.size'] = 5\n",
    "plt.rcParams['ytick.minor.width'] = 1\n",
    "\n",
    "\n",
    "#load in the database simulation (also a good check that you are reading the right database file!)\n",
    "sim = db.get_simulation('cosmo25')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "83a20c89-c3ce-4512-b8e6-74959d1a46e1",
   "metadata": {},
   "source": [
    "<h3> Finding Black Holes in the Database </h3>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "69207c78-9a53-4526-a19b-0b2803c59c95",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# many galaxies/halos have supermassive black holes associated with them in the simulation\n",
    "# If you look at the available properties for a halo, you may see several \"BH\" and \"BH_central\" properties pop up\n",
    "\n",
    "halo = db.get_halo('cosmo25/%7779/1') # get halo 1, the most massive halo in the simulation at timestep number 7779\n",
    "halo.keys()                           #prints out all of the available properties"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "887499d3-a917-44b2-ad47-82f076dae8d6",
   "metadata": {},
   "outputs": [],
   "source": [
    "# That's a lot of black holes! Turns out, this was one of our major results from Romulus. \n",
    "# More massive halos/galaxies host a ton of \"wandering\" supermassive black holes. The most massive halo in this simulation, therefore, hosts the most.\n",
    "# But, even smaller galaxies can host multiple black holes. Let's look at a lower mass galaxy more similar to the Milky Way\n",
    "\n",
    "halo = db.get_halo('cosmo25/%7779/44') #get halo 44 from timestep number 7779\n",
    "\n",
    "# We look at all of the black holes in the halo by asking for the property \"BH\"\n",
    "\n",
    "print(halo['BH'])\n",
    "\n",
    "#Note that this property is a LIST. More specifically, it is a list of black hole objects.\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2eefa4f1-bc03-4787-bdb1-0ab6c80e2c19",
   "metadata": {},
   "outputs": [],
   "source": [
    "# If we want to examine only one BH, we can provide an index\n",
    "\n",
    "print(halo['BH'][0]) #print out the first black hole object associated with halo number 48\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4a05f3ba-a037-42bb-a014-1e6308d7ef93",
   "metadata": {},
   "outputs": [],
   "source": [
    "# You may notice there are two kinds of BHs associated with a galaxy/halo. \n",
    "# \"BH\" list ALL the black holes within the entire dark matter halo (i.e. the BHs in the galaxy and anything gravitationally bound to it)\n",
    "# The \"BH_central\" property excludes any BH that is associted with SUBSTRUCTURE, or halos/galaxies that are orbiting the galaxy we care about.\n",
    "# This will often be the BHs that we care about. These are the ones associated with our target galaxy, rather than some other galaxy orbiting around it.\n",
    "\n",
    "print(halo['BH_central'])\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0e17af9f-f577-4581-846b-896fbdea68c9",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Because of how they are defined, BH_central is a SUBSET of of BH (every object in \"BH_central\" is also bound to the halo, so is also included in \"BH\")\n",
    "\n",
    "print(len(halo['BH']), len(halo['BH_central'])) #if we check their numbers with len(), indeed there are fewer \"BH_central\" objects.\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b257459f-f44f-43da-9759-d76ec5cb7509",
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "# Now, each black hole is itself an object stored in the database. In fact, they are treated similar to halos/galaxies.\n",
    "# Just like a halo/galaxy, black holes store various properties about themselves\n",
    "# You can list off these available properties the same way you would with a galaxy/halo\n",
    "\n",
    "bh = halo['BH_central'][0] # make a new variable \"bh\" to save the first central BH in halo number 48\n",
    "bh.keys()                  # list the available properties associated with this BH\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5c1eec46-dcc6-43c7-bcc9-50aa53818fff",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Some important properties of the black hole objects you'll want to keep in mind:\n",
    "\n",
    "# The ID number associated with the BH. It seems like a random large number and, really, it mostly is. But this is how we keep track of them\n",
    "print(\"BH id: \", bh.halo_number)\n",
    "\n",
    "# The mass\n",
    "print(\"BH Mass: \", bh['BH_mass'])\n",
    "\n",
    "# The instantaneous accretion rate (i.e. how fast it is growing by consuming nearby gas)\n",
    "print(\"BH instantaneous accretion (or growth) rate (in solar masses per year)\", bh['BH_mdot'])\n",
    "\n",
    "# The accretion rate averaged over the last timestep (~1.7 Million Years)\n",
    "print(\"Average BH Accretion rate over the last 1.7 Myr (in solar masses per year): \", bh['BH_mdot_ave'])\n",
    "\n",
    "# The distance from the center of the black hole's host galaxy\n",
    "print(\"Distance to galaxy center (in kpc): \", bh['BH_central_distance'] )\n",
    "\n",
    "# The Host halo of the black hole\n",
    "print(\"BH Host Halo:\", bh['host_halo'])\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ce0bea75-ed43-430f-9b90-363433ad512e",
   "metadata": {},
   "source": [
    "<h3> Black Holes are LINKED to Their Host Galaxies </h3>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "520fe690-e173-4e9b-a47a-aae1c8a53289",
   "metadata": {},
   "source": [
    "Note above that the 'host_halo' property of a black hole returns the host galaxy/halo object within the database, similar to how the \"BH_central\" and \"BH\" properties of a galaxy/halo returns the black hole object. We can often utilize this fact when gathering more information from the simulation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "91e137a1-6945-45e0-b6af-3ad791ba84f2",
   "metadata": {},
   "outputs": [],
   "source": [
    "# The 'host_halo' property returns a halo/galaxy object. \n",
    "\n",
    "bh_host = bh['host_halo']\n",
    "halo44 = db.get_halo('cosmo25/%7779/44')\n",
    "\n",
    "# These two objects are literally the same thing because 'host_halo' and get_halo are returning the exact same object here.\n",
    "\n",
    "print(bh_host, halo44)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "532044ba-b122-4df1-b247-0d9e3c444bc1",
   "metadata": {},
   "outputs": [],
   "source": [
    "# We can even ask python if the two results are the same!\n",
    "\n",
    "# Note the double == here denotes a logical statemt: \"is bh_host exactly equal to/the same as halo44?\"\n",
    "print(bh_host==halo44) #this will print True if both objects are exactly the same\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7650354f-f2fb-41c5-990f-fbf1343a044f",
   "metadata": {},
   "source": [
    "<h3> Gathering Properties of All Black Holes in a Timestep </h3>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "498c9649-50fb-43a5-8d7a-8199bd46eed1",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Because black holes behave just like halos/galaxies in the database, we can use many of the same techniques to look at all the black holes at once\n",
    "\n",
    "#Here we gather data on the ID number, mass, accretion rate, and distance to host center for every black hole in step 7779\n",
    "\n",
    "#first get the step we want\n",
    "step = db.get_timestep('cosmo25/%7779')\n",
    "\n",
    "# Then ask for the properties.\n",
    "# The datebase will search ALL objects for the properties we ask for and give all the results\n",
    "# But only black holes have the ones we are asking for\n",
    "# Similarly, when we asked for Mstar or Mvir, only halos/galaxies had those properties\n",
    "\n",
    "bhid, mass, mdot, dist = step.calculate_all('halo_number()', 'BH_mass', 'BH_mdot', 'BH_central_distance')\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "76fba637-3e1c-4f77-b640-c74e2c90602b",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Now we can plot the proerties we got!\n",
    "# First, a histogram of BH masses\n",
    "# Note, it is useful here to use log10 because masses can range from 10^6 to 10^10 solar masses (over four orders of magnitude!)\n",
    "# Also, typically, we care about the shape of the histogram and less about the total numbers.\n",
    "# Our simulation has a finite volume. To some extend, the total numbers we get are a product of the size. A bigger simulation will have more black holes\n",
    "# But, more fundamental to our simulation is the RELATIVE numbers. So, it is useful to normalize your plot\n",
    "# So, rather than a histogram (Number of Objects vs Some Property) we normalize by both the size of our bins (x-axis ranges) and the total number of objects\n",
    "# This is done using the density=True keyword.\n",
    "\n",
    "plt.hist(np.log10(mass), range=[6,10], bins=20, density=True, histtype='step', lw=3) # lw is another useful keyword. It sets the line width to be thicker when > 1.\n",
    "plt.ylabel('Probability Distribution')                 # this is what we cann the value when normalized in this way\n",
    "plt.xlabel('BH Mass [M$_{\\odot}$]')                    # An example of a useful math command. \\odot creates the symbol we use for \"solar mass\"\n",
    "plt.yscale('log')                                      # Set the y-axis to be log-scale. \n",
    "                                                            # Try commenting this line out with a '#' at the start to see how it looks!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "76a8ed87-7ac4-46aa-8b79-147a4b3e488b",
   "metadata": {},
   "outputs": [],
   "source": [
    "# We can also make scatter plots like before\n",
    "\n",
    "#here c denotes the color and \"k\" stands for black. Some other colors: \"r\" = Red, \"g\" = green, \"b\" = blue, \"o\" = orange, \"c\" = cyan\n",
    "plt.scatter(np.log10(mass), np.log10(mdot), marker='o', c='k')\n",
    "plt.xlabel('BH Mass [M$_{\\odot}$]')\n",
    "plt.ylabel('Accretion Rate [M$_{\\odot}/yr$]')\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ed931594-ccbd-4d47-8b0c-28845758349c",
   "metadata": {},
   "outputs": [],
   "source": [
    "# We can utilize the connection between BHs and their hosts to gather properties of host galaxies for all black holes\n",
    "# For exampe, here we are extracting the BH ID, BH mass, and host galaxy stellar mass\n",
    "# The general format for asking for such a property is \"host_halo.halo_property\".\n",
    "# When calculate_all() sees a property like this, it knows to find the property \"halo_property\" for each \"host_halo\" of each BH.\n",
    "# Note that because a single galaxy/halo can host multiple BHs, doing this will result in multiple BHs having the same halo/galaxy properties\n",
    "\n",
    "bhid, bhmass, host_mstar = step.calculate_all('halo_number()', 'BH_mass', 'host_halo.Mstar')\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "04c27fa1-2436-4d81-9ae7-6af2d035a990",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Now we can plot another histogram, this time using the host galaxy stellar mass\n",
    "# What this tells us is the environment in which most BHs exist\n",
    "# If we were to select a BH at random, in what galaxy are we most likely to find it?\n",
    "\n",
    "plt.hist(np.log10(host_mstar), range=[7,12], bins=20, density=True, histtype='step', lw=3) # lw is another useful keyword. It sets the line width to be thicker when > 1.\n",
    "plt.ylabel('Probability Distribution')                 # this is what we cann the value when normalized in this way\n",
    "plt.xlabel('Host Galaxy Stellar Mass [M$_{\\odot}$]')                    # An example of a useful math command. \\odot creates the symbol we use for \"solar mass\"\n",
    "plt.yscale('log')                                      # Set the y-axis to be log-scale. \n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "37c7d94e-fd5b-4d8b-92d2-ead7bb747635",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Now we can re-create a famous relationship between BH mass and galaxy mass\n",
    "plt.scatter(np.log10(host_mstar), np.log10(bhmass), marker='o', c='k')\n",
    "plt.ylabel('BH Mass [M$_{\\odot}$]')\n",
    "plt.xlabel('Galaxy Stellar Mass [M$_{\\odot}$]')\n",
    "\n",
    "# You might notice that there doesn't seem to be a very good relation at all... \n",
    "# That's because we are plotting using the BHs. But, multiple BHs can exist in the same galaxy\n",
    "# A better way of plotting would be to treat the BH as a property of each *galaxy* \n",
    "# so that each would have one BH mass associated with it\n",
    "# To do this, we will need to use the fact that BHs and galaxies are LINKED"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6d9cab60-4575-47b8-a95b-e2b930ddacd0",
   "metadata": {},
   "source": [
    "<h3> Utilizing Links to Collect BH Properties from Galaxies </h3>\n",
    "\n",
    "Now we can go the other direction and find properties of black holes for each galaxy. Typically, we will want to collect information on either the most massive BH, or the most central, or the one that is growing the fastest. This section shows you how to do this. This will use the \"live calculation\" function link() which is described in some detail <a href=https://github.com/pynbody/tangos/blob/master/docs/live_calculation.md> here</a>.\n",
    "\n",
    "The general syntax to extract a linked property of a halo is `my_halo.calculate('link(my_link_name, link_property, \"max/min\", constraint1, constraint2, ...constrantN)')`\n",
    "\n",
    "In our case for black holes `my_link_name` will (typically) be `BH_central`, `link_property` will be the property you want to select the min/max value of (e.g. `BH_mass`) and constraints1,2,...,N will be logical expressions to further make your selection. For example, stating `BH_central_distance<1` will only select BHs within 1 kpc of the galaxy center."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "83bf45fb-74df-4294-804c-2409adccbc1d",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Ok here is an example.\n",
    "# Here we are finding the central BH with the largest mass that is also within 0.7 kpc of the galaxy's center\n",
    "# The result of this calculation is a BH object itself, which we can extract information from, like it's mass and accretion rate\n",
    "\n",
    "my_halo = db.get_halo('cosmo25/%7779/4')\n",
    "my_bh = my_halo.calculate('link(BH_central,BH_mass, \"max\", BH_central_distance<0.7)')\n",
    "print(my_bh)\n",
    "print(\"My BH's mass: \", my_bh['BH_mass'])\n",
    "print(\"My BH's accretion rate: \", my_bh['BH_mdot'])\n",
    "print(\"My BH's id number: \", my_bh.halo_number)\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cdb2a638-d509-41a5-88b6-dbc7d7c8767b",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Now we can apply these links to collect BH masses for all of our galaxies.\n",
    "# In this case, we want to find the most MASSIVE black hole in each galaxy\n",
    "# We also want to select BHs near the center to avoid including any wandering BHs\n",
    "# we can ask for the property of a linked object by adding `.property_name` after link()\n",
    "# the .BH_mass at the end of the link function will collect BH_mass from each linked BH selected to have the maximum BH mass\n",
    "# and within 0.7 kpc from the galaxy center\n",
    "\n",
    "mstar, mbhcen = sim.timesteps[-4].calculate_all('Mstar','link(BH_central,BH_mass, \"max\", BH_central_distance<0.7).BH_mass')\n",
    "plt.scatter(np.log10(mstar), np.log10(mbhcen), marker='o', c='k')\n",
    "plt.ylabel('BH Mass [M$_{\\odot}$]')\n",
    "plt.xlabel('Galaxy Stellar Mass [M$_{\\odot}$]')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3203c46d-3411-4337-ad06-92a263e2141a",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Now let's look at some different quantities and see how accretion rate relates to star formation rate\n",
    "# For the galaxy, we use the at() function to find the value of the histogram SFR_encl_25Myr at 10 kpc\n",
    "# For the BH, we now collect BH_mdot_ave and also select for the maximum accretion rate instead of mass.\n",
    "\n",
    "sfr_10, mdotbhcen = sim.timesteps[-4].calculate_all('at(10,SFR_encl_25Myr)','link(BH_central,BH_mdot_ave, \"max\", BH_central_distance<0.7).BH_mdot_ave')\n",
    "ok_to_use = (sfr_10 > 0) & (mdotbhcen>0 ) #only select ones with non-zero values for both quantities\n",
    "plt.scatter(sfr_10[ok_to_use], mdotbhcen[ok_to_use], marker='o', c='k') #plot while applying the selection criteria above\n",
    "plt.yscale('log')\n",
    "plt.xscale('log')\n",
    "plt.ylabel('BH Accretion rate [M$_{\\odot}/yr$]')\n",
    "plt.xlabel('Galaxy Star Formation Rate [M$_{\\odot}/yr$]')\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "60a8bd4d-acb9-4bf7-a92c-987a34183a02",
   "metadata": {},
   "outputs": [],
   "source": [
    "# We can make the same plot but colored by different galaxy or BH properties. In this case, I color by BH mass\n",
    "\n",
    "sfr_10, mdotbhcen, mbh = sim.timesteps[-4].calculate_all('at(10,SFR_encl_25Myr)',\n",
    "                                                         'link(BH_central,BH_mdot_ave, \"max\", BH_central_distance<0.7).BH_mdot_ave',\n",
    "                                                         'link(BH_central,BH_mdot_ave, \"max\", BH_central_distance<0.7).BH_mass')\n",
    "ok_to_use = (sfr_10 > 0) & (mdotbhcen>0 ) #only select ones with non-zero values for both quantities\n",
    "plt.scatter(sfr_10[ok_to_use], mdotbhcen[ok_to_use], marker='o', c=np.log10(mbh[ok_to_use]),cmap='viridis') #plot while applying the selection criteria above\n",
    "plt.colorbar(label='log BH Mass [M$_{/odot}$]') #now that we are coloring our points, provide a colorbar!\n",
    "plt.yscale('log')\n",
    "plt.xscale('log')\n",
    "plt.ylabel('BH Accretion rate [M$_{\\odot}/yr$]')\n",
    "plt.xlabel('Galaxy Star Formation Rate [M$_{\\odot}/yr$]')\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "efecdb74-f5cb-4ffe-8ef1-4db2202ad29b",
   "metadata": {},
   "source": [
    "<h3> Things to Try on your Own </h3>\n",
    "\n",
    "1) re-make the Mbh-Mstar and Mdot - SFR plots, coloring by different properties. Some interesting ones to try: \n",
    "    - MColdGas (totdal cold gas mass of the galaxy)\n",
    "    - ColdGasMass_encl (cold gas mass enclosed by different radii. use the at() function for this)\n",
    "    - gas_mass_profile (similar to above but for total gas mass; again use at())\n",
    "    - star_mass_profile (similar to above but for total stellar mass)\n",
    "    - Mstar (total stellar mass of the galaxy)\n",
    "    - BH_mass\n",
    "    - BH_mdot\n",
    "    \n",
    "2) look back at the tangos_intro notebook. based on some of the properties you've made scatter plots for in this notebook and in that one (e.g. mdot-sfr, mstar-mhalo, etc) pick some galaxies to look at in more detail. Look at their images and their time evolution and star formation histories."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c7baa5ca-9279-4b34-b1e0-abb1fc803233",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}