Primordial Black Hole Merger Rate in f(R) Gravity

Primordial black holes (PBHs) are known as one of the potential candidates for dark matter. They are expected to have formed due to the direct gravitational collapse of density fluctuations in the early Universe. In this regard, examining the merger rate of PBHs within modified theories of gravity can offer a deeper insight into their abundance. In this work, we delve into the calculation of the merger rate of PBHs within the theoretical framework of f(R) gravity. Our analysis reveals an enhancement in the merger rate of PBHs compared to that obtained from general relativity. Additionally, modulating the field strength f R0 induces shifts in the PBH merger rate, presenting a potential observational signature of modified gravity. We also explore the upper bounds on the abundance of PBHs obtained from f(R) gravity models by comparing the results with gravitational-wave and observational data. The results indicate that in certain regions not excluded by benchmarking data, the parameter space for these upper bounds may be considered reliable.


INTRODUCTION
Gravitational waves (GWs), as cosmological phenomena, have been a central area of investigation for several decades.Consequently, the detection of GWs has introduced a novel framework for analyzing astrophysical and cosmological phenomena.One of the main sources of GWs is the merging of binary black holes, which produces distinguishable signals in GW detectors (Mandel & Broekgaarden 2022;Olsen et al. 2022).In recent years, numerous binary black hole mergers have been identified through the collaborative efforts of the Laser Interferometer Gravitational-Wave Observatory (LIGO)-Virgo-KAGRA detectors (Abbott et al. 2016a(Abbott et al. ,b,c, 2020a,b),b).This accomplishment has opened a fresh avenue for comprehending compact objects distributed throughout the Universe.
Notably, the binary black hole mergers identified by the LIGO-Virgo-KAGRA observatories are related to those of stellar mass, i.e., those with a mass less than 100 M ⊙ .This revelation holds the potential to offer significant insights into the distribution of black hole s fakhry@sbu.ac.ir masses across the Universe, as well as their mechanisms of formation.The exact processes responsible for the creation of black holes involved in binary merger events detected by the LIGO-Virgo-KAGRA detectors are yet to be conclusively determined.These black holes could originate from the collapse of massive stars (potentially through diverse channels) (Fishbach et al. 2021;Rodriguez et al. 2021), or they might have arisen from the direct collapse of primordial density fluctuations during the early stages of the Universe (Bird et al. 2016;Sasaki et al. 2016;Clesse & García-Bellido 2017).
It is captivating to note that the GW information gathered by the LIGO-Virgo-KAGRA detectors corresponds closely to the idea of mergers involving PBHs, which are among the potential explanations put forth to account for dark matter, see, e.g., (Spergel & Steinhardt 2000;Alcock et al. 2000;Liebling & Palenzuela 2012;Roszkowski et al. 2018;Braine et al. 2020;Di Giovanni et al. 2020, 2021;Stasenko & Belotsky 2023).The formation of PBHs is theorized to occur due to nonlinear peaks within initial density fluctuations.In fact, when these fluctuations surpass a certain threshold value, they may directly collapse into PBHs (Zel'dovich & Novikov 1967;Hawking 1971;Volonteri & Bellovary 2012).Various studies have explored the development of nonlinear density fluctuations on scales larger than the horizon, as well as the critical amplitude of curvature perturbations required for black hole formation (Niemeyer & Jedamzik 1999;Shibata & Sasaki 1999;Polnarev & Musco 2007;Young et al. 2014;Allahyari et al. 2017).Additionally, PBHs can cover a broad spectrum of masses, setting them apart from astrophysical black holes (Sasaki et al. 2018).
Due to their stochastic spatial distribution, PBHs within the confines of dark matter halos could potentially encounter each other or other compact objects, forming binary systems that emit GWs and ultimately undergo merger.In recent years, the detection of GWs stemming from the coalescence of binary black holes by the LIGO-Virgo-KAGRA detectors has sparked renewed fascination with PBHs (Abbott et al. 2019(Abbott et al. , 2021;;The LIGO Scientific Collaboration et al. 2021).Given the mass range of black holes engaged in binary mergers, which often surpasses the spectrum of masses observed in typical astrophysical black holes, it becomes plausible that these black holes originated in early Universe.Considering PBHs as conceivable candidates for dark matter and acknowledging the possibility of their clustering within dark matter halos, the structural properties of these halos could potentially influence the merger rate of PBHs.Thus, it has been suggested that certain quantities, such as the halo mass function, halo concentration parameter, and halo density profile, could enhance the precision of theoretical models and their prognostications concerning the merger rate of PBHs.For more details, see our previous works (Fakhry et al. 2021(Fakhry et al. , 2022a,b;,b;Fakhry & Del Popolo 2023;Fakhry et al. 2023a,b).
Several models have been suggested to offer an elaborate depiction of the halo mass function (Press & Schechter 1974;Sheth & Tormen 1999;Sheth et al. 2001;Jenkins et al. 2001;Reed et al. 2003;Warren et al. 2006;Reed et al. 2007;Del Popolo et al. 2017), and halo concentration parameter (Prada et al. 2012;Okoli & Afshordi 2016;Ludlow et al. 2016;Okoli 2017) within the context of GR.However, a fundamental inquiry arises regarding whether the principles of GR and the ideas derived from it provide an adequate framework for explaining the formation and evolution of large-scale structures like galactic halos or not.Nowadays, there is a growing consensus that, while GR has achieved impressive feats in accurately predicting numerous events, it falls short when it comes to explaining large-scale phenomena such as dark sectors of the Universe, see, e.g., (Sherwin et al. 2011;Yang & Xu 2014;Huterer & Shafer 2018;Di Valentino et al. 2020;Abdalla & Marins 2020;Ghodsi Y. et al. 2022;Brandenberger 2022;Ghodsi Yengejeh et al. 2023).Therefore, to precisely characterize these dark sectors of the universe, numerous endeavors have been carried out to extend the Einstein-Hilbert action, resulting in what are referred to as modified theories of gravity.
One of the renowned modified theories of gravity is known as f (R) gravity, wherein a function of the Ricci scalar is substituted for the Ricci scalar itself within the Einstein-Hilbert action.To some extent, f (R) gravity has demonstrated the ability to address the characterization of dark sectors of the Universe.However, some studies have shown that general f (R) gravity faces challenges in simultaneously satisfying cosmological constraints and solar system tests of gravity (Gu 2011;Guo 2014;Negrelli et al. 2020).In this context, it has been demonstrated that certain categories of f (R) gravity theories offer a viable solution to this matter to a satisfactory degree.On of these models, known as the Hu-Sawicki model (Hu & Sawicki 2007), introduces a scalar field, known as the scalaron, which is derived from the modification of gravitational theory.An interesting property of the scalaron in this model is its chameleonic nature, whereby it becomes light in low-density regions and heavy in high-density regions of matter.The Hu-Sawicki model is one of the viable models within the framework of f (R) gravity and has been used to explain dark sectors of the Universe with higher accuracy, see, e.g., (Martinelli et al. 2009) and references therein.
As a result, enhancing the precision of theoretical forecasts concerning the dark sectors of the Universe holds the promise of yielding a more comprehensive picture of the abundance of PBHs, which are being considered as potential candidates for dark matter.Consequently, our aim is to compute the merger rate of PBHs within the framework of f (R) gravity.In this respect, the outline of the work is as follows: In Sec. 2, we describe the theoretical framework of f (R) gravity.In Sec. 3, we elucidate the concepts of dark matter halo models and essential factors, such as the halo density profile, the halo concentration parameter, and the halo mass function, within the theoretical contexts of both GR and f (R) gravity.Moreover, in Sec. 4, we specify the merger rate of PBHs in f (R) gravity and compare it with the corresponding findings from GR.We also constrain the abundance of PBHs in the framework of f (R) gravity while considering various masses of PBHs.Finally, in Sec. 5, we discuss the results and summarize the findings.of f (R) gravity is defined by (Lombriser et al. 2012) in which R stands for the Ricci scalar, κ represents the Einstein gravitational constant, g is the determinant of the metric tensor g µν , and S m denotes the action associated with matter fields ψ m .It is important to observe that by selecting f = −2Λ, GR can be reinstated along with a cosmological constant.It is worth emphasizing that different models can be discerned by the specific f (R) function chosen.Modifying the metric through variation yields the adjusted field equations for f (R) gravity.
(2) where G µν = R µν − 1/2Rg µν is defined as the Einstein tensor, and R µν represents the Ricci curvature tensor.The additional scalar degree of freedom, termed the scalaron and denoted by f R ≡ df (R)/dR, emerges from the theory and characterizes the alterations in forces.T µν stands for the energy-momentum tensor.The Levi-Civita-like connection is denoted as ∇ µ , and the d'Alembert operator is abbreviated as □ ≡ ∇ α ∇ α .It is evident that the aforementioned field equations involve derivatives up to the fourth order.Consequently, these equations can be interpreted as the standard field equations for GR supplemented by an additional scalar field f R .Furthermore, taking the trace of the modified field equations permits the derivation of the equation of motion for the scalaron where ρ m is the matter density of the Universe.Given the quasistatic approximation for a flat Universe and scales sufficiently below the horizon, one can neglect the time derivative of f R in the field equations.Hence, Eq.
(3) can be simplified as in which ρ m denotes the density of matter in the Universe.Assuming a quasistatic approximation for a flat Universe and focusing on scales significantly smaller than the horizon, one can disregard the temporal change of f R within the field equations.As a result, Eq. ( 3) can be rendered more straightforward as follows: which addresses the Newtonian potential Φ within the framework of f (R) gravity.Ensuring highly precise forecasts akin to those of GR within the solar system holds substantial importance, particularly when adjusting the gravitational theory to account for the observed accelerated expansion of the Universe on large scales.On the other hand, f (R) gravity within Jordan frame can be converted into its counterpart in Einstein frame through a conformal transformation applied to the metric (Lombriser et al. 2012) where A 2 (ϕ) is the conformal factor and a tilde demonstrates quantities in Einstein frame, and The integration of Eq. ( 8) yields the subsequent formulation for the scalar field Variation of the action with respect to ϕ results in: where α = d ln A/dϕ and V eff represents the effective potential dictating the behavior of ϕ.It is worth mentioning that T = A 4 (ϕ)T .In the quasistatic regime, one can disregard temporal derivatives in Eq. ( 12), leading to the scalar field equation In (Hu & Sawicki 2007), Hu and Sawicki proposed a general class of broken power-law models that can satisfy the aforementioned requirements.The specific functional form of f (R) suggested by Hu and Sawicki has the following form where m 2 = κρ m0 /3 defines the characteristic mass scale, and ρm0 = ρ(ln a = 0) represents the presentday background density of matter.The parameters c 1 , c 2 , and n > 0 are dimensionless free parameters that require specific determination to reproduce the expansion history and satisfy solar-system test via the chameleon mechanism (Khoury & Weltman 2004;Navarro & Van Acoleyen 2007;Faulkner et al. 2007).
It is important to select the sign of f (R) in a manner that ensures the second derivative of f (R) adheres to the subsequent condition: This requirement guarantees the stability of the solution in high-density regions when R greatly exceeds m 2 .Additionally, the presence of a nonzero and positive second derivative of the functional expression f (R) ensures the alignment of cosmological tests derived from the model with those originating from GR In regions characterized by significantly higher curvature compared to m 2 , the functional form of Hu-Sawicki model can be expanded as follows: lim Although the Hu-Sawicki model does not include a true cosmological constant, at constant c 1 /c 2 , the limiting case of c 1 /c 2 2 → 0 behaves as a cosmological constant at both large-scale and local expriments.Additionally, the finite c 1 /c 2 2 leads to the constant value of curvature, which remains unchanged as the matter density changes.As a result, by this choice one can have a class of models that are able to accelerate the expansion of the Universe similar to the behavior of the standard model of cosmology.Hence, one can rewrite Eq. ( 17) as follows: where R0 is the present-day background curvature, and f R0 ≡ f R ( R0 ).Additionally, by demanding similarity to standard model of cosmology as where ρΛ can be interpreted as the background energy density of dark energy.Eventually, from Eq. ( 10), one can obtain the following relations where 1+n) .The field strength f R0 controls the strength of the modification and is constrained by cosmological and solar-system tests (Desmond & Ferreira 2020).Different values of f R0 have been considered in the literature, ranging from 10 −4 to 10 −8 , e.g., (Mirzatuny & Pierpaoli 2019).In this work, we mainly focus on the Hu-Sawicki model with n = 1 and |f R0 | = 10 −4 , 10 −5 , and 10 −6 (labeled as f 4, f 5, and f 6, respectively).

Halo Density Profile
In the conventional understanding of cosmology, dark matter halos are regarded as non-linear structures that have been dispersed throughout the Universe due to the development and evolution of hierarchical structures.Initially, density fluctuations in the early Universe might have surpassed a critical threshold, leading to their collapse under the influence of self-gravitational forces, thereby making them capable of forming dark matter halos, see, e.g., (Genel et al. 2010;Ishiyama 2014;Del Popolo & Fakhry 2023).Physically speaking, these circumstances can be described by a dimensionless parameter called density contrast, which is derived from the excursion sets theory.The density contrast is defined as δ(r) ≡ [ρ(r) − ρ]/ρ, where ρ(r) is the density of the overdense region at any given point r, and ρ is the average density of the background.
In addition, the cosmological and structure formation models depend on the characteristics of the inner regions of dark matter halos.These halos' masses are governed by a radius-dependent function known as the density profile.To establish a reliable criterion for predicting the distribution of dark matter within galactic halos, various techniques such as spectroscopic observations of gravitational lensing, x-ray temperature maps, and stellar dynamics in galaxies have been utilized (Reed et al. 2005).Over the past decades, both analytical approaches and numerical simulations have been employed to derive an appropriate density profile that aligns with the observed data (Einasto 1965;Jaffe 1983;de Zeeuw 1985;Dehnen 1993;Navarro et al. 1996).One of the density profiles proposed based on N -body simulations within the framework of cold dark matter models is known as the Navarro, Frenk, and White (NFW) density profile (Navarro et al. 1996) where ρ s and r s are the scaled density and radius that vary from halo to halo.However, using analytical methods, Einasto discovered an alternative and appropriate definition for the density profile (Einasto 1965), which is as follows: whre α is the shape parameter for the Einasto density profile.
It is important to highlight that for both of the aforementioned expressions, one possesses which implies that the logarithmic slope of the density distribution must be −2 at the scaled radius.
In fact, the scaled density can be specified as ρ s = ρ crit δ c , where ρ crit represents the critical density of the Universe at a specific redshift z, while δ c denotes the linear threshold for overdensities, which relies on the concentration parameter C according to the following relation The concentration parameter essentially characterizes the central density of galactic halos, which is defined as the ratio between the virial radius of the halo, r vir , and its scale radius, r s .The halo virial radius covers a volume within which the average halo density is 200 to 500 times the critical density of the Universe.According to numerical simulations and analytical investigations, to have accurate predictions, the concentration parameter needs to dynamically evolve with mass and redshift (Prada et al. 2012;Ludlow et al. 2016;Okoli & Afshordi 2016;Okoli 2017).This aligns with the dynamics associated with the merging history of dark matter halos and their developmental pathways.This alignment stems from the fact that smaller halos have already virialized, resulting in a higher degree of concentration compared to the larger ones.Nevertheless, it has been demonstrated that within the context of f (R) gravity, the concentration parameter is influenced not only by mass and redshift but also by the scalar parameter f R0 (Mitchell et al. 2019).In this study, we utilize the concentration parameter formulated in Okoli & Afshordi (2016) for our computations within the framework of GR, while adopting the concentration parameter obtained from Mitchell et al. (2019) for our calculations involving f (R) gravity.In the next section, we will delve into the statistical characteristics of dark matter halos concerning the halo mass function.

Halo Mass Function
The presence of dark matter halos provides a practical and fundamental way to examine the non-linear gravitational collapse in the Universe.Hence, gaining a proper statistical perspective on the mass distribution of these halos can enhance our understanding of the physics governing them.Consequently, the halo mass function has been proposed as a way to describe the mass distribution of these structures within a given volume.In simple terms, the halo mass function explains the masses of these structures that have densities higher than a certain threshold, unaffected by the Universe's expansion, and destined to collapse.As the Universe expands, the density contrast can grow to a critical point, surpassing linear regimes and entering nonlinear regimes.At this stage, overdensities detach from the expansion of the Universe, enter the turnaround phase, and undergo collapse, leading to the formation of structures.
In the standard model of cosmology based on GR, the halo mass function can be derived analytically using the framework of excursion set theory, which models the density field as a stochastic process across scales.The fundamental premise of excursion set theory is the spherical collapse model, which determines the threshold overdensity required for collapse of a spherical perturbation (Press & Schechter 1974).In the case of an Einstein-de Sitter Universe and a spherical-collapse halo model, the threshold overdensity can be calculated as follows (Nakamura & Suto 1997) (26) where Ω m represents the density parameter of matter content.Consequently, one can approximate the threshold overdensity as 1.686 within a narrow redshift range.
A key feature of the present analysis is that δ sc is approximately independent of halo mass M due to Birkhoff's theorem, which states that the evolution of a spherical density profile is oblivious to external influences.
In Jenkins et al. (2001), a suitable definition of the differential mass function has been introduced to specify different fits for dark matter halos In this context, n(M, z) represents the number density of halos with a mass M at redshift z, ρ m denotes the cosmological matter density, and g(σ) is multiplicity function that relies on the geometry of overdensities at the time of collapse.The function σ(M, z) signifies the linear root mean square fluctuation of overdensities on mass M and at redshift z, precisely defined as where W (k, M ) represents the Fourier spectrum of the top-hat filter, depending on the mass M and wavenumber k.Additionally, P (k, z) stands for the redshiftdependent power spectrum of the density fluctuations.Numerous investigations have been carried out to estimate the halo mass function through analytical methods and numerical simulations.The objective of these studies is to find the most accurate representation for cosmic observations.By incorporating a homogeneous and isotropic collapse in the standard excursion sets theory, a straightforward analytical expression for the multiplicity function emerges as follows (Press & Schechter 1974) which is called the Press-Schechter (PS) mass function.The above calculation pertains to the collapse of spherically-symmetric overdensities, but the real Universe involves a much more intricate scenario.The collapse dynamics are not spherical; instead, they are triaxial, and small overdense regions need additional matter to collapse due to significant influences from the surrounding shear field (Sheth et al. 2001).To address this complexity, the excursion set approach adopts ellipsoidal collapse, which introduces a stochastic barrier, prompting the investigation of a general barrier.
In the case of the standard model of cosmology, a straightforward Gaussian distribution for the barrier B, with a mean value B that linearly changes with the variance S ≡ σ 2 (M, z), proves to be highly accurate in replicating the N -body halo mass function (Corasaniti & Achitouv 2011a,b;Achitouv & Corasaniti 2012).Moreover, this barrier aligns with the measured overdensity required for collapse in the initial conditions (Achitouv et al. 2013), and it offers the advantage of having an exact solution for the Markovian multiplicity function.Consequently, the precise solution for a fixed and linear barrier can be expressed as follows: where a is defined as a = 1/(1 + D B ), and D B is a parameter that characterizes the diffusive nature of the barrier B(S) (Kopp et al. 2013).
As evident from Eq. ( 29), at a fixed redshift, the mass function relies solely on the halo mass through σ(M ), and significant changes are not anticipated to occur.This mass function represents a simple model proposed for the formation of dark matter halos, known as the spherical collapse model, which often aligns with observational data.However, it exhibits quantitative deviations from numerical results at certain mass limits (Jenkins et al. 2001).
To address this issue, one can approximate the precise solution for a generic barrier by expanding it in higherorder derivatives of the barrier.One notably successful improvement was presented by Sheth and Tormen, which is based on a more realistic model and provides a better fit to simulation results (Sheth & Tormen 1999;Sheth et al. 2001).Their approach adopts an ellipsoidal collapse model with dynamical threshold density fluctuations, in contrast to the nearly global threshold used in the PS model.Actually, Sheth and Tormen introduced the idea of considering a dynamically varying threshold overdensity for ellipsoidal collapses, denoted as δ ec , which provides a more realistic depiction of the halo mass function.By assuming zero prolateness, they calculated this quantity as where γ = 0.47, β = 0.615, and defining ν ≡ δ sc /σ(M ).
It becomes evident that this quantity depends not only implicitly on the redshift but also on the mass of the structure, and it is referred to as the moving barrier.
Under this assumption, one can obtain the halo mass function for ellipsoidal collapse, also known as the Sheth-Tormen (ST) mass function, which is given by where a = 0.3222, b = 0.707, and p = 0.3.It is anticipated that this mass function will exhibit greater sensitivity to changes in redshift compared to the PS mass function.
However, in modified gravity theories, the story becomes far more complex.For example, in f (R) gravity models, which introduce an additional scalar degree of freedom, Birkhoff's theorem is violated.This makes the collapse process dependent on the environment, which is known as the chameleon screening mechanism, wherein dense regions suppress the enhanced gravitational forces to some extent.Consequently, the threshold of overdensities becomes dependent on halo mass, redshift, and specific shape of the model parameter f R0 , which alters the scalar field gradients (Kopp et al. 2013) where To account for these intricacies, the halo mass function in f (R) gravity must be derived through numerical solutions of the fully nonlinear modified Einstein equations, using initial conditions based on peaks theory rather than simplistic analytical approximations.The resulting δ c (z, M, f R0 ) can then be incorporated into the excursion set framework, along with a drifting diffusive barrier to model the stochasticity of realistic triaxial collapse.In this way, the complex environmental couplings between the scalar field gradient and local density configurations are captured.This enables accurate predictions of the halo mass function in f (R) gravity, which can exhibit significant deviations from the GR case and thus provide a signature to constrain these modified gravity models.Under there assumptions, the multiplicity function for f (R) gravity is given by (Kopp et al. 2013) where g f (R),sk (σ) and g GR,sk (σ) represent multiplicity functions applying the sharp-k filter for f (R) gravity and GR, respectively, and their definitions are provided as follows: and The value of a ≃ 0.714 has been considered in the method developed in Kopp et al. (2013).In above relations, collapse barriers have the following form where β characterizes the linear drift of the barrier height as a function of the variance S. The value of β needs to be determined by calibrating to N -body simulations or models of ellipsoidal collapse.Typical values are β ≃ 0.1 − 0.2 for a standard model of cosmology.In the context of f (R) gravity, the drift parameter β would likely be slightly modified from its GR value, since the dynamics of aspherical collapse can be altered.However, as a first approximation, one can assumes the GR value of β in deriving the f (R) halo mass function.Moreover, g GR,sx (σ) refers to the multiplicity function calculated in the context of GR using a sharp-x filtering method, which establishes a connection between halo mass and variance.Its complete expression is (Kopp et al. 2013) where In the above relations, κ ≃ 0.65 S/δ 2 sc , Γ(x, y) is the incomplete gamma function, and Erfc(x) is the complementary error function.By following all the outlined procedures, the halo mass function can be derived within the framework of the f (R) gravity.Clearly, this mass function relies not only on mass and redshift but also on the field strength f R0 .Consequently, variations in f R0 are anticipated to lead to changes in the halo mass function in f (R) gravity.

MERGER RATE OF PBHS
As previously mentioned, PBHs are distributed randomly throughout the Universe, which means they can encounter each other and form binary systems.In this section, our objective is to compute the merger rate of PBHs within the context of f (R) gravity and compare the results with those obtained from GR.To achieve this, we have previously discussed the essential tools for modeling dark matter halos.
Let us consider two PBHs with masses m 1 and m 2 and a relative velocity at large separation v rel = |v 1 −v 2 | that coincidentally encounter each other within a dark matter halo.Under these conditions and according to Keplerian mechanics, the maximum gravitational radiation is emitted at the periastron r p .If the emitted gravitational energy is greater than the kinetic energy of the system, PBHs will become gravitationally bound and form binary systems.As a consequence of this situation, there is an upper limit to the periastron distance (Peters 1964) where G represents the gravitational constant and c stands for the velocity of light.Moreover, within the Newtonian limit, a connection between the impact parameter and the periastron is indicated by the following relation (Sasaki et al. 2018) When considering scenarios where the tidal forces of surrounding black holes rarely result in head-on collisions, it can be expected that the gravitational interaction between the two PBHs leads to the formation of a binary system with maximum eccentricity.On the other hand, under the strong gravitational focusing limits, it can be roughly implied that the tidal forces of the surrounding black holes cannot disturb the orbital parameters of the formed binaries.
Note that PBH binaries forming within dark matter halos have merger times ranging from hours to kiloyears, with dissipative two-body encounters resulting in much shorter merger times than the age of the Universe.On the other hand, binaries formed through nondissipative three-body encounters typically have longer merger times, making them less likely to significantly contribute to the population of BH mergers observed by LIGO-Virgo-KAGRA detectors.Nevertheless, there might be specific circumstances where the merger rate of binaries formed through three-body encounters is likely significant.In this work, we mainly focus on those binaries that can be formed through two-body encounters.Therefore, the cross-section for binary formation can be computed using the following relation (Quinlan & Shapiro 1989;Mouri & Taniguchi 2002) (52) As a result, when Eq. ( 50) is inserted into Eq.( 52), it leads to an explicit expression for the cross-section for the binary formation This relation is derived using strong limit gravitational focusing, where the tidal forces from neighboring compact objects on the binary system can be disregarded, i.e., r p ≪ b.
In this study, we focus on events that align with the sensitivity of LIGO-Virgo-KAGRA detectors.Therefore, our analysis is limited to cases where m 1 = m 2 = M PBH .With this constraint in mind, the binary formation rate within a galactic halo can be calculated as follows (Bird et al. 2016) where 0 < f PBH ≤ 1 indicates the proportion of PBHs contributing to dark matter.Additionally, ρ halo denotes the halo density profile, and the angle bracket represents an average over the distribution of PBH relative velocities within the galactic halo.Furthermore, the halo virial mass, representing the mass enclosed within the virial radius, is determined using the following equation which can be computed, taking into account the NFW and Einasto density profiles.The velocity dispersion of dark matter particles in galactic halos is another essential element in computing the merger rate of PBHs.In this regard, it has been revealed that the velocity dispersion profiles in f (R) gravity are similar to those in the standard model of cosmology (He et al. 2015).Thus, it is feasible to utilize, with reasonable accuracy, the velocity dispersion relation of dark matter particles obtained in Prada et al. (2012) for both f (R) gravity and GR where v max represents the maximum velocity within a radius of r max .Additionally, we require that the Maxwell-Boltzmann distribution, truncated at the virial velocity, governs the probability distribution function of relative velocities among PBHs in the galactic halos ) where J 0 is determined by fulfilling the requirement that 4π In the realm of PBH binary formation, two distinct mechanisms come into play, coexisting harmoniously in different epochs of the Universe (Sasaki et al. 2018).This study delves into the formation of PBH binaries within dark matter halos during the late-time Universe.Theoretical forecasts from this process even propose the intriguing possibility that the vast majority of dark matter might consist of PBHs (Bird et al. 2016;Fakhry et al. 2021Fakhry et al. , 2022a;;Fakhry & Del Popolo 2023).
However, during the early Universe, the initial clustering of PBHs could lead them to break away from the Hubble flow, giving rise to binary formations (Ali-Haïmoud et al. 2017).These PBH binaries, emitting GWs constantly, gradually reduce in size and eventually merge.Nonetheless, the merger process is not without its challenges, as tidal forces from surrounding PBHs might disrupt some of these binaries before they come to merge (Kavanagh et al. 2018;Raidal et al. 2019).
The merger time for PBH binaries formed in the early Universe hinges on their orbital parameters, which, due to the random distribution of PBHs during that era, can vary significantly at the time of their binary formation.Consequently, a diverse array of scenarios unfolds: some PBH binaries have already merged, others are destined to merge in the current Universe, and still others await their fate in the future.It is in this context that LIGO-Virgo-KAGRA observations can be elucidated, showcasing how PBHs constitute only a tiny fraction of dark matter (Hall et al. 2020;Hütsi et al. 2021;Chen et al. 2022).It is crucial to note that the journey of PBHs in the cosmic symphony continues to enthrall scientists, paving the way for captivating discoveries yet to come.Hence, both mechanisms remain valid for describing black hole binary mergers, and the predictions concerning the abundance of PBHs are still under scrutiny by the LIGO-Virgo-KAGRA detectors.In this context, we minimally assume the contribution of PBHs to dark matter as f PBH = 1 when calculating the merger rate of PBHs.I also set the mass of PBHs to be M PBH = 30 M ⊙ .However, we will later extend our analysis to encompass the changes in both mass and fraction of PBHs.
In Fig. 1, we have presented the merger rate of PBHs within each halo as a function of the halo mass for three models of f (R) gravity, i.e., f 4, f 5 and f 6, and compared them with the merger rate obtained for GR while considering the Einasto density profile.It is clear from the figure that the merger rate of PBHs within each halo, in comparison to the outcome derived from GR, is higher for all the examined models under f (R) gravity.This can be attributed to the effect of nonlinear dynamics considered for density fluctuations, δ c (z, M, f R0 ), and formation and evolution of halo structures in f (R) gravity, which appears in all stages of calculations.An additional intriguing finding emerges as we observe a gradual reduction in the influence of the field strength-dependent dynamics in density fluctuations, moving from f 4 to f 6.This can be attributed to the fact that the merger rate of PBHs, as dark matter candidates, is directly shaped not only by the dynamics of density fluctuations but also by the field strength f R0 .
Furthermore, the LIGO-Virgo-KAGRA detectors are designed to record and process the accumulated merger rate of black holes.Therefore, it is necessary to calculate the overall merger rate of PBHs per unit volume and time.To achieve this, one needs to combine the halo mass function, dn/dM vir , with the merger rate per halo, Γ(M vir ), The significance of the upper limit of integration in determining the merger rate of PBHs can be easily dismissed.This is due to the decreasing nature of the halo mass function, which causes the contribution of PBH merger rate to diminish exponentially as the halo mass increases.This behavior aligns with the hierarchical dynamics of halo formation, where low-mass halos have higher dark matter density than high-mass halos because they have already undergone virialization.As a result, the lower limit of integration holds greater importance in this analysis.Referring to arguments presented in Fakhry et al. (2021Fakhry et al. ( , 2022a)), when calculating  the merger rate of PBHs with a mass of 30 M ⊙ , the lower limit for halo mass needs to be set at M c ≃ 400 M ⊙ .This implies that signals from dark matter halos with masses below 400 M ⊙ are expected to be negligible.Additionally, in this analysis, we employ Eq. ( 32) to represent the halo mass function for GR, while we use Eq. ( 40) for the halo mass function in the context of f (R) gravity.
In Fig. 2, we have displayed the merger rate of PBHs per unit time and volume in three f (R) gravity models, i.e., f 4, f 5 and f 6, and compared them with the corresponding results obtained from GR while taking into account the Einasto density profile.These calculations are performed for the present-time Universe.As it is evident from the inset figure, smaller-mass halos continue to play a more significant role in driving the merger rate of PBHs compared to the larger halos.This is a direct consequence of the higher density of dark matter particles within subhalos, as previously discussed, creating a higher concentration in subhalos compared to the host halos.It is important to highlight that the merger rate of PBH binaries can be precisely quantified through the integration across the area below the curves.Furthermore, it can be deduced that across all the considered models within f (R) gravity, the merger rate of PBHs surpasses that obtained from GR. Also, the direct dependence of the merger rate of PBHs with the field strength f R0 is evident.Evidently, as the field strength weakens, the merger rate of PBHs for f (R) gravitaty progressively converge towards the findings obtained through GR.This implies that by constraining the value of field strength and comparing the predictions from the present analysis with GW data, it is possible to potentially constrain the abundance of PBHs, which in turn introduces a novel method.
One of the primary purposes of this analysis is to compute the merger rate of PBHs in the distant Universe, allowing for a meaningful comparison with the merger events detected by the LIGO-Virgo-KAGRA observatories.This endeavor aims to offer theoretical forecasts based on the gravitational models proposed in this research, shedding light on the forthcoming landscape of GW phenomena and the advancement of these detection instruments.Through our analysis, we aim to illustrate the redshif evolution of the merger rate of PBHs.On the other hand, with the present sensitivity of GW detectors, they can detect merger events occurring within a comoving volume of 50 Gpc 3 , equivalent to a redshift of approximately z ∼ 0.75 (Abbott et al. 2021;The LIGO Scientific Collaboration et al. 2021).This raises an intriguing question of computing the redshift evolution of the merger rate of PBHs.It should be emphasized that Eq. ( 58), which incorporates the halo mass function and the concentration parameter, is sensitive to changes in redshift.To address this objective, we have presented in Fig. 3 the redshift evolution of the merger rate of PBHs for three f (R) gravity models, i.e., f 4, f 5 and f 6.Additionally, we have compared these results with the corresponding outcomes for GR while considering both NFW and Einasto density profiles.The direct relation between the merger rate PBHs and redshift is evident.This can be explained by considering the impact of hierarchical dynamics and the structure of halo mergers.This suggests that during earlier periods char- .Total merger event rate of PBHs with respect to their fraction for three models of f gravity, i.e., f 4, f 5, and f 6, while taking into account the PBH masses to be MPBH = 10, 30, and 100 M⊙.The shaded band represents the total merger rate of black holes, as estimated by the LIGO-Virgo-KAGRA detectors during the second half of the third observing run (O3b), within the range of (17.9-44)Gpc −3 yr −1 .Both NFW and Einasto density profiles have been incorporated.
acterized by higher redshifts, there may have been a greater abundance of subhalos.Consequently, in higher redshifts, PBHs merged at a more accelerated rate compared to the present-time Universe.The findings indicate that the redshift evolution of the merger rate of PBHs within f (R) gravity models exceeds that obtained from the framework of GR.This implies that, assuming the credibility of f (R) gravity, the merger rate of PBHs will be amplified over time.Furthermore, it is evident that the impact of the field strength f R0 on enhancing the merger rate of PBHs will endure during the late-time Universe.
Hereafter, our attention will be focused on the expected PBH fraction, f PBH , which comes from the f (R) gravity models.The issue of the fraction of PBHs has been a significant concern from the outset of the development of the PBH scenario.Furthermore, a crucial constraint placed on PBHs pertains to their abundance within the Universe during late times.Nowadays, the abundance fractions of many mass ranges of PBHs have been strongly constrained using observational and nonobservational methods (Carr et al. 2021(Carr et al. , 2023)).However, a specific mass interval of PBHs, referred to as asteroid mass PBHs, i.e., those with masses around 10 −17 M ⊙ ≤ M PBH ≤ 10 −12 M ⊙ (Montero-Camacho et al. 2019;Smyth et al. 2020;Ray et al. 2021;Ghosh & Mishra 2022), has yet to be definitively constrained, and could potentially account for a substantial portion of the dark matter composition.On the othe hand, comparing the PBH merger rate obtained from each theoretical model with the estimated one via the LIGO-Virgo-KAGRA detectors can potentially be one of the best references for validating that model.In Fig. 4, based on NFW and Einasto density profiles, we have plotted the total merger rate of PBHs as a function of f PBH for the models of f (R) gravity and compared the results with those extracted from GR.Additionally, the shaded band represents the total merger rate of black holes estimated by the LIGO-Virgo-KAGRA detectors during the latest observing run, O3b, which is (17.9-44)Gpc −3 yr −1 (The LIGO Scientific Collaboration et al. 2021). is evident from the figure that the merger rate of PBHs for both gravitational models is inversely proportional to their masses but directly proportional to their fractions.This is due to the fact that the number density of PBHs changes inversely with their masses.It can also be inferred that, compared to the corresponding results derived from GR, the models of f (R) gravity satisfy the constraints stemming from GW data for relatively lower values of the fraction of PBHs.In addition, the direct effect of field strength f R0 from f (R) gravity in imposing stringent constraints on the fraction of PBHs is evident in the present analysis.
Up until now, we have assumed in our analysis that the mass of the involved PBHs is 30 M ⊙ and that they can contribute maximally to dark matter.However, it is of interest to calculate the merger rate of PBHs based on different assumptions regarding their fractions and masses.To this end, in Fig. 5, we have depicted the total merger rate of PBHs as a function of 10, M ⊙ ≤ M PBH ≤ 100, M ⊙ and f PBH for f (R) gravity.In this calculation, NFW and Einasto density profiles have been incorporated, considering three values of field strength: f 4, f 5, and f 6.Once again, it is evident that the merger rate of PBHs is inversely proportional to their masses.In other words, the smaller the mass of PBHs participating in the merger event, the higher their number density per unit volume, and consequently, their total merger rate.As a result, the theoretical framework of f (R) gravity imposes more stringent constraints on the fraction of PBHs if smaller masses of PBHs are considered.Additionally, it can be deduced from the figure that the field strength value f R0 is also an actively contributing factor in constraining the fraction of PBHs.This establishes a direct relationship between the field strength and the stringency of the constraints on the fraction of PBHs.In this regard, the most stringent constraints can be obtained from f (R) models of gravity with field strengths f 4, f 5, and f 6, respectively.We have also quantified the overall results of our analysis for the merger rate of PBHs within the frameworks of f (R) gravity and GR in Table 1.The results show that the merger rate of PBHs, while considering the Einasto density profile, are on average about 60% higher than the results obtained from the NFW density profile.Furthemore, it can be inferred that the most stringent constraint in this analysis is obtained from the models of f (R) gravity, while considering the Einasto density profile and M PBH ≃ 10 M ⊙ , which, despite all theoretical uncertainties, enters the LIGO-Virgo-KAGRA sensitivity band for f PBH ≳ 0.1.

CONCLUSIONS
Primordial black holes are considered one of the most mysterious phenomena in astrophysics, and fundamental questions about their nature continue to be raised.As PBHs are expected to interact solely through gravity, and considering that a substantial collection of black holes demonstrates characteristics of perfect fluids on significantly large scales, PBHs present themselves as natural candidates for dark matter.On the other hand, due to their random distribution in the Universe, PBHs have the possibility of encountering each other, forming binaries, and eventually merging through the continu-Table 1.Total merger rate of PBHs in the context of GR and f (R) gravity with field strengths f 4, f 5, and f 6 as a function of different PBH masses, i.e., MPBH = 10, 20, 30, 50, and 100 M⊙, while considering NFW and the Einasto density profiles, at the present-time Universe (z = 0).ous propagation of GWs in the medium of dark matter halos.The dynamics of PBHs as dark matter candidates are expected to be influenced by the local and statistical properties of dark matter halos.However, a fundamental challenge arises as to whether dark matter halo models based on GR are good enough to accurately predict the merger rate of PBHs.To address this question, in this study, we have calculated the merger rate of PBHs within the framework of f (R) gravity and compared it with the corresponding results obtained from GR.To accomplish this task, we have initially established an appropriate framework for dark matter halo models that suits both GR and f (R) gravity.We have introduced the definition of key parameters including the halo density profile, the halo concentration parameter, and the halo mass function.We have also discussed the field strength-dependent dynamical conditions induced by f (R) gravity and introduced the appropriate density contrast parameter, δ c (z, M, f R0 ), and the mass function, g f (R),sx (σ), for f (R) gravity.

PBH Mass (M⊙) Density Profile RGR(Gpc
Furthermore, we have investigated the encounter conditions of randomly distributed PBHs within the context of dark matter halos.Under these assumptions, and considering M PBH = 30 M ⊙ and f PBH = 1, we have calculated the merger rate of PBHs per halo while considering three models of f (R) gravity, and compared the results with those obtained under GR.The results indicate that, in comparison to the outcome derived from GR, the merger rate of PBHs within each halo is higher for all the examined models under f (R) gravity.This is because of the field strength-dependent dynamics of density fluctuations, δ c (z, M, f R0 ), which affects the formation and evolution of halo structures under f (R) gravity.Additionally, we have witnessed a gradual reduction in the influence of field strenght-dependent dynamics, while moving from f 4 to f 6.
Based on the PBH scenario and suitable halo mass functions, we have calculated the merger rate of PBHs per unit time and volume for three models of f (R) gravity, and qualitatively and quantitatively compared them with the corresponding results obtained from GR.The results demonstrate that smaller-mass halos continue to exert a more significant influence on the merger rate of PBHs compared to larger halos.This phenomenon directly stems from the higher density of dark matter particles within subhalos, leading to a greater concentration within subhalos than in the host halos.Moreover, it can be deduced that across all the examined models under f (R) gravity, the merger rate of PBHs surpasses what is obtained from GR.The direct connection between the merger rate of PBHs and the field strength is also evident.
The potential for binary PBH formation throughout the age of the Universe, stemming from their random distribution, serves as strong motivation to investigate the evolution of the PBH merger rate as a function of redshift.In light of this, we have specified the redshiftevolution of the PBH merger rate for three models of f (R) gravity, and compared these findings with that obtained from GR. Consequently, the results demonstrate a direct correlation between the total merger rate of PBHs and the redshift in both models.In simpler terms, PBHs have exhibited a greater tendency to form binaries at higher redshifts compared to the present-time Universe.Furthermore, the direct correlation between the merger rate of PBHs and redshift becomes evident.This phenomenon can be elucidated by considering the influence of hierarchical dynamics of halo structures.The findings indicate that the redshift-dependent evolution of the PBH merger rate within f (R) gravity models surpasses that obtained from the framework of GR.This suggests that, if we assume the validity of f (R) gravity, the merger rate of PBHs will increase over time.Furthermore, it becomes evident that the influence of the field strength in enhancing the PBH merger rate will persist throughout the late-time Universe.
Lastly, we have computed the PBH merger rate for GR and three f (R) gravity models as a function of their fraction with M PBH = 30 M ⊙ , and then compared these results with the black hole mergers estimated by the LIGO-Virgo-KAGRA detectors during the latest observing run, i.e., (17.9-44)Gpc −3 yr −1 .The results indicate that the merger rate of PBHs is inversely proportional to their masses, yet directly proportional to their fractions in both gravitational models.This phenomenon arises due to the inverse variation of PBH number density with their masses.Furthermore, it can be inferred that, when compared to the corresponding results derived from GR, the f (R) gravity models satisfy the constraints imposed by GW data, particularly for lower values of the PBH fraction.Moreover, the significant impact of the field strength f R0 , originating from f (R) gravity, in imposing stringent constraints on the PBH fraction becomes evident in our analysis.I have also provided a separate calculation for the merger rate of PBHs as a function of their mass and fraction in the context of f (R) gravity.The findings indicate that, when the Einasto density profile is taken into account, the merger rate of PBHs is, on average, roughly 60% higher than the findings from the NFW density profile.Fur-thermore, it can be deduced that the f (R) gravity models, taking into account the Einasto density profile and M PBH ≃ 10 M ⊙ , yield the most stringent constraints.This mass range penetrates the LIGO-Virgo-KAGRA sensitivity region for f PBH ≳ 0.1 despite all theoretical uncertainties.
It should be noted that constraints on PBHs are subject to a variety of uncertainties, encompassing various gravitational frameworks, conditions that may have been imposed on the structures during collapse and formation, specific processes that can impact the growth or evaporation of PBHs (such as accretion and merger history), uncertainties arising from black hole formation scenarios, and their contribution to LIGO-Virgo-KAGRA mergers, as we have observed.While the presence of these factors might lead to computational errors, the development of future instruments and a deeper understanding of these unidentified processes may ultimately result in more stringent constraints on the fraction of PBHs.

Figure 3 .
Figure 3. Redshift evolution of the total merger rate of PBHs per unit source time and per unit comoving volume for (top) NFW and (bottom) Einasto density profiles.Solid, dashed, and dot-dashed lines exhibit this relation for dark matter halos in f (R) gravity, with |fR0| = 10 −4 , 10 −5 , and 10 −6 , respectively.While dotted lines indicate it for dark matter halos in GR.

Figure 4 .
Figure4.Total merger event rate of PBHs with respect to the PBH fraction for (top) NFW and (bottom) Einasto density profiles.Solid, dashed, and dot-dashed lines exhibit this relation for dark matter halos in f (R) gravity, with |fR0| = 10 −4 , 10 −5 , and 10 −6 , respectively.While the dotted line indicates it for dark matter halos in GR.The mass of involving PBHs has been considered as MPBH = 30M⊙.Also, The shaded band represents the total merger rate of black holes estimated by the LIGO-Virgo-KAGRA detectors during the second half of the third observing run (O3b), i.e., (17.9-44)Gpc −3 yr −1 .
Figure5.Total merger event rate of PBHs with respect to their fraction for three models of f gravity, i.e., f 4, f 5, and f 6, while taking into account the PBH masses to be MPBH = 10, 30, and 100 M⊙.The shaded band represents the total merger rate of black holes, as estimated by the LIGO-Virgo-KAGRA detectors during the second half of the third observing run (O3b), within the range of (17.9-44)Gpc −3 yr −1 .Both NFW and Einasto density profiles have been incorporated.
The merger rate of PBHs within each halo as a function of halo mass for Einasto density profile.Solid, dashed, and dot-dashed lines exhibit this relation for dark matter halos in f (R) gravity, with |fR0| = 10 −4 , 10 −5 , and 10 −6 , respectively.While the dotted line indicates it for dark matter halos in GR.