docs/GW_events.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "817490f7-38eb-4c38-863b-6b9492f455d8",
   "metadata": {},
   "source": [
    "# Gravitational wave event rates"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "459bd779",
   "metadata": {},
   "source": [
    "This document outlines the default settings for calculating detectable gravitational wave event rates (per year) in [LeR](https://ler.readthedocs.io/en/latest/index.html).\n",
    "\n",
    "**Parameters Considered:**\n",
    "\n",
    "* $z_s$: Source redshift.\n",
    "* $\\theta$: Gravitational wave source properties. $m_1$ (mass of the heavier one), $m_2$ (mass of the lighter one), $D_l$ (luminosity distance), $\\iota$ (inclination-angle), $\\phi$ (phase-of-coalescence), $\\psi$ (polarization-angle), $ra$ (right-ascension), $dec$ (declination), $t$ (time-of-coalescence).\n",
    "* $\\rho(z_s,\\theta)$: Network optimal SNR.\n",
    "* $\\rho_{th}$: SNR threshold.\n",
    "* $R_o^U(z_s)$: Normalized merger-rate distribution (source-frame).\n",
    "* $\\mathcal{N}^U$: Normalizing factor.\n",
    "* $P$: Prior distribution.\n",
    "* $\\Theta$: Detectability condition.\n",
    "\n",
    "Default cosmology: LambdaCDM(H0=70, Om0=0.3, Ode0=0.7). But, it can be changed by the user at LeR initialization."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1e0a7f27",
   "metadata": {},
   "source": [
    "## Compact Binaries Merger Rate Density\n",
    "\n",
    "The merger rate density is a crucial parameter for understanding the population of compact binary mergers, including binary black holes (BBH), binary neutron stars (BNS), and neutron star-black hole (NSBH) systems. It quantifies the number of mergers occurring per unit volume per unit time. The merger rate density is typically expressed in units of $\\text{Gpc}^{-3} \\text{yr}^{-1}$ and varies with redshift ($z$).\n",
    "\n",
    "For available models for merger rate density, see the example section.\n",
    "\n",
    "### Merger Rate Density Formula\n",
    "\n",
    "The default merger rate density in LeR follows the model presented in [WIERDA et al. 2021](https://arxiv.org/pdf/2106.06303.pdf). This model is a functional fit to the population I/II star merger-rate density, normalized to the local merger rate density. It extends from the M10 model to the Belczynski et al. (2017) model, which incorporates the metallicity dependence of the star formation rate, influenced by pair-instability supernova (PSN) and pair-instability pulsation supernova (PPSN) effects.\n",
    "\n",
    "The merger rate density in the detector frame is given by:\n",
    "\n",
    "\\begin{equation}\n",
    "\\mathcal{R}_m(z_s) = \\frac{\\mathcal{R}_O(b_4+1)e^{b_2 z_s}}{b_4+e^{b_3 z_s}} \\text{Gpc}^{-3}\\text{yr}^{-1} \\tag{1}\n",
    "\\end{equation}\n",
    "\n",
    "- Where $z_s$ is the redshift of the source.\n",
    "- The local merger rate $\\mathcal{R}_O$ and fitting parameters $b_2$, $b_3$, and $b_4$ are specific to each model.\n",
    "\n",
    "The following table provides the local merger rate $\\mathcal{R}_O$ for different compact binary models (refer to [GWTC-3 IV A](https://journals.aps.org/prx/pdf/10.1103/PhysRevX.13.011048)):\n",
    "\n",
    "| Model | $\\mathcal{R}_O$ ($\\text{Gpc}^{-3} \\text{yr}^{-1}$) |\n",
    "|-------|------------------------------------|\n",
    "| BNS   | $105.5^{+190.2}_{-83.9}$           |\n",
    "| BBH   | $23.9^{+14.9}_{-8.6}$              |\n",
    "| NSBH  | $45^{+75}_{-33}$                   |\n",
    "\n",
    "### Merger Rate Density at Source Redshift\n",
    "\n",
    "The merger rate density at source redshift $z_s$ in the source frame, denoted by $R(z_s)$, is given by:\n",
    "\n",
    "\\begin{equation} \n",
    "R(z_s) = \\frac{\\mathcal{R}_m(z_s)}{1+z_s} \\frac{dV_c}{dz_s} \\tag{3}\n",
    "\\end{equation}\n",
    "\n",
    "- The factor $(1+z_s)$ in the denominator accounts for time dilation.\n",
    "- $\\frac{dV_c}{dz_s}$ represents the differential comoving volume.\n",
    "- This formula converts the merger rate density from the detector frame to the source frame.\n",
    "\n",
    "### Normalization of Merger Rate Distribution\n",
    "\n",
    "The normalizing factor of the merger rate distribution with respect to $z_s$, denoted by $\\mathcal{N}^U$, is given by:\n",
    "\n",
    "\\begin{equation}\n",
    "\\mathcal{N}^U = \\int_{z_{\\text{min}}}^{z_{\\text{max}}} dz_s \\frac{R_m(z_s)}{1+z_s} \\frac{dV_c}{dz_s} \\tag{4}\n",
    "\\end{equation}\n",
    "\n",
    "The normalized merger-rate distribution, $R_o^U(z_s)$, is expressed as:\n",
    "\n",
    "\\begin{equation} \n",
    "R_o^U(z_s) = \\frac{1}{\\mathcal{N}^U} \\frac{R_m(z_s)}{1+z_s} \\frac{dV_c}{dz_s} \\tag{5}\n",
    "\\end{equation}\n",
    "\n",
    "Distribution of $R_o^U(z_s)$ wrt $z_s$ is shown below for BBH.\n",
    "\n",
    "<p align=\"center\">\n",
    "  <img src=\"_static/merger_rate_bbh.png\" alt=\"zs\">\n",
    "</p>\n",
    "\n",
    "\n",
    "### Sampling Method\n",
    "\n",
    "Sampling of the source redshift $z_s$ is based on the inverse transform sampling method. An interpolation of the inverse cumulative distribution function (CDF) is created and stored for future use, making the sampling process faster and more efficient.\n",
    "\n",
    "**Note:** The merger rate density models and parameters are subject to updates based on new observational data and theoretical developments."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "af41ebf5",
   "metadata": {},
   "source": [
    "## Annual Rate of Detectable Gravitational Wave Events\n",
    "\n",
    "Continuing from the discussion on Compact Binaries Merger Rate Density, we now focus on the annual rate of detectable gravitational wave events. This rate, denoted as $R_U$, quantifies the expected number of gravitational wave events that can be observed by detectors within a year. The expression for $R_U$ is given by:\n",
    "\n",
    "\\begin{equation} \n",
    "R_U = \\mathcal{N}^U\\int dz_s R_o^U(z_s)\\left\\{\\Theta[\\rho(z_s,\\theta)-\\rho_{th}] P(\\theta) d\\theta \\right\\} \\tag{2}\n",
    "\\end{equation}\n",
    "\n",
    "Key aspects of this equation include:\n",
    "\n",
    "- **SNR Operator,** $\\mathbf{\\Theta[\\rho(z_s,\\theta)-\\rho_{th}]}$**:** This is a step function that determines whether the signal-to-noise ratio (SNR) $\\rho(z_s,\\theta)$ of an event exceeds the threshold SNR $\\rho_{th}$. If $\\rho(z_s,\\theta) \\geq \\rho_{th}$, the event is considered detectable.\n",
    "- The SNR threshold $\\mathbf{\\rho_{th}}$ is typically chosen based on the sensitivity and design specifications of the gravitational wave detector. This strategy is valid as a proxy in the Gaussian noise regime.\n",
    "- **Prior Distribution,** $\\mathbf{P(\\theta)}$**:** This normalized distribution represents the prior knowledge or assumptions about the distribution of gravitational wave parameters $\\theta$.\n",
    "- **Normalized Merger-Rate Distribution,** $\\mathbf{R_o^U(z_s)}$**:** This function represents the merger rate density (in the source frame) normalized with respect to the redshift $z_s$.\n",
    "- **Gravitational Wave Parameters,** $\\mathbf{\\theta}$**:** The set of parameters includes the masses of the binary components ($m_1$ and $m_2$), the inclination angle ($\\iota$), the polarization angle ($\\psi$), the right ascension ($ra$), the declination ($dec$), the phase of coalescence ($\\phi_c$), and the time of coalescence ($t_c$).\n",
    "- **Luminosity Distance,** $\\mathbf{\\mathcal{D}_L}$**:** The luminosity distance is derived from the source redshift $z_s$ and cosmological parameters, providing a measure of the distance to the source of the gravitational waves.\n",
    "\n",
    "The integration over $z_s$ and $\\theta$ in the equation accounts for the contributions from all detectable events across different redshifts and with varying gravitational wave parameters. This comprehensive approach enables the estimation of the annual rate of detectable gravitational wave events, which is crucial for understanding the observational capabilities of gravitational wave detectors and planning future observational strategies.\n",
    "\n",
    "Note: Refer to the section [GWSNR](https://gwsnr.readthedocs.io/en/latest/) for details on the calculation of the optimal signal-to-noise ratio (SNR)."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2324d5b3",
   "metadata": {},
   "source": [
    "## Compact binaries mass distribution"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "aec197c9",
   "metadata": {},
   "source": [
    "### BBH Mass Distribution: PowerLaw+Peak Model\n",
    "\n",
    "The PowerLaw+Peak model is a parametric model used to describe the mass distribution of binary black hole (BBH) mergers. It combines a power-law component with a Gaussian peak to capture both the broad distribution and potential local enhancements in the mass spectrum of BBHs. The model is defined by the following parameters:\n",
    "\n",
    "- $m_{\\text{minbh}}$: Minimum black hole mass in the distribution.\n",
    "- $m_{\\text{maxbh}}$: Maximum black hole mass in the distribution.\n",
    "- $\\alpha$: Slope of the power-law component.\n",
    "- $\\mu_g$: Mean of the Gaussian peak.\n",
    "- $\\sigma_g$: Standard deviation of the Gaussian peak.\n",
    "- $\\lambda_{\\text{peak}}$: Fractional contribution of the Gaussian peak to the overall distribution.\n",
    "- $\\delta_m$: Smoothing parameter for the transition between the power-law and Gaussian components.\n",
    "- $\\beta$: Exponent for the mass ratio distribution.\n",
    "\n",
    "For analytical expressions and details on the PowerLaw+Peak model, refer to the [GWTC-2 2021](https://arxiv.org/pdf/2010.14533.pdf) (Appendix B3). The default values for the parameters are based on [GWTC-3 2021](https://arxiv.org/pdf/2111.03604.pdf) (page-23). [LeR](https://ler.readthedocs.io/en/latest/index.html) uses this model through [gwcosmo](https://lscsoft.docs.ligo.org/gwcosmo/).\n",
    "\n",
    "Here is the default parameters and resultant distribution:\n",
    "\n",
    "<table>\n",
    "<tr>\n",
    "<td>\n",
    "\n",
    "| Parameter            | Value  |\n",
    "|----------------------|--------|\n",
    "| $m_{\\text{minbh}}$   | 4.98   |\n",
    "| $m_{\\text{maxbh}}$   | 112.5  |\n",
    "| $\\alpha$             | 3.78   |\n",
    "| $\\mu_g$              | 32.27  |\n",
    "| $\\sigma_g$           | 3.88   |\n",
    "| $\\lambda_{\\text{peak}}$ | 0.03 |\n",
    "| $\\delta_m$           | 4.8    |\n",
    "| $\\beta$              | 0.81   |\n",
    "\n",
    "</td>\n",
    "<td>\n",
    "\n",
    "![mass_bbh](_static/bbh_mass.png)\n",
    "\n",
    "</td>\n",
    "</tr>\n",
    "</table>\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bfd94ccf",
   "metadata": {},
   "source": [
    "### BNS mass distribution: Bimodal Gaussian Model\n",
    "\n",
    "The mass distribution of neutron stars in binary neutron star (BNS) systems can be characterized by a bimodal distribution, which can be modeled using a double Gaussian distribution. This bimodal distribution is thought to reflect the different evolutionary paths that neutron stars can take. For more details, refer to [Alsing et al. 2018](https://arxiv.org/pdf/1805.06442.pdf) and [M. Farr et al. 2020](https://arxiv.org/pdf/2005.00032.pdf).\n",
    "\n",
    "The probability density function for the neutron star mass, $m$, is given by:\n",
    "\n",
    "\\begin{equation}\n",
    "p(m) = wN (\\mu_L, \\sigma_L) + (1 − w)N (\\mu_R, \\sigma_R)\n",
    "\\end{equation}\n",
    "\n",
    "where:\n",
    "- $w$: Fraction of neutron stars in the low-mass peak.\n",
    "- $\\mu_L$ and $\\mu_R$: Mean masses of the low-mass and high-mass peaks, respectively.\n",
    "- $\\sigma_L$ and $\\sigma_R$: Standard deviations of the low-mass and high-mass peaks, respectively.\n",
    "\n",
    "Neutron stars that form from the collapse of massive stars in binary systems can either remain isolated or undergo mass transfer from their companion stars. Neutron stars that remain isolated typically have a lower mass distribution, with a peak at around 1.35 $M_{\\odot}$. On the other hand, neutron stars that undergo mass transfer from their companion stars can have a higher mass distribution, with a peak at around 1.8 $M_{\\odot}$.\n",
    "\n",
    "Here is the default parameters and resultant distribution:\n",
    "\n",
    "<table>\n",
    "<tr>\n",
    "<td>\n",
    "\n",
    "| Parameter   | Value            |\n",
    "|-------------|------------------|\n",
    "| $w$         | 0.643            |\n",
    "| $\\mu_L$     | 1.352 $M_{\\odot}$|\n",
    "| $\\sigma_L$  | 0.08 $M_{\\odot}$ |\n",
    "| $\\mu_R$     | 1.88 $M_{\\odot}$ |\n",
    "| $\\sigma_R$  | 0.3 $M_{\\odot}$  |\n",
    "| $m_{\\text{max}}$ | 2.3 $M_{\\odot}$ |\n",
    "\n",
    "</td>\n",
    "<td>\n",
    "\n",
    "![mass_bns](_static/bns_mass.png)\n",
    "\n",
    "</td>\n",
    "</tr>\n",
    "</table>\n",
    "\n",
    "Each normal distribution is independently truncated and normalized in the range [1, 2.3] $M_{\\odot}$, ensuring that the neutron star masses are within the physically plausible range.\n",
    "\n",
    "This model has been found to be a good fit for the observed mass distribution of neutron stars in binary systems.\n",
    "\n",
    "Note: Both primary and secondary masses are sampled from the same distribution. However, if the secondary mass is greater than the primary mass, their values are exchanged to maintain the convention that the primary mass is the larger of the two."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "90a54ee5",
   "metadata": {},
   "source": [
    "## Sampling Methods for Compact Binary Sources\n",
    "\n",
    "The following table summarizes the sampling methods and ranges for various parameters of compact binary sources:\n",
    "\n",
    "| Parameter | Unit                 | Sampling Method        | Range               |\n",
    "|-----------|----------------------|------------------------|---------------------|\n",
    "| $z_s$     | None                 | Merger rate $R_o^U(z_s)$ | [0, 10]            |\n",
    "| $m_1, m_2$ | $\\mathcal{M}_{\\odot}$ | PowerLaw+Peak model   | [$m_{\\text{min}}$, $m_{\\text{max}}$] |\n",
    "| $ra$      | Radian               | Uniform                | [0, $2\\pi$]         |\n",
    "| $dec$     | Radian               | Cosine                 | [$-\\pi/2$, $\\pi/2$] |\n",
    "| $\\iota$   | Radian               | Sine                   | [0, $\\pi$]          |\n",
    "| $\\psi$    | Radian               | Uniform                | [0, $\\pi$]          |\n",
    "| $\\phi_c$  | Radian               | Uniform                | [0, $2\\pi$]         |\n",
    "| $t_c$     | Second               | Uniform (one year)     | [$t_{\\text{min}}$, $t_{\\text{max}}$] |"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "47581a9c",
   "metadata": {},
   "source": [
    "## Statistical Form of Detectable Merger Rate (Monte Carlo Integration)\n",
    "\n",
    "The detectable merger rate, $R_U$, can be statistically estimated using Monte Carlo integration. This can be done using either a step function or the probability of detection ($P_{\\text{det}}$).\n",
    "\n",
    "- **With Step Function:**\n",
    "\n",
    "\\begin{equation} \n",
    "\\begin{split}\n",
    "R_U &= \\mathcal{N}^U \\left< \\Theta[\\rho(z_s,\\theta)-\\rho_{th}]\\right>_{z_s,\\theta \\in \\text{Sampled}[z_s,\\theta]} \\\\\n",
    "&= \\frac{\\mathcal{N}^U}{N} \\sum_i \\Theta[\\rho(z_s^i,\\theta_i)-\\rho_{th}]\n",
    "\\end{split} \\tag{6}\n",
    "\\end{equation}\n",
    "\n",
    "- **With $P_{\\text{det}}$:**\n",
    "\n",
    "\\begin{equation} \n",
    "\\begin{split}\n",
    "R_U &= \\mathcal{N}^U \\left<  P_{\\text{det}}(z_s,\\theta)\\right>_{z_s,\\theta \\in \\text{Sampled}[z_s,\\theta]} \\\\\n",
    "&= \\frac{\\mathcal{N}^U}{N} \\sum_i P_{\\text{det}}(z_s^i,\\theta_i)\n",
    "\\end{split} \\tag{7}\n",
    "\\end{equation}\n",
    "\n",
    "In both cases, the average is taken over a sample of $(z_s, \\theta)$, where $N$ is the total number of samples. The step function $\\Theta[\\rho(z_s,\\theta)-\\rho_{th}]$ evaluates to 1 if the SNR $\\rho(z_s,\\theta)$ is above the threshold $\\rho_{th}$, and 0 otherwise. The probability of detection, $P_{\\text{det}}(z_s,\\theta)$, provides a smoother transition from detectable to non-detectable events at the SNR threshold.\n",
    "\n",
    "The table below presents the merger rate values for two different detector configurations, considering spin-less systems for simplicity. The signal-to-noise ratio (SNR) operator is implemented using a step function with a threshold of 8, and the waveform approximant used is IMRPhenomD.\n",
    "\n",
    "| Detector Configuration | BBH (popI-II) Merger Rate (yr<sup>-1</sup>) | BNS Merger Rate (yr<sup>-1</sup>) | Ratio BBH:BNS |\n",
    "|-------------------|---------|---------|-------|\n",
    "| [L1, H1, V1] (O4) | 440.9   | 3.2     | 137.8 |\n",
    "| [CE, ET] (3G)     | 91942.8 | 38810.4 | 2.4   |     "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "aa94cd90",
   "metadata": {},
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.10.14"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}