Gravitational wave events sampling (Analytical formulation)

Initial setup

Note: I will interchangeably use terms like observable events and detectable events.

Define all the parameters involved.

Source parameters: \(\theta \in \{\) \(m_1\) (mass of the heavier one), \(m_2\) (mass of the lighter one), \(\iota\) (inclination-angle), \(\phi\) (phase-of-coalescence), \(\psi\) (polarization-angle), \(ra\) (right-ascension), \(dec\) (declination) \(\}\) and \(z_s\) : red-shift of the source.

Given \(d N^U_{obs}(z_s)\) is the number of detectable gravitational wave (GW) events from sources at red-shift \(z_s\) in a spherical shell of thickness \(d z_s\), then, the rate of observing GWs is given by,

\[\begin{split}\begin{equation} \begin{split} \mathcal{R}_U &= \int_{z_{min}}^{z_{max}} \frac{d N^U_{obs}(z_s)}{d t} \\ &= \int_{z_{min}}^{z_{max}} \frac{d^2 N^U_{obs}(z_s)}{d t \;d\mathcal{V}_c} \frac{d\mathcal{V}_c}{d\Omega dz_s} d\Omega dz_s \\ &= \int_{z_{min}}^{z_{max}} \frac{d^2 N^U_{obs}(z_s)}{\left(d t/(1+z_s)\right) \;d\mathcal{V}_c}\frac{1}{(1+z_s)}\frac{d\mathcal{V}_c}{d\Omega dz_s} d\Omega dz_s \\ &= \int_{z_{min}}^{z_{max}} \frac{d^2 N^U_{obs}(z_s)}{d\tau \;d\mathcal{V}_c}\frac{1}{(1+z_s)} \frac{d\mathcal{V}_c}{d\Omega dz_s} d\Omega dz_s \\ &= \int_{z_{min}}^{z_{max}} \frac{d^2 N^U_{obs}(z_s)}{d\tau \;d\mathcal{V}_c}\frac{1}{(1+z_s)} \frac{d\mathcal{V}_c}{d\Omega dz_s} \Bigg\{ \int_{all\; sky} d\Omega \Bigg\}dz_s \\ &= \int_{z_{min}}^{z_{max}} \frac{d^2 N^U_{obs}(z_s)}{d\tau \;d\mathcal{V}_c}\frac{1}{(1+z_s)} \Bigg(\frac{d\mathcal{V}_c}{d\Omega dz_s} 4\pi \Bigg)dz_s \\ &= \int_{z_{min}}^{z_{max}} \frac{d^2 N^U_{obs}(z_s)}{d\tau \;d\mathcal{V}_c}\frac{1}{(1+z_s)} \frac{dV_c}{dz_s} dz_s \\ \end{split} \tag{1} \end{equation}\end{split}\]

Note that if \(d\tau=dt/(1+z_s)\) is considered the proper time and it can be converted from the time at detector frame \(dt\) using the time-dilation factor \((1+z_s)\). Consequently, \(\frac{d^2 N^U_{obs}(z_s)}{d t \;d\mathcal{V}_c}\) and \(\frac{d^2 N^U_{obs}(z_s)}{d \tau \;d\mathcal{V}_c}\) are the observed merger rate density at detector-frame and source-frame respectively. We want to use the \(R^U_{obs}(z_s)=\frac{d^2 N^U_{obs}(z_s)}{d \tau \;d\mathcal{V}_c}\) for our analysis as most observational papers and the output of theoretical predictions are in the source-frame. \(\frac{dV_c}{dz_s}dz_s\) is considered a spherical-shell volume element in co-moving coordinates at red-shift \(z_s\). So, the rate expression simplifies to integrating (density) \(\times\) (time dilation effect) over the shell volume element. For more information on the volume element refer to this `page <>`__ of the documentation.

Note: \(\frac{dV_c}{dz_s}\) is the differential co-moving volume at red-shift \(z_s\) and you can get the value by using astropy cosmology package for a given cosmology.

Now we have,

\[\begin{split}\begin{equation} \begin{split} \mathcal{R}_U &= \int_{z_{min}}^{z_{max}} \frac{d^2 N^U_{obs}(z_s)}{d \tau\;d\mathcal{V}_c} \frac{1}{(1+z_s)} \frac{dV_c}{dz_s} dz_s \\ &= \int_{z_{min}}^{z_{max}} R^U_{obs}(z_s) \frac{1}{(1+z_s)} \frac{dV_c}{dz_s} dz_s \end{split} \tag{2} \end{equation}\end{split}\]

Equation 2 is exactly same to the of equation A2 from WIERDA et al. 2021.

We want to re-write it in terms intrinsic merger rate distribution \(R(z_s)\). \(R(z_s)\) is the merger rate density distribution regardless of whether the source is detectable or not.

\[\begin{split}\begin{equation} \begin{split} \mathcal{R}_U &= \int_{z_{min}}^{z_{max}} R^U_{obs}(z_s) \frac{1}{(1+z_s)} \frac{dV_c}{dz_s} dz_s \\ &= \int_{z_{min}}^{z_{max}} R(z_s) P(obs|z_s) \frac{1}{(1+z_s)} \frac{dV_c}{dz_s} dz_s \\ \end{split} \tag{3} \end{equation}\end{split}\]

\[\begin{split}\begin{equation} \begin{split} \mathcal{R}_U &= \mathcal{N}^U \int_{z_{min}}^{z_{max}} P(z_s) P(obs|z_s) dz_s \\ \end{split} \tag{4} \end{equation}\end{split}\]

\(P(obs|z_s)\) is the probability of observing a GW event given that the source is at red-shift \(z_s\). Here, it is considered marginalized over all the gravitational wave parameters except the red-shift, i.e., \(P(obs|z_s) = \int_{\theta} d\theta P(obs|z_s, \theta) P(\theta)\).

In Eqn(4) the merger rate density becomes the prior distribution \(P(z_s)\) for the red-shift, and it is normalized over the red-shift range. So the normalizing factor is given by,

\[\begin{equation} \begin{split} \mathcal{N}^U = \int_{z_{min}}^{z_{max}} R(z_s) \frac{1}{(1+z_s)} \frac{dV_c}{dz_s} dz_s \end{split} \tag{5} \end{equation}\]

Therefore,

\[\begin{equation} \begin{split} P(z_s) = \frac{R(z_s)}{\mathcal{N}^U} \frac{1}{(1+z_s)} \frac{dV_c}{dz_s} \end{split} \tag{6} \end{equation}\]

This is done because we want to sample the red-shift values from this prior distribution, which later will be used in the monte-carlo integration to estimate the rate.

Rate equation (statistical formulation)

Let’s now consider all the gravitational waves source parameters \(\theta \in \{m_1,m_2,z_s,\iota,\phi,\psi,ra,dec,\text{spins}\}\),

\begin{equation} \begin{split} \mathcal{R}_U &= \mathcal{N}^U \int_{z_{min}}^{z_{max}} P(z_s) P(obs|z_s) dz_s \\ & \text{consider} \int \rightarrow \int_{z_s}\int_{\theta} \rightarrow \int_{z_s}\int_{m_1}\int_{m_2}\int_{\iota}\int_{\phi}\int_{\psi}\int_{ra}\int_{dec}\int_{\text{spins}} \\ \mathcal{R}_U &= \mathcal{N}^U \int P(z_s) P(obs|z_s, \theta) P(\theta) d\theta dz_s \end{split} \tag{8} \end{equation}

Here, I have assumed \(P(\theta)\) is a normalized multidimensional pdf. All parameters in \(\theta\), other than \(m_1\) and \(m_2\), are sampled independently. The pdf \(P(\theta)\) is a product of individual pdfs for each parameter.

Final expression for the rate of observing GW events is given by,

\begin{equation} \begin{split} \mathcal{R}_U &= \mathcal{N}^U \bigg< P(obs|\theta, z_s) \bigg>_{z_s\in P(z_s), \theta\in P(\theta)} \\ P(\theta) &= P(m_1,m_2)P(\iota)P(\phi)P(\psi)P(ra)P(dec)P(\text{spins}) \end{split} \tag{4} \end{equation}

The order of sampling and rate calculation steps in LeR are listed below.

Sample \(z_s\) from \(P(z_s)\).
Sample \(\theta\) from \(P(\theta)\).
Calculate SNR with gwsnr
Apply the SNR threshold and check whether the event is detectable or not.
Calculate rate of GW events, \(\mathcal{R}_U\) using equation 4.

[ ]: