Primordial Black Hole–Neutron Star Merger Rate in Modified Gravity

In this work, we investigate the merger rate of primordial black hole–neutron star (PBH-NS) binaries in two widely studied modified gravity (MG) models: Hu–Sawicki f(R) gravity and the normal branch of Dvali–Gabadadze–Porrati gravity. In our analysis, we take into account the effects of MG on the halo properties, including halo mass function, halo concentration parameter, halo density profile, and velocity dispersion of dark matter particles. We find that these MG models, due to their stronger gravitational field induced by an effective fifth force, predict enhanced merger rates compared to general relativity. This enhancement is found to be redshift-dependent and sensitive to model parameters and PBH mass and fraction. Assuming a PBH mass range of 5–50 M ⊙, we compare the predicted merger rate of PBH-NS binaries with those inferred from LIGO–Virgo–KAGRA observations of gravitational waves (GWs). We find that the merger rates obtained from MG models will be consistent with the GW observations if the abundance of PBHs is relatively large, with the exact amount depending on the MG model and its parameter values, as well as PBH mass. We also establish upper limits on the abundance of PBHs in these MG frameworks while comparing them with the existing non-GW constraints, which can potentially impose even more stringent constraints.


INTRODUCTION
For many years, gravitational waves (GWs) have been a focal point in the field of cosmology, offering a novel approach to explore a variety of astrophysical and cosmological events.The merger of compact binary systems is considered as a potential source of GW emission (Mandel & Broekgaarden 2022).In recent years, GW detectors have documented numerous events resulting from compact binary mergers (e.g., Abbott et al. 2016aAbbott et al. ,b,c, 2020a,b),b).Compact binaries are typically classified into three categories: binary black holes (BBH), black holeneutron star (BH-NS) binaries, and binary neutron stars (BNS).Among these, most recorded GWs are attributed to BBH mergers in the mass range of 10-100 M ⊙ (e.g., Abbott et al. 2019Abbott et al. , 2021aAbbott et al. , 2023)).While the genesis s fakhry@sbu.ac.ir ma shiravand@sbu.ac.ir m farhang@sbu.ac.ir of blach holes (BHs) has been thoroughly investigated, the precise process of their formation remains an active area of study (e.g., Raidal et al. 2019;Bouffanais et al. 2021;Nitz & Wang 2021a,b;Mandel & Farmer 2022;Lin et al. 2023).They might have originated from stellar collapse through various channels (Fishbach et al. 2021;Rodriguez et al. 2021), resulting in what is known as astrophysical black holes (ABHs).Alternatively, they could be remnants from the early Universe, where density fluctuation-induced gravitational collapse could have led to the formation of primordial black holes (PBHs) (Zel'dovich & Novikov 1967;Hawking 1971).
PBHs, with their diverse masses (e.g., Sasaki et al. 2018), are thought to stand out from ABHs and are considered potential dark matter candidates (for other possible candidates, see, e.g., Spergel & Steinhardt 2000;Alcock et al. 2000;Liebling & Palenzuela 2012;Roszkowski et al. 2018;Di Giovanni et al. 2020, 2021).Over the years, robust observational techniques have been used to constrain the abundance of PBHs across various mass scales.These observations have provided reliable frameworks for studying the early Universe at small scales (e.g., Carr et al. 2017;Lehmann et al. 2018;Carr et al. 2021).Furthermore, the presumption of their involvement in merger events linked to GW detectors allows for the acquisition of robust limitations on the abundance of PBHs.In spite of numerous theoretical ambiguities, PBHs of stellar mass (those in line with GW observations) might play a substantial role in the composition of dark matter (Bird et al. 2016;Clesse & García-Bellido 2017).Nevertheless, the hypothesis of binary PBH genesis in the early Universe implies that the contribution of such Black Holes to the structure of dark matter ought to be minimal to align with the observations of the LIGO-Virgo-KAGRA (LVK) detectors (Sasaki et al. 2016).
Moreover, the LVK detectors are able to identify a distinct kind of GW event referred to as BNS merger events (Ye et al. 2020).BNS systems can dynamically originate in high-density areas like star clusters.These events emit not only GWs but also electromagnetic waves, rendering them an important topic for multi-messenger astronomy research.It is important to highlight that distinguishing between NS candidates and solar-mass BHs poses a challenge without separate analysis of electromagnetic signals and setting suitable upper bounds on the tidal deformability parameter of BNS (Abbott et al. 2017a,b).This issue has been tackled in prior research, such as the studies conducted by Abbott et al. (2017a,b).Currently, the field of multi-messenger astronomy, which explores this problem, is leading the research in this domain (e.g., Tsai et al. 2021;Dasgupta et al. 2021).
Furthermore, BH-NS binary mergers can make a significant contribution to the data collected by the LVK detectors (Ruiz et al. 2021).Such events can yield valuable insights for multi-messenger observations and can generate both GWs and electromagnetic signals during the merger phase (Barbieri et al. 2020).Typically, these events involve a post-merger phase where the residual matter from the NS is accreted by the BH, leading to a luminous event (Fragione 2021).BH-NS binary mergers are of particular interest as they provide unique perspectives into the nuclear equation of state of NSs and the accretion processes of BHs (see, e.g., Hinderer et al. 2019;Zappa et al. 2019;Fragione et al. 2021;Tiwari et al. 2021).Moreover, they can aid in constraining the spin and abundance of BHs.Recently, the LVK detectors reported the first two direct events related to BH-NS mergers.The estimated mass components of these events were (8.9 +1.2 −1.5 M ⊙ , 1.9 +0.3 −0.2 M ⊙ ) and (5.7 +1.8  −2.1 M ⊙ , 1.5 +0.7 −0.3 M ⊙ ) (Abbott et al. 2021b).The formation process of BH-NS binaries is still fraught with many uncertainties.As previously mentioned, one possible assumption is that the contributing BHs are of primordial origin.However, unlike the formation of PBH binaries, the formation of PBH-NS binaries is associated with the late-time Universe, taking place immediately after the formation of cosmological and astrophysical structures.As a result, it is anticipated that more precise models describing dark matter halos will significantly influence the density and velocity distribution of PBHs involved in these events (see, e.g., Fakhry et al. 2021Fakhry et al. , 2022aFakhry et al. ,b, 2023c;;Fakhry & Del Popolo 2023;Fakhry et al. 2023b;Fakhry 2024).Various methods have been utilized to elaborate on halo characteristics within the context of GR (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;Del Popolo & Fakhry 2023;Fakhry et al. 2023a).However, a critical question arises as to whether the principles of GR and their derived concepts adequately explain the formation and evolution of large-scale structures such as galactic halos.Currently, there is growing consensus that while GR has been successful in accurately predicting numerous events, it falls short in fully elucidating large-scale phenomena, especially within the dark sectors of the Universe (for more details on see, e.g., Yang & Xu 2014;Huterer & Shafer 2018;Ghodsi Y. et al. 2022;Ghodsi Yengejeh et al. 2023).Consequently, considerable efforts have been devoted to generalizing the Einstein-Hilbert action, known as modified gravity (MG) (see.e.g., Nojiri & Odintsov 2007;Clifton et al. 2012;Saridakis et al. 2021).
To fulfill the observational criteria related to gravity, such as restoring GR in extensively tested highdensity scenarios and aligning with the expansion history (Ishak 2019), MG theories necessitate additional screening mechanisms (e.g., Aviles et al. 2019).The Chameleon and Vainshtein effects are two commonly studied mechanisms (Vainshtein 1972;Jain et al. 2013).The operation of these screening mechanisms is the most influential factor distinguishing different MG models, and their interaction and response with the Large Scale Structure and the cosmological environment play a pivotal role in governing the effects induced by MG theories.In this study, we limit our investigation to two widely recognized MG models that adhere to the screening mechanism.The first model, referred to as Hu-Sawicki f (R) gravity, incorporates additional scalar fields and their interaction with matter, allowing for nonlinear functions of the Ricci scalar (Hu & Sawicki 2007).The second model, known as the normal branch of the Dvali-Gabadadze-Porrati (nDGP) model, explores the possibility that gravity propagates in extra dimensions, distinguishing it from other conventional forces (Dvali et al. 2000).Both of these complex MG theories share a common characteristic: the emergence of a fifth force on cosmological scales, which arises due to the existence of extra degrees of freedom.In such theories that modify GR on large scales, the evolution of perturbations deviates from what is predicted by the standard model of cosmology.This modification introduces a fifth force that is expected to influence the formation of structures in the Universe (see, e.g., Carloni et al. 2008;Hu et al. 2016;Corasaniti et al. 2020;Mitchell et al. 2021).A notable feature of the extra degrees of freedom in the aforementioned MG models is their adaptability depending on the environment, which is known as their chameleon nature.Specifically, they are light in areas of low density but heavy in high-density regions of matter (Brax et al. 2013;Joyce et al. 2016;Burrage & Sakstein 2018).
In light of these discussions, it is plausible to anticipate that MG theories may result in alterations to the theoretical constructs of dark matter halos.Hence, the objective of this research is to calculate the merger rate of PBH-NS binaries within the context of Hu-Sawicki f (R) and nDGP gravitational models.In this regard, the structure of the study is organized as follows: Section 2 is dedicated to outlining the theoretical underpinnings of MG models.In Section 3, we explain the principles of dark matter halo models and key elements such as the halo density profile, the halo concentration parameter, and the halo mass function, all within the theoretical frameworks of MG models.Furthermore, in Section 4, we delineate the merger frequency of PBH-NS binaries in MG models and compare it with the analogous results from GR.In Section 5, we deliberate on the results.Lastly, in Section 6, we encapsulate the findings and draw conclusions.

MODIFIED GRAVITY MODELS
In this section, we focus on two fundamental models discussed in the field of MG literature: the Hu-Sawicki f (R) model and the nDGP model.These models embody the screening mechanisms known as the Chameleon and Vainshtein classe.

Hu-Sawicki f (R) Model
The Hu-Sawicki f (R) model involves introducing a non-linear modification function, denoted as f (R), to the conventional Einstein-Hilbert action (Hu & Sawicki 2007) where R is the Ricci scalar, κ is the Einstein gravitational constant, g is the metric determinant, and L m is the matter Lagrangian.By selecting f = −2Λ, it is worth noting that one can restore GR while also incorporating a cosmological constant.By employing a conformal transformation, Eq. ( 5) can be converted to a scalar-tensor theory involving the scalaron, denoted as f R ≡ df (R)/dR.This scalaron serves as the degree of freedom introduced by MG.Selecting an appropriate functional form for f (R) is crucial to ensure the occurrence of cosmic acceleration in the late-time Universe while still adhering to the limitations imposed by the solar system tests (Brax et al. 2008).A class of broken power-law models, referred to as the Hu-Sawicki model, is proposed in Hu & Sawicki (2007) that effectively addresses the aforementioned situation.Specifically, this model can be described as follows where m 2 = κρ m0 /3 represents the characteristic mass scale, and ρm0 denotes the background density of matter in the present-time Universe.Also, c 1 , c 2 , and n > 0 are dimensionless free parameters that need to be determined in a specific manner to accurately recover the expansion history and pass solar-system tests, using the chameleon mechanism.Note that the stability of the solution in high-density areas, i.e., where R ≫ m 2 , has to be maintained.In addition, cosmological experiments based on the f (R) model should align with experiments based on GR.To fulfill this requirement, one must have f RR = d 2 f /dR 2 > 0. Hence, one can expand the Hu-Sawicki model as follows lim Despite the absence of an actual cosmological constant in the Hu-Sawicki model, the model exhibits characteristics resembling a cosmological constant in both largescale and local experiments.Moreover, the finite value of c 1 /c 2 results in a constant curvature that remains unaffected by variations in matter density.Consequently, this choice allows for a class of models capable of accelerating the expansion of the Universe, which is similar to the behavior observed in the standard model of cosmology.Therefore, Eq. (3) can be reformulated as where R0 denotes the present-day background curvature, and f R0 ≡ f R ( R0 ) is the field strength1 .Cosmological and solar-system tests have limited the field strength f R0 (Desmond & Ferreira 2020).In this regard, various values of f R0 have been explored in the literature, spanning from 10 −4 to 10 −8 (Mirzatuny & Pierpaoli 2019).In our study, we concentrate on the Hu-Sawicki model with n = 1 and |f R0 | = 10 −4 , 10 −5 , and 10 −6 , denoted as f4, f5, and f6, respectively.

nDGP Model
The normal branch Dvali-Gabadadze-Porrati (nDGP) model of gravity is a modified gravity theory proposed in Dvali et al. (2000).The nDGP model describes the Universe as a four-dimensional brane embedded in a fivedimensional Minkowski space.The model assumes an action comprising two terms as where R 5 , g 5 , and κ 5 represent the Ricci scalar, metric determinant, and Einstein gravitational constant of the fifth dimension, respectively.Moreover, r c = (κ 5 /2κ) denotes the crossover distance, which signifies a specific scale below which GR transforms into a fourdimensional model.For scales bigger than r c , the second term in the action becomes dominant, resulting in expected deviations from GR.This model is composed of two branches: the "normal" branch (nDGP), provided here, and the selfaccelerating branch wich is called sDGP.We concentrate on the former since it is free of ghost instabilities (e.g., Lombriser et al. 2009).At larger scales, gravity becomes stronger, while at smaller scales, gravity behaves like GR due to Vainshstein screening.Therefore, one can investigate the nDGP model coupled with dark energy to match the desired ΛCDM expansion history.This path remains compelling because of prior investments in simulations.In this case, the only adjustable parameter to be constrained is n = H 0 r c (where H 0 stands for the Hubble constant), with values between 1 and 5 being extensively studied.Note that GR will be recovered if n → ∞, which corresponds to a large gradients of gravitational force in Vainshtein screening.The nDGP model has been the subject of various studies and simulations to understand its implications for structure formation and cosmology, and its properties were investigated through numerical simulations and comparisons with observational data (see, e.g., Hernández-Aguayo et al. 2021;Fiorini et al. 2022;Naidoo et al. 2023).In this work, we adopt three permissible values of n = 1, 2, and 5, and we denote the models associated with these values as nDGP(1), nDGP(2), and nDGP(5), respectively.

DARK MATTER HALO MODELS
In this section, we provide a brief overview of dark matter halo models derived from MG models.These models are used as initial inputs in calculating the merger rate of binary systems consisting of PBHs and NSs.

Halo Density Profile
Within the framework of cosmological perturbation theory, dark matter halos are identified as dynamical structures that exist in the nonlinear regime.These halos exhibit a density distribution that can be characterized by a function of radius, known as the halo density profile.In other words, a halo density profile describes the distribution of dark matter within a galactic halo and can be used to predict the dark matter distribution within the host halo.
Over the years, researchers have used various analytical methods and numerical simulations to develop halo density profiles that accurately represent observational data, particularly related to the rotation curves of galaxies (see, e.g., Einasto 1965;Jaffe 1983;de Zeeuw 1985;Hernquist 1990;Dehnen 1993;Navarro et al. 1996).The Navarro-Frenk-White (NFW) profile is one of the most commonly used density profiles for dark matter halos (Navarro et al. 1996).The NFW profile is a twoparameter functional form that describes the halo mass distribution as a function of distance from the halo center.Another density profile that is obtained from analytical models is the Einasto profile, which is a generalization of a power law distribution with functionality of distance and shape parameter (Einasto 1965).For both of the mentioned profiles, it is required that the logarithmic slope of the density distribution is −2 at the scale radius r s .

Halo Concentration Parameter
The concentration parameter is a dimensionless quantity that fundamentally represents the central density of galactic halos, determined by the ratio of the halo virial radius (r vir ) to its scale radius (r s ).The halo virial radius encompasses a volume where the average density is 200 to 500 times the critical density of the Universe.The concentration parameter holds significance in cosmology as it describes the shape of the halo density profile and the abundance of subhalos.Numerical simulations and analytical studies suggest that for precise predictions, the concentration parameter must dynamically change with halo mass and redshift (e.g., Prada et al. 2012;Dutton & Macciò 2014;Ludlow et al. 2016;Okoli & Afshordi 2016).This corresponds to the dynamics related to the merger history of dark matter halos and their evolutions.This correlation arises from the fact that smaller halos have already reached a state of virialization, leading to a higher concentration compared to host galactic halos.In this study, we use the concentration parameter defined in Okoli & Afshordi (2016) for calculations within the framework of GR.We also utilize the concentration parameters derived from Mitchell et al. (2019) and Mitchell et al. (2021) for computations related to the Hu-Sawicki f (R) model of gravity and the nDGP model of gravity, respectively.

Halo Mass Function
The halo mass function plays a vital role in characterizing halos according to their mass.It offers an extensive insight into the distribution of dark matter halos with respect to mass.The density contrast is denoted as δ(x) ≡ [ρ(x) − ρ]/ρ, where ρ(x) signifies the density of a region denser than the average at point x, and ρ represents the mean density of the background.Essentially, the halo mass function serves as a crucial instrument for categorizing cosmological structures based on their mass, particularly identifying structures that surpass a specific density contrast threshold, which triggers their collapse and subsequent halo formation.
The precise forecasting of the halo mass function serves as a crucial evaluation of cosmological models.In this regard, the Press-Schechter (PS) method provides an analytical framework for modeling the statistics of hierarchical structure formation originating from initial density fluctuations (Press & Schechter 1974).This method allows the formation of dark matter halos to be represented through stochastic processes and Gaussian density fields that are smoothed by a filter function.The PS formalism falls short in predicting halo abundances.For this purpose, several approaches have been carried out to overcome this issue2 .
The PS model was initially developed in Del Popolo & Gambera (1998), where it was found that the collapse threshold is not a constant value but rather, it is dependent on the mass.Consequently, the collapse threshold can be characterized as where α = 0.585, and β = 0.46.In addition, ν = (δ c /σ) 2 , where σ(M, z) is the linear root-mean-square fluctuation of overdensities on a comoving scale with mass M at redshift z.Using the excursion set formalism, one can demonstrate that Eq. ( 6) corresponds to the following barrier where α = 0.585, β = 0.46, α 2 = 0.4, β 2 = 0.02, α 3 = 0.45, and β 3 = 0.29.Also, Ω Λ ≃ 0.692 is the density parameter of cosmological constant.
In the framework of the excursion set theory, the mass function is defined as the comoving number density of halos within a particular mass interval (Bond et al. 1991): where f (ν) is known as multiplicity function that denotes the distribution of the first crossing.
Also, the constants are assigned specific values, namely A 2 = 0.93702 and a = 0.707.By substituting Eq. ( 10) into Eq.( 9), a precise semianalytical mass function can be derived, referred to as the mass function of GR henceforth.
Nonetheless, in MG theories, it is necessary to rescale the mass function.This is because, in those models that introduce an additional scalar degree of freedom, Birkhoff's theorem is violated.This results in the collapse process being contingent on the environment, a phenomenon known as the chameleon screening mechanism.In this mechanism, dense regions partially mitigate the amplified gravitational forces.Consequently, the threshold of surplus locations depends on the halo mass, redshift, and the particular configuration of the model parameters, which alters the gradients of the scalar field.
In Gupta et al. (2022), a rescaled halo mass function for the Hu-Sawicki f (R) gravity is presented based on N -body simulations: where X ≡ ln(σ −1 ), which naturally depends on M vir and z.Also, based on varoius models of Hu-Sawichi f (R) gravity, the best values of A, B, and C are provided in Table 1.Furthermore, in Mitchell et al. (2021), a suitable form for the halo mass function in nDGP models is proposed, employing dark matter only simulations: where 14) It is important to highlight that the aforementioned parameters represent the best-fit values, derived via the method of unweighted least squares.

PBH-NS MERGER RATE
Consider a scenario where a PBH with mass m 1 and a NS with mass m 2 are situated within an isolated dark matter halo.Assume they abruptly cross paths on a hyperbolic trajectory, and their relative velocity at a Given this premise, two-body scattering suggests that significant gravitational radiation is emitted at the nearest physical separation, referred to as the periastron.If the radiated gravitational energy surpasses the system's kinetic energy, the compact objects will become gravitationally bound, resulting in a binary formation.The upper limit of the periastron is dictated by this condition (Peters 1964): where G represents the gravitational constant and c denotes the speed of light.Furthermore, within the framework of the Newtonian approximation, the impact parameter can be defined as (Sasaki et al. 2018): Furthermore, when the strong limits of gravitational focusing are set, that is, r p ≪ b, the tidal forces exerted by the nearby compact objects on the binary can be disregarded.Consequently, the cross-section for the binary formation is given by (Quinlan & Shapiro 1989;Mouri & Taniguchi 2002) Therefore, the rate at which binaries form within each galactic halo conforms to the following equation (Bird et al. 2016): where 0 < f PBH ≤ 1 represents the fraction of PBHs determining their contribution to dark matter, ρ h denotes the halo density profile, which can be modeled using the NFW or Einasto profiles.The angle bracket is an average over the PBH relative velocity distribution within the galactic halo.Additionally, ρ NS represents the density profile of NSs, characterized by the following spherically-symmetric form: The parameters ρ * NS and r * NS represent the characteristic density and radius of NSs, respectively.Eq. ( 19) illustrates that to determine the distribution of NSs, two quantities must be specified.Firstly, the defining radius of a NS, r * NS ≃ (0.01-0.1)r s , is considered (Sasaki et al. 2022).Secondly, the defining density of NSs, ρ * NS , is established by normalizing the NS distribution to their projected population within a specific galaxy.In order to achieve this, we employ the time-constant form of the initial Salpeter stellar mass-function, denoted as χ(m * ) ∼ m −2.35 * .Moreover, we posit that all stars within the mass range of M NS ≃ (8-20) M ⊙ undergo supernova explosions, leading to the formation of NSs.As a result, the number count of NSs within a galaxy that has a stellar mass M * is determined as follows: where χ(m * )m * is normalized to unity.It is crucial to determine the galactic stellar mass by defining the correlation between stellar mass and halo mass, represented as M * (M vir ).To achieve this, we employ the stellar mass-halo mass relationship outlined in Behroozi et al. (2013), under the assumption that the highest concentration of NSs is located in the central region of the galactic halo.
The calculation of the merger rate of PBH-NS binaries also depends on dark matter particle velocity dispersion within galactic halos.Recent studies suggest that MG model velocity dispersion profiles bear a close resemblance to those in the conventional cosmological model (He et al. 2015).Consequently, it is justifiable to apply the dark matter particle velocity dispersion relation, as derived in Prada et al. (2012) to our analysis.Furthermore, we set the probability distribution func-tion of relative velocities among PBHs in the galactic halos to be regulated by the Maxwell-Boltzmann distribution, which is truncated at the virial velocity (Fakhry et al. 2021(Fakhry et al. , 2022a)).
To calculate the total merger rate of PBH-NS binaries, one can multiply the halo mass function, dn(M )/dM , by the rate of binary formation within each halo, Γ(M ), and integrate over a minimum mass of halos: The above formula indicates that the ultimate result is minimally influenced by the upper limit.This can be attributed to the fact that the halo mass function incorporates a term that decreases exponentially, rendering the merger rate to decrease as the halo mass increases.Conversely, the lower limit holds significant importance.This concurs with the dynamics of hierarchical structure formation, given that the smallest halos are anticipated to exhibit greater density than the host halos.Nonetheless, dynamical relaxation processes lead to the evaporation of the smallest halos by ejecting objects.Furthermore, merger and accretion processes lose their effectiveness during the dark energy dominated era.Consequently, the smallest halos could have entirely evaporated in the past 3 Gyr, leading to the potential disregard of signals from some halo mergers (Bird et al. 2016).Hence, the evaporation timeframe of dark matter halos guarantees their presence in the presenttime Universe.In our prior studies, we have established that the evaporation period of dark matter halos with a standard mass of 400 M ⊙ , hosting PBHs with a mass of 30 M ⊙ , is approximately 3 Gyr (Fakhry et al. 2021(Fakhry et al. , 2022a)).As a result, one can consider dark matter halos with the aforementioned typical mass as the minimum mass in our analysis.However, dark matter halos hosting PBHs lighter than 30 M ⊙ can possess slightly lower masses than 400 M ⊙ to persist today, while those hosting PBHs heavier than 30 M ⊙ can have slightly higher masses than 400 M ⊙ to fulfill the conditions.

RESULTS AND DISCUSSIONS
Fig. 1 illustrates the merger rate of PBH-NS binaries within each halo as a function of halo mass for Hu-Sawicki f (R) and nDGP gravities with different parameter values, compared with GR predictions.We have employed the NFW density profile and assumed that PBHs have a mass of 5 M ⊙ and NSs have a mass of 1.4 M ⊙ for this calculation.Moreover, we have considered that PBHs account for all of the dark matter, that is, f PBH = 1.In this figure, calculations have been performed for the permissible boundaries of the characteristic radius of neutron stars, i.e., r * NS = 0.01, and 0.1 r s .As expected, the respective outcomes for any other potential values lie between these two limits.The left subfigure indicates that the merger rate of PBH-NS binaries in each halo is greater for all Hu-Sawicki f (R) models than for GR.This is attributed to the nonlinear effects of f (R) gravity on density fluctuations and halo formation and evolution, which are incorporated in all stages of the calculations.Another notable finding is that the impact of the field strength on density fluctuations diminishes progressively from f4 to f6.This is because the merger rate of PBH-NS binaries is dependent not only on density fluctuations but also on the field strength f R0 .A similar result can be found by looking at the right subfigure.Accordingly, the merger rate of PBH-NS binaries in each halo in nDGP models is much higher than the results obtained for GR.It is also clear that the enhancement of the merger rate is highest for the nDGP(1) model, moderate for the nDGP(2) model and lowest for the nDGP(5) model.
In Fig. 2, we have shown the merger rates of PBH-NS binaries per unit time and volume in Hu-Sawicki f (R) and nDGP models, and compared them with the GR results using NFW density profile.These calculations are conducted for the present-time Universe.We have applied Eqs. ( 10), ( 12), and (13) for the halo mass functions of GR, Hu-Sawicki f (R), and nDGP models, respectively.We have also utilized the relations derived in (Ludlow et al. 2016;Gupta et al. 2022;Mitchell et al. 2021) as the concentration parameters correspondingly.The merger rate of PBH-NS binaries decreases as the minimum mass of dark matter halos increases.This means that the binary merger rate is more affected by smaller halos than larger ones.These results support the dynamics of hierarchical structures, which imply that the mass of halos and the density of dark matter particles are inversely related.Therefore, the subhalos have more dark matter particles and higher concentration than the host halo.It is noteworthy that the merger rate of PBH-NS binaries can be precisely determined by integrating across the area under the curves.Our analysis reveals that for any Hu-Sawicki f (R) model, the merger rate of PBH-NS binaries is higher than that of GR.Accordingly, the merger rate of PBH-NS binaries clearly depends on the field strength f R0 .As the field strength becomes smaller, the merger rate approaches the GR result.Based on this argument, one can utilize the value of the field strength and GW data to compare theoretical predictions with observations and estimate the number of PBHs.This approach holds potential as a new technique.
A similar pattern for the nDGP models can also be obtained.The difference between the merger rate of PBH-NS binaries in nDGP gravity and GR is more pronounced for nDGP(1) than for nDGP(2) and nDGP(5) respectively.This means that the nDGP models have different effects on the merger rate of PBH-NS binaries depending on the parameter n = H 0 r c , which controls the strength of the extra dimension.We have also quan-tified the amplification in the merger rate of PBH-NS binaries for both gravitational models.The two models yield different reinforcement processes.Nevertheless, the growth of the merger rate for both models is inversely related to the minimum mass of halos, leading to the most substantial growth occurring in low-mass halos.Furthermore, it is apparent that the enhancement in the merger rate of PBH-NS binaries is slightly greater in the nDGP models compared to the results obtained from Hu-Sawicki f (R) models.
The specific goal of this analysis is to illustrate the merger rate of PBH-NS binaries throughout the evolution of the Universe, facilitating comparison with the merger events detected by the LVK observatories.This effort is grounded in gravitational models and aims to offer theoretical forecasts, elucidating the prospective trajectory of GW phenomena and the progression of these detection instruments..The current sensitivity of GW detectors allows the detection of merger events within a comoving volume of 50 Gpc 3 , corresponding to a redshift of roughly 0.75 z.This presents a challenging issue in calculating the redshift evolution of the merger rate of PBH-NS binaries.The sensitivity to changes in redshift is emphasized, particularly in Eq. ( 21 the halo mass function and the concentration parameter.In Fig 3, we have depicted the redshift evolution of the merger rate of PBH-NS binaries for f (R) and nDGP models, contrasting the results with those for GR while taking into account NFW and Einasto density profiles.The direct relation between the merger rate of PBH-NS binaries and redshift is evident, which can be attributed to the influence of hierarchical dynamics and halo merger structure.The findings suggest that the redshift evolution of the merger rate PBH-NS binaries in the mentioned gravitational models is higher than that obtained from the GR framework, potentially leading to an amplified binary merger rate while going back in time.Furthermore, one has to highlight the effect of model parameters, like f R0 and n = H 0 r c , on the increase of the merger rate of PBH-NS binaries during the late-time Universe.Furthermore, the merger rates derived from all models, considering the Einasto density profile, display an increase of approximately 60% compared to those obtained when considering the NFW density profile.
Fig. 4 presents the merger rate of PBH-NS binaries for both Hu-Sawicki f (R) and nDGP models of gravity as a function of PBH mass and the two characteristic radii of NS, compared with GR.We have only provided the results for the NFW density profile, knowing the proportionality of the results between the two profiles.Also given the linear scaling of the merger rates with PBH abundance, we have only shown the results for f PBH = 1.In all the models considered, the merger rate of PBH-NS binaries has an inverse relationship with the mass of PBHs.This is because the number density of PBHs changes inversely with their mass within a certain volume.Smaller PBHs are more likely to participate in PBH-NS binary formations than larger ones.This is be-cause the formation of PBH-NS binaries is influenced by factors such as the mass range of PBHs and the density contrast boost factor (Tsai et al. 2021).
In our previous discussions, we assumed PBHs can provide the maximum contribution to dark matter.Nevertheless, it is interesting to determine the merger rate of PBH-NS binaries, considering various assumptions about PBH fractions and masses.To this end, in fig.5, we have depicted the upper limit of the cumulative merger rate of PBH-NS binaries as a function of the mass spectrum of stellar mass PBHs for both Hu-Sawicki f (R) and nDGP models.This calculation incorporates both NFW and Einasto density profiles.It is once again clear that the merger rate of PBH-NS binaries is inversely related to the mass of PBHs.
Moreover, the findings indicate that the Hu-Sawicki f (R) and nDGP theoretical models enforce more stringent constraints on the fraction of PBHs, particularly when smaller masses of PBHs contribute to binary systems.The figure also suggests that the value of the field strength f R0 plays a significant role in constraining the fraction of PBHs.This reveals a direct relation between the field strength and the severity of the constraints on a fraction of PBHs.In this context, the most rigorous constraints can be derived from the f (R) gravity with the field strength f4, f5, and f6 respectively.Similar results are inferred from nDGP models, wherein the most stringent constraints are imposed on the abundance of PBHs from nDGP(1), nDGP(2), and nDGP(5) models, respectively.
Additionally, it can be deduced that the most precise constraint in this study is derived from the f (R) models, specifically when the Einasto density profile and M PBH ≃ 5 M ⊙ are considered.A similar outcome is obtained in the nDGP model with the same PBH mass and density profile, particularly when the nDGP(1) model is taken into account.Despite the presence of theoretical uncertainties, if the proportion of PBHs comprises roughly over 10% of the total dark matter content, the findings align with the LVK sensitivity band.
Nonetheless, it is important to highlight that the constraints derived from the merger rate of PBH-NS binaries in these models signify the upper limits allowed by GW detectors.This is attributed to the fact that the BH component in the mergers, which originate from astrophysical sources, can also be detected by the LVK detectors.Our analysis suggests that, to annually observe a minimum of O(10), O(1), and O(0.1) merger events of PBH-NS binaries with a mass of (M PBH = 5 M ⊙ , M NS = 1.4 M ⊙ ) in the comoving volume of 50 Gpc 3 , the anticipated upper bounds on the abundance of PBHs in both gravitational models are estimated to be f PBH ∼ O(10 −1 ), O(10 −2 ), and O(10 −3 ), respectively.
However, there exist robust constraints from non-GW data on the abundance of PBHs within the mass range of 5-50 M ⊙ .In Figure 6, we have shown a comparative analysis between the outcomes of f (R) and nDGP models concerning the upper bounds on the abundance of annual detections of PBH-NS merger events are expected in the comoving volume of 50 Gpc 3 , with MPBH = 5M⊙, MNS = 1.4M⊙.Non-GW constraints on PBH dark matter, including gravitational microlensing constraints from quasars (ML Quasar), constraints from the disruption of wide binaries (WB), constraints from the disruption of ultra-faint dwarfs (UFD), X-ray constraints related to accreting PBHs (X-ray I), and the corresponding one in the Milky Way (X-ray II), and constraints from modification of the cosmic microwave background (CMB) spectrum due to accreting PBHs, including particle dark matter accretion, as well as the FIRAS data are provided as benchmarks for this analysis.
PBHs and the non-GW data related to the examined mass range.The non-GW data include constraints from gravitational microlensing of quasars (ML Quasar) and OGLE-II+III, disruption of wide binaries (WBs), disruption of ultrafaint dwarf galaxies (UFDs), X-ray constraints associated with accreting PBHs (X-ray I) and its Milky Way counterpart (X-ray II), alterations of the cosmic microwave background (CMB) spectrum due to accreting PBHs (including particle dark matter accretion), and FIRAS data.The merger rates of PBH-NS binaries within the context of nDGP models exhibit a marginal increase compared to those predicted by f (R) models.Consequently, when comparing the two theoretical frameworks, it becomes apparent that the nDGP models impose slightly more rigorous constraints on the abundance of PBHs than those inferred from f (R) models.
The expected upper bounds on the PBH fraction from both gravitational models are potentially detectable about (i): O(10), (ii): O(1), or (iii): O(0.1) merger events per year of (M PBH = 5M ⊙ , M NS = 1.4M ⊙ ) in the comoving volume detectable by LVK detectors.Upon comparison, it is clear that the relevant constraints for case (i) in both gravitational models are less stringent than those derived from non-GW data.That is why we have not provided their results in the figure.However, other cases for the upper limit on PBH abundance may remain more stringent than the constraints derived from non-GW data.In this regard, case (ii) imposes less stringent constraints than those obtained from X-ray I, X-ray II, and CMB observations (including particle dark matter accretion) for 20 M ⊙ ⩽ M PBH ⩽ 100 M ⊙ , and is more stringent than those constraints obtained from non-GW data for M PBH < 20 M ⊙ .Case (iii) is also less stringent than the constraints derived from CMB observations (including particle dark matter accretion) for 50 M ⊙ ⩽ M PBH ⩽ 100 M ⊙ , and is more stringent than those constraints obtained from non-GW data for M PBH < 50 M ⊙ .This suggests that within the parameter space not excluded by non-GW data, it may be feasible to rely on upper limits concerning the expected abundances of PBHs in both gravitational models.

CONCLUSIONS
In this study, we have elucidated the merger rate of PBH-NS binaries within the framework of two extensively-studied MG models: Hu-Sawicki f (R) gravity and the nDGP gravity.By employing the key parameters of dark matter halos within MG models, including the halo density profile, halo concentration parameter, and halo mass function, we have demonstrated an enhancement in the PBH-NS merger rate compared to GR.This enhancement stems from an extra scalar degree of freedom in these MG theories, which generates a fifth force that enhances gravitational interactions on cosmological scales.
The results also demonstrate that the merger rate increases as the field strengths, characterized by the MG parameters in the two sets of models, increase.Specifically, this includes the field strength, denoted as f R0 , for f (R) gravity, and n = H 0 r c for nDGP gravity.Moreover, lower mass PBHs display higher merger rates due to their increased number density for a given mass fraction.Integration the merger rate across the halo weighted by the halo mass function reveals that the total merger rate is significantly higher in MG models compared to GR.Furthermore, an analysis of the redshift evolution suggests that this deviation from GR is more pronounced at higher redshifts.
Upper limits on the abundance of PBHs have been set based on sensitivity bands of LVK detectors.In this context, the predicted merger rates of PBH-NS binaries, as derived from MG theories, impose a constraint on the fraction of PBHs, denoted as f PBH ≳ 0.1.This im-plies that the theoretical predictions from MG theories could be useful in determining the abundance of PBHs in the Universe.However, additional studies in multimessenger astronomy could provide more information on these constraints and help improve our understanding of the complexities involved.
Moreover, compared with the existing non-GW constraints on the abundance of PBHs, our analysis within the MG models could potentially impose more rigorous constraints on the fraction of PBHs with masses below 20 M ⊙ to be f PBH ≲ 0.01, and for masses below 50 M ⊙ to be f PBH ≲ 0.001.Additionally, slightly more stringent constraints can be discerned for nDGP models in comparison to f (R) models.These results highlight how MG theories can produce accurate predictions of the merger rate of PBH-NS binaries and the abundance of PBHs.Future observations of GW events could confirm these predictions and elucidate gravitational phenomena beyond GR.
In conclusion, we have demonstrated the impact of large-scale gravitational variations on the merger rate of PBH-NS binaries and constraints on the abundance of PBHs.Despite many uncertainties, the merger rates predicted by MG models emphasize the potential of applying this effect to GW observations to assess the validity of these theories against GR.This implies that upcoming research in this area should leverage the supplementary effects of MG theories to expand the theoretical foundations.However, it should be considered that effective impacts may extend beyond the factors contemplated in this study, and the degeneracy in the parameter space should be considered as a potential uncertainty.

Figure 1 .
Figure 1.The merger rate of PBH-NS binaries within each halo as a function of halo mass in the present-time Universe, considering Hu-Sawicki f (R) gravity (left) and nDGP gravity (right), for various parameter values of the model.GR predictions are over-plotted for comparison.Solid and dashed lines correspond to r * NS = 0.01rs, and 0.1rs, respectively.The halo density is assumed to follow the NFW profile.We have assumed MPBH = 5 M⊙, MNS = 1.4 M⊙, and fPBH = 1.

Figure 2 .
Figure 2. Top: The total merger rate of PBH-NS binaries as a function of halo mass in the present-time Universe, for Hu-Sawicki f (R) (left) and nDGP gravities (right), and with the NFW density profile for the halos, compared to GR predictions.Solid and dashed lines correspond to r * NS = 0.01rs, and 0.1rs, respectively.Bottom: The relative enhancement of the merger rate from MG predictions to those from GR as a function of halo mass.We have assumed MPBH = 5 M⊙, MNS = 1.4 M⊙, and fPBH = 1.

Figure 3 .
Figure 3. Upper limit of the merger event rate of PBH-NS binaries as a function of redshift, considering Hu-Sawicki f (R) gravity (left) and nDGP gravity (right).These results are compared with those obtained in GR.Both NFW and Einasto density profiles are employed.We have assumed MPBH = 5 M⊙, MNS = 1.4 M⊙, and fPBH = 1.

Figure 4 .
Figure 4. Total merger rate of PBH-NS binaries as a function of different PBH masses, considering Hu-Sawicki f (R) gravity (left) and nDGP gravity (right).These results are compared with those obtained in GR.Solid and dashed lines correspond to r * NS = 0.01rs, and 0.1rs, respectively.The NFW density profile is employed.We have assumed MNS = 1.4 M⊙, and fPBH = 1.

Figure 5 .
Figure5.Upper limit on the total merger event rate of PBH-NS binaries with respect to the fraction and mass of PBHs, considering Hu-Sawicki f (R) gravity (left) and nDGP gravity (right).The shaded gray bands represent the total merger rate BH-NS binaries, as estimated by the LIGO-Virgo-KAGRA detectors during their latest observing run, within the range of (7.8 − 140) Gpc −3 yr −1 .The shaded red and blue bands show the results for the NFW and Einasto density profiles, respectively.We have assumed MNS = 1.4 M⊙.

Figure 6 .
Figure6.The expected upper bounds on PBH fraction as a function their mass, for Hu-Sawicki f (R) (black lines) and nDGP (brown lines) gravities.The dashed and dotdashed lines represent scenarios for which O(1) and O(0.1) annual detections of PBH-NS merger events are expected in the comoving volume of 50 Gpc 3 , with MPBH = 5M⊙, MNS = 1.4M⊙.Non-GW constraints on PBH dark matter, including gravitational microlensing constraints from quasars (ML Quasar), constraints from the disruption of wide binaries (WB), constraints from the disruption of ultra-faint dwarfs (UFD), X-ray constraints related to accreting PBHs (X-ray I), and the corresponding one in the Milky Way (X-ray II), and constraints from modification of the cosmic microwave background (CMB) spectrum due to accreting PBHs, including particle dark matter accretion, as well as the FIRAS data are provided as benchmarks for this analysis.