Primordial Black Hole Merger Rate in Ellipsoidal-Collapse Dark Matter Halo Models

Dark matter halos are considered as nonlinear cosmological structures whose masses can be modeled as a function of radius called density profile. In recent years, the analytical models and numerical simulations have provided a clearer picture of the properties and behavior of these structures. In addition, many studies have been performed to extract a convenient density profile that can provide an acceptable description for the observational data dehnen ; Hernquist:1990be ; Jaffe:1983iv ; deZeeuw:1985sk ; Plummer:1911zza ; Navarro:1995iw ; ein .

Let us mention two of the most commonly used density profiles as follows. The first one is the Navarro-Frenk-White (NFW) profile, which is extracted from the $N$-body simulations and is well compatible with most of the rotation curve data Navarro:1995iw . The relation of this density profile is

The other one, which is derived from analytical models, is the Einasto profile that is also well consistent with the observational data ein , and has the form

In these relations, $\rho_{s}$ and $r_{s}$ are the scaled parameters that vary from halo to halo, and $\alpha$ is the shape parameter for the Einasto profile. It should be noted that for both of the above forms, one has

i.e. the logarithmic slope of the density distribution is $-2$.

On the other hand, the halo density profile can be described in terms of two parameters, namely the mass and concentration. The halo concentration describes the central density of halos and is defined as

where $r_{\rm vir}$ is a viral radius considered as a radius within which the average halo density reaches $200$ to $500$ times the critical density of the universe. Also, the $N$-body simulations show that the concentration parameter is a decreasing function of the halo mass and is a redshift dependent function at fixed mass afshordi ; prada ; Dutton:2014xda ; ludlow . We will discuss about the mass distribution of dark matter halos in the next subsection.

The existence of dark matter halos provides a convenient and fundamental framework to study nonlinear gravitational collapse in the universe. Hence, having a proper statistical view of the mass distribution of these halos can improve our understanding of the physics governing those. With this argument, a function called the halo mass function has been proposed, which describes the mass distribution of these halos within a given volume. In other words, the halo mass function describes the mass of those structures whose overdensities exceed the threshold overdensity, separate from the expansion of the universe, and collapse. Furthermore, in the standard cosmology, one may define a linear quantity called density contrast as $\delta(x)\equiv[\rho(x)-\bar{\rho}]/\bar{\rho}$, where $\rho(x)$ is the local density at any point $x$ and $\bar{\rho}$ is the mean background energy density. As noted earlier, this quantity may grow to the critical point while the universe expands, exceeds linear regimes, and enters into nonlinear regimes. This situation occurs when the overdensities separate from the expansion of the universe, enter the turnaround phase and collapse. That is, the structures are formed at this stage. For an Einstein-de Sitter universe and a spherical-collapse halo model, the threshold overdensity has been calculated to be $\delta_{\rm sc}=1.686$. This threshold depends only on the redshift value and is independent of all local quantities such as mass and radius Lukic:2007fc . Hence, it can be considered as a fixed threshold in a narrow redshift range. However, this threshold varies for ellipsoidal-collapse halo models, which we will discuss later.

In Ref. Jenkins:2000bv , in order to specify various fits for dark matter halos, an appropriate definition of the differential mass function has been presented as

where $n(M)$ is the number density of halos with mass $M$, $\rho_{m}$ is the cosmological matter density, and $g(\sigma)$ depends on the geometry of overdensities at the collapse time. The function $\sigma(M,z)$ is the linear root mean square fluctuation of overdensities on mass $M$ and redshift $z$, which is defined as

In this relation, $W(k,M)$ is the Fourier spectrum of the top-hat filter which depends on mass $M$ and wavenumber $k$, and $P(k,z)$ is the power spectrum of the fluctuations.

There is a wide range of studies that have been conducted to extrapolate the halo mass function based on analytical approaches and numerical simulations. The purpose of these studies is to provide the best fit for the cosmic observations. The first model for the dark matter halo mass function, assuming a homogeneous and isotropic collapse, was presented by Press and Schechter (P-S) ps as

which is called the P-S mass function. This formalism is based on the assumption that every astrophysical or cosmological object is formed via a gravitational collapse of overdensities. Moreover, although the final collapse is a nonlinear process, it is assumed that, in the early universe, the density fluctuations had been very small and resulted in a linear approximation. As is clear from relation (7), at a fixed redshift, the mass function depends only on the mass of halos via $\sigma(M)$, and it is expected that no significant change can be observed. This mass function has been proposed as the simplest model for the formation of dark matter halos, i.e. a spherical-collapse model, and, in many cases, is consistent with the observational data. Nevertheless, it quantitatively deviates from the numerical results at some mass limits Jenkins:2000bv . Therefore, some improvements have been made to address this issue. One of the most successful improvements was provided by Sheth and Tormen (S-T), which is based on a more realistic model and fits simulation results better Sheth:1999mn ; st . Their formalism was based on an ellipsoidal-collapse model with dynamical threshold density fluctuations, in contrast to an almost global threshold in the P-S model.

As mentioned earlier, the threshold overdensity for spherical-collapses, $\delta_{\rm sc}$, has been introduced as a global value. It means that, in about certain redshifts, all structures with overdensities higher than such a threshold can collapse. S-T have proposed the idea that dynamically considering the threshold overdensity for the ellipsoidal-collapses, $\delta_{\rm ec}$, can provide a more realistic picture of the halo mass function. With this assumption and considering prolateness to be zero st , they have extracted this quantity as

with $\gamma=0.47$, $\beta=0.615$ and $\nu\equiv\delta_{\rm sc}/\sigma(M)$. It is clear that this quantity not only implicitly depends on the redshift, but also on the mass of the structure, and is called the moving barrier. With this assumption, one can find the halo mass function for the ellipsoidal-collapse, which is also called the S-T mass function, to be

with $F=0.3222$, $a=0.707$ and $p=0.3$. This mass function is expected to be more sensitive than the P-S mass function with redshift changes. Thus, we now have all the tools that one needs to study PBHs in dark matter halos. In this regard, in the next section, we will talk about the probability of encountering PBHs, their binary formation, and their merger rate within a certain volume and time interval.

The PBHs are a special type of black holes that are formed in the early universe due to the direct collapse of density fluctuations or equivalently curvature perturbations. The PBHs were not only able to form binaries when the universe had been dominated by radiation but also could encounter each other in the late time universe due to their random distribution.

More important is that these types of black holes have been considered as candidates for dark matter for over $40$ years. In this regard, it is believed that the gravitational wave events recorded by the LIGO detectors can be described by the PBH scenario at typical mass $30~{}M_{\odot}$, if these black holes could be considered as a component of dark matter. Hence, we have focused on the stellar-mass PBHs that may reside in the galactic halos and have potentially been proposed as a dark matter candidate. Indeed, the recent detections of stellar-mass black hole mergers have drawn much attention to study their origins, and the PBH scenario is one of the plausible theories.

The presence of PBHs with random distributions in dark matter halos gives those a chance to form binaries via the close encounter and emitting gravitational waves. In particular, the numerical simulations Moore:1999nt and analytical approaches Kamionkowski:2008vw to investigate the subhalos resided in the larger halos show that the dark matter in the subhalos has lower velocity dispersion and higher density than the host halos bird . Hence, the subhalos are likely to have the largest contribution to the formation of PBH binaries. Thus, the probability of PBH merger in the subhalos is more significant than in the host halos bird .

Suppose two PBHs with masses $m_{i}$ and $m_{j}$ and relative velocity at the large separation $v_{\rm rel}=|v_{i}-v_{j}|$ accidentally encounter each other in a dark matter halo. At the closest physical separation (i.e., at periastron), due to the maximum scattering amplitude, a significant gravitational radiation occurs. Regarding the emission of the gravitational radiation $E_{\rm rad}$, the Keplerian mechanics implies that the time-averaged radiated energy to be peter:1964

where $G$ is the gravitational constant, $c$ is the velocity of light, $a$ is the semi-major axis of the orbit, $e$ is the eccentricity of the orbit and $r_{\rm p}=a(1-e)$ is the periastron. As the head-on collision rarely happens, it can be assumed that near the periastron, due to the emission of the maximum gravitational radiation, the trajectory is an ellipse with a maximum eccentricity (i.e. $e=1$). Under these assumptions, after one orbital period $T=2\pi\sqrt{a^{3}/[G(m_{i}+m_{j})]}$, the radiated gravitational energy is

If $\Delta E_{\rm rad}>E_{\rm kin}$, where $E_{\rm kin}=\mu v_{\rm rel}^{2}/2$ is the kinetic energy and $\mu=m_{i}m_{j}/(m_{i}+m_{j})$ is the reduced mass, then the PBHs will form binary systems. This condition provides a maximum for the periastron as

which means that $r_{\rm p}<r_{\rm p,max}$ leads to the binary formation. In addition, in the Newtonian limit, the impact parameter $b$ is calculated as a function of $r_{\rm p}$ Sasaki:2018dmp , namely

In other words, if the value of the impact parameter for the encounter is less than $b(r_{\rm p,max})$, one can expect to form a PBH binary. Therefore, the cross-section for the binary formation is obtained to be quinlan ; Mouri:2002mc

where the strong limit of the gravitational focusing, i.e. $r_{\rm p}\ll b$, has been applied.

Our focus is on the merger rate of the PBHs that are compatible with the LIGO sensitivity, i.e. $\sim\left(30~{}M_{\odot}-30~{}M_{\odot}\right)$ events in the galactic halos. Accordingly, we have normalized those masses to $30~{}M_{\odot}$ with their relative velocities as the average velocities of dark matter halos, i.e. $200~{}\rm km/s$.

By inserting Eq. (12) into Eq. (14) and assuming $m_{i}=m_{j}=M_{\rm pbh}$ and $v_{\rm rel}=v_{\rm pbh}$, one can reach an explicit form of the cross-section related to the normalized mass and velocity of the PBHs as

With these considerations, the merger rate of PBHs within each halo can be calculated via the formula bird ; Nishikawa:2017chy

where $\rho_{\rm halo}(r)$ is the halo density profile that can be considered to be the NFW or the Einasto density profile, and $\langle\xi v_{\rm pbh}\rangle$ represents an average over the relative velocity distribution of PBHs in the galactic halos.

Moreover, the mass located within the virial radius of the halo, the virialized mass, can be found by

By inserting relation (1) into relation (17) and integrating, one can find the virialized mass for the NFW density profile as

Similarly, by considering relation (2), the virialized mass for the Einasto density profile dan ; retana can be obtained as

In this relation, $l(C,\alpha)$ is a function of concentration and shape parameter and has the following form

where $\Gamma(x,y)=\int_{0}^{y}t^{x-1}e^{-t}dt$ is the incomplete Gamma function.

To calculate the halo velocity dispersion, one can use its relation to the maximum velocity in an $r_{\rm max}$ radius, which has been introduced in Ref. prada as

In this work, we assume that the relative velocity distributions of PBHs in a halo are random and follow the Maxwell-Boltzmann statistics. Hence, one can write the velocity probability distribution function as

where $A_{0}$ is determined when the condition $4\pi\int_{0}^{v_{\rm vir}}P(v)v^{2}dv=1$ is satisfied and a cutoff is considered at the halo virial velocity.

It is clear from Eqs. (16), (18) and (19) that, in order to calculate the merger rate in each halo, the mass-concentration relation, $C(M_{\rm vir})$, has to be determined. For this purpose, according to the initial conditions governing the dark matter halos during the collapse, various results can be found. In Ref. bird , the merger rate has been performed by using two famous spherical concentration-mass relations of Prada, et al. prada and Ludlow, et al. ludlow .

In this research, we have employed the ellipsoidal-collapse concentration-mass relation introduced in Ref. afshordi that we refer to it as Okoli-Afshordi (O-A) concentration-mass relation. Furthermore, for the Einasto density profile, we have chosen the value of the shape parameter presented in Ref. Klypin:2014kpa . Also, we have set the mass of PBHs to be $30~{}M_{\odot}$. In Fig. 1, we have indicated the merger rate of PBHs per halo as a function of halo mass by considering the O-A concentration-mass relation as an ellipsoidal model, and the Ludlow concentration-mass relation for a spherical model obtained in Ref. bird . The results show that the merger rate of PBHs grows with increasing halo mass for both models. A noteworthy point is that in smaller mass halos for the case of the O-A model the merger rate is almost one order of magnitude larger than the corresponding result for the Ludlow model.

In the following, we propose to determine the effect of these changes on the total merger rate of PBHs in a given volume and time interval.