  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
(nDGP) 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, PBH mass and fraction. Assuming PBH mass range of $5-50
M_{\odot}$, 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 with the existing non-GW constraints, which can potentially
impose even more stringent constraints.
