Neutron-rich materials released from mergers of binary neutron stars (BNS) and black hole-neutron star (BH-NS) binaries undergo rapid neutron capture process (r-process), which has long been proposed to account for most of the elements heavier than iron in the universe (Lattimer & Schramm 1974 Symbalisty & Schramm 1982) . Li & Paczyński (1998) first predicted that the radioactive decays of these r-process nuclei can provide a heating source for powering a new type of transient days to week following a BNS or BH-NS merger. The detailed research by Metzger et al. (2010) suggested that the luminosity of such an radioactivity-powered transient is around several times 10 41 erg s −1 , which is a factor of ∼ 10 3 higher than that of a typical nova, so that such a transient is called a "kilonova". The characteristics of the lightcurves, peak luminosities, timescales, and spectra of the kilonova emission have been widely investigated theoretically (Kulkarni 2005; Rosswog 2007; Roberts et al. 2011; Kasen et al. 2013 Kasen et al. , 2015 Kasen et al. , 2017 Tanaka & Hotokezaka 2013; Yu et al. 2013; Grossman et al. 2014; Metzger & Piro 2014; Metzger & Fernández 2014; Perego et al. 2014; Wanajo et al. 2014; Just et al. 2015; Martin et al. 2015; Fernández & Metzger 2016; Metzger 2017; Kawaguchi et al. 2020) .
In the past decade, numerical simulations revealed that the physical properties of the kilonova ejecta could have a large diversity. During the dynamical timescale of the onset of BNS and BH-NS mergers, some highly neutron-rich matter would be directly ejected by angular momentum via tidal force (e.g., Hotokezaka et al. 2013; Bauswein et al. 2013; Foucart et al. 2013; Kyutoku et al. 2013 Kyutoku et al. , 2015 . The lanthanide-rich "red" tidal tails would expand outwards in the equatorial plane due to angular momentum transport by hydrodynamic processes. There also exist some of the less neutron-rich matter unbound from a BNS or a BH-NS merger. For BNS mergers, a lanthanide-free "blue" component mainly resides in the polar regions, due to heating by the shocks at the contact interface between two merging NSs that operate on the dynamical time (e.g., Oechslin & Janka 2006; Wanajo et al. 2014; Radice et al. 2016; Sekiguchi et al. 2016) or by neutrino irradiation from the remnant NS (e.g., Metzger & Fernández 2014; Perego et al. 2014; Martin et al. 2015; Kasen et al. 2015) . The amount of blue-component ejecta driven by neutrino irradiation is significantly affected by the survival time of the merger remnant NS Fujibayashi et al. 2017; Gill et al. 2019 ). In addition, the initial high-temperature remnant disk around the remnant NS or BH can release thermal neutrinos to form a handful of lanthanide-free ejecta toward the polar direction (e.g., Fernández et al. 2015; Just et al. 2015; Fujibayashi et al. 2020) . It is predicted that the observed kilonova emission from a BNS and BH-NS merger would be a combination of a "blue" and a "red" component. While both components could be simultaneously visible for the line of sight near to the polar direction, the blue emission is possible to be blocked by the red component for observers with an equatorial line of sight. Except for these two components, it could form an intermediate-opacity equatorial-dominated ejecta by viscous heating and angular momentum transport of the remnant disk (e.g., Fernández & Metzger 2013; Just et al. 2015; Siegel & Metzger 2017) . Furthermore, the central compact remnant after a merger could provide additional energy, e.g., the spin energy of a magnetar Metzger & Piro 2014) or the fall-back accretion energy for a BH (Rosswog 2007; Kyutoku et al. 2015; Ma et al. 2018) , to enhance the brightness of the kilonova. Overall, the comprehensive characteristics of a kilonova depend on many factors, including the central remnant after the merger, physical properties of the ejected material, and potential additional energy injection. The asymmetric geometry of the ejecta also provides viewing-angle-dependent properties of kilonova lightcurves and polarization properties Bulla et al. 2019 Bulla et al. , 2020 Kawaguchi et al. 2020; Korobkin et al. 2020; Zhu et al. 2020; Darbha & Kasen 2020; Kóbori et al. 2020) .
Despite the clear predictions of kilonova emission from BNS and BH-NS mergers, its low luminosity and fast evolution compared with supernova emission makes difficult to detect kilonovae. Traditional optical time-domain survey projects have been designed to discover ordinary supernovae or tidal disruption events whose rise and decay timescales in brightness are long, e.g. weeks or even months. Without target-of-opportunity observations, only some kilonovae candidates could be recorded in the such time-domain survey projects. Recently, thanks to improved cadence of supernova searches and the operation of some new large survey projects, some fast-evolving transients, which had lightcurves lasting days to weeks and brightness in the range of theoretically predicted kilonova brightness, have been discovered (e.g., Drout et al. 2014; Prentice et al. 2018; Rest et al. 2018; McBrien et al. 2020) . However, none of these events have been firmly identified as kilonovae.
If short gamma-ray bursts (sGRBs) also originate from BNS or BH-NS mergers (Paczynski 1991; Paczynski & Rhoads 1993; Eichler et al. 1989; Narayan et al. 1992) , then another plausible way to constrain kilonova models is via follow-up observations of nearby sGRBs on timescales of hours to a week. In the past few years, there have been several kilonova candidates detected in superposition with the decaying sGRB afterglows (Berger et al. 2013; Tanvir et al. 2013; Fan et al. 2013; Gao et al. 2015 Gao et al. , 2017b Jin et al. 2015 Jin et al. , 2016 Jin et al. , 2020 Yang et al. 2015; Gompertz et al. 2018; Ascenzi et al. 2019; Rossi et al. 2020 ). However, due to the scarce and ambiguous observational data, these kilonova candidates also could not be firmly confirmed.
There are some limitations of searching kilonovae using time-domain survey observations or sGRBs follow-up observations. For time-domain survey observations, it is easy to miss the optimal observational time window for searching kilonovae due to their faint nature and short emission timescale. The lack of spectral observations or color evolution for potential kilonovae candidates make them hard to identify. For sGRBs follow-up observations, on the other hand, gamma-ray detectors usually have poor spatial positioning capabilities (e.g., Connaughton et al. 2015) . Current gamma-ray detectors also do not achieve an all-sky coverage. The beaming effect of sGRBs (Rhoads 1999; Sari et al. 1999 ) makes only a fraction of events being triggered. Also the putative kilonova emission is likely outshone by the luminous on-axis jet afterglow. In view that BNS and BH-NS mergers are targeted gravitational wave (GW) events for ground-based interferometer detectors, the strategy of optical follow-up observations for GW triggers to search kilonovae has been studied in detail (e.g., Metzger & Berger 2012; Cowperthwaite & Berger 2015; Gehrels et al. 2016; Cowperthwaite et al. 2019) . Compared with sGRBs, GW radiation is essentially isotropic from the source (Creighton & Anderson 2011) . Also GW detectors such as the Advanced Laser Interferometer Gravitational Wave Observatory (LIGO) have all-sky coverage and can detect sources out to distances of ∼ 300 and ∼ 650 Mpc for BNS and BH-NS mergers, respectively (Cutler & Thorne 2002) . GW sources will be localized to ∼ 10 deg 2 by the network of the the Advanced LIGO, Advanced Virgo and KAGRA GW detectors (Abbott et al. 2018a) . Therefore, an optimal searching strategy for kilonovae would be to take advantage of target-of-opportunity follow-up observations of GW triggers.
On 2017 August 17, a BNS merger event (GW170817) for the first time was observed in GWs by the Advanced LIGO and the Advanced Virgo detectors (Abbott et al. 2017a) . Independent gamma-ray observations discovered a sGRB (GRB 170817A; Abbott et al. 2017b; Goldstein et al. 2017; Savchenko et al. 2017; Zhang et al. 2018) , and follow-up observations discovered an ultraviolet-optical-infrared transient (AT 2017gfo; Abbott et al. 2017c; Andreoni et al. 2017; Arcavi et al. 2017; Chornock et al. 2017; Coulter et al. 2017; Covino et al. 2017; Cowperthwaite et al. 2017; Díaz et al. 2017; Drout et al. 2017; Evans et al. 2017; Hu et al. 2017; Kasliwal et al. 2017; Kilpatrick et al. 2017; Lipunov et al. 2017; McCully et al. 2017; Nicholl et al. 2017; Pian et al. 2017; Shappee et al. 2017; Smartt et al. 2017; Soares-Santos et al. 2017; Tanvir et al. 2017; Troja et al. 2017; Utsumi et al. 2017; Valenti et al. 2017 ) and a broadband afterglow from radio to X-rays consistent with off-axis jet emission (Alexander et al. 2017; Haggard et al. 2017; Hallinan et al. 2017; Margutti et al. 2017; Troja et al. 2017 Troja et al. , 2020 D'Avanzo et al. 2018; Dobie et al. 2018; Lazzati et al. 2018; Lyman et al. 2018; Piro et al. 2019; Ghirlanda et al. 2019) . The ultraviolet-optical-infrared transient was well consistent with the theoretical predictions of kilonova emission (Cowperthwaite et al. 2017; Perego et al. 2017; Tanaka et al. 2017 Tanaka et al. , 2018 Villar et al. 2017; Gao et al. 2017a; Ai et al. 2018; Kawaguchi et al. 2018; Li et al. 2018; Wanajo 2018; Yu et al. 2018; Ren et al. 2019; Wu et al. 2019) . The multi-messenger observations of this BNS merger officially opened the era of GW-led astronomy, and provided a smoking-gun evidence for the long-hypothesized origin of sGRBs and kilonovae.
Many open questions regarding BNS and BH-NS mergers, such as event rates, forming channels, NS equation of state (EoS), system parameters distributions, kilonova luminosity function, and diversity of kilonova, can be only answered based on statistical analyses of a large sample of GW and EM detections. Searching for more GW and EM data are conducive to the further understanding of the nature of BNS and BH-NS mergers. Over the course of the third observing run (O3) of the LIGO Scientific Collaboration and Virgo Collaboration (LVC), several compact binary merger candidates have been announced. Even though many efforts of follow-up observations for GW triggers during O3 Antier et al. 2020a,b; Coughlin et al. 2020a,b; Gompertz et al. 2020; Kasliwal et al. 2020; Page et al. 2020 ) have been carried out, no confirmed EM counterpart, especially kilonova, has been reported since the multi-messenger observations of GW170817 1 . The incomplete coverage of the error boxes of the GW events could partially be the reason of non-detection, but it is also likely that previous and present follow-up observations are not deep enough to search for EM counterparts (Sagués Carracedo et al. 2020; Coughlin et al. 2020a) . Figure 1 shows the timeline of GW observations (above the line), i.e. the starting times of the past LIGO and LVC observing runs and future networks of 2nd, 2.5th, and 3rd generation detector (see Section 4.1.2 for details), and EM observations (below the line), i.e. the starting times of large survey telescopes including the Panoramic Survey Telescope and Rapid Response System (PS1, Kaiser et al. 2010) , the Zwicky Transient Facility (ZTF, Graham et al. 2019) , the Multi-channel and their effects on the EM counterparts (Section 2.3). We note that the more massive BH-NS systems have stronger GW radiation intensities and also would merge sooner than less massive BH-NS systems. Therefore the redshift and system parameter distributions for the observed BH-NS merger events should in principle a result of the convolution of both system parameters and the intrinsic redshift distribution. In this paper, we ignore the possible redshift evolution of BH-NS system parameters to allow us to separately discuss BH-NS systems parameter distributions and the redshift distibution f (z). We describe our concrete models in the following.
The number density per unit time for BH-NS mergers at a given redshift z can be estimated as
whereρ 0 is the local BH-NS merger rate, and f (z) is the dimensionless redshift distribution factor. The comoving volume element dV (z)/dz in Equation (1) is
where c is the speed of light, and D L is the luminosity distance which is expressed as
There are many previous studies for the local BH-NS merger rateρ 0 . For recent binary population synthesis simulations, by considering the constraints of observed compact star binary samples and supernova rates, O'Shaughnessy et al. (2008) and Abadie et al. (2010) predicted that the BH-NS merger event rate density isρ 0 = 34.8 +1125 −34.2 Gpc −3 yr −1 in the local universe. By sampling from compact binary merger catalogs derived from N -body and populationsynthesis simulations, Santoliquido et al. (2020) (see also Rastello et al. 2020) presented the local event rate density asρ 0 = 38 +32 −24 Gpc −3 yr −1 for BH-NS mergers in young star clusters, andρ 0 = 45 +45 −32 Gpc −3 yr −1 for BH-NS mergers in the field. Observationally, the lack of BH-NS merger events detected during the first and second observing runs by LIGO and LVC gave the upper limit of the local BH-NS event rate asρ 0 < 610 Gpc −3 yr −1 . Antier et al. (2020b) systematically analyzed the GW candidates discovered during third run by LVC and estimated that the event rates of BH-NS and MassGap mergers are ∼ 27 Gpc −3 yr −1 and ∼ 3.3 Gpc −3 yr −1 , respectively. Based on these theoretical simulation results and observational constraints, we set the local event rate density of BH-NS mergers aṡ ρ 0 ≈ 35 Gpc −3 yr −1 hereafter.
Compact binary mergers occur with a delay time scale with respect to star formation history. The delay time distribution models include the Gaussian delay model (Virgili et al. 2011) , log-normal delay model (Wanderman & Piran 2015) , and power-law delay model (Virgili et al. 2011) . The observations of sGRB prefer to support either a Gaussian delay model (Virgili et al. 2011) or a log-normal model (Wanderman & Piran 2015) . Hereafter, we adopt log-normal delay model as our merger delay model. For the best fitting parameters consistent with sGRB data, we derive f (z) (see Appendix A) in Equation (1) 
where η = −5.51. By a method of Monte Carlo simulation based on Equation (1), one can generate a group of BH-NS events following the f (z) redshift distribution for the log-normal delay model.
In this section, we briefly introduce the model of kilonova emission from BH-NS mergers which has dependence on the latitudinal viewing angle θ view and longitudinal viewing angle ϕ view . Three different ejecta components, i.e. tidal dynamical ejecta, neutrino-driven wind ejecta, and viscosity-driven wind ejecta, are considered. More detailed descriptions of such a model as well as the related phenomenology have been presented in Paper I.
Tidal dynamical ejecta. Inferred from numerical relativity simulations, the tidal dynamical ejecta, directly formed by the tidal forces, is highly anisotropic with the mass mainly distributed around the equatorial plane and shaped like a crescent Kawaguchi et al. 2015; Brege et al. 2018) . Due to the lack of weak-interaction process, such dynamical ejecta should be lanthanide-rich which has a relatively low electron fraction Y e 0.2 (e.g., Foucart et al. 2014; Kyutoku et al. 2018) . The mass distribution angle range of the ejecta in the longitudinal direction is typically ϕ d ≈ π, while the half-opening angle in the latitudinal direction is almost in the range of θ d ≈ 10 • − 20 • where we set ϕ d = 15 • as a constant hereafter. Following Kawaguchi et al. (2016) , we simply assume that the mass distributions of dynamical ejecta in the latitudinal and longitudinal direction are both nearly homogeneous, i.e.,
Based on the simulation results from Kyutoku et al. (2015) ; Kawaguchi et al. (2016) , the mass velocity profile has approximately a flat distribution with a rapid decay at low and high velocities, i.e.,
For a dynamical ejecta with a mass of M d , its density can be obtained under these assumptions of the mass distribution. We find the lower positions at the half maximum of the velocity profile are nearly at 0.1 c, so we set the minimum velocity of the dynamical ejecta as v min,d 0.1 c. Summarized from numerical relativity simulation results (see Section 2.3.1), the root-mean-square (rms) velocity of the dynamical ejecta could be approximately expressed as a function of the mass ratio between BH mass and NS mass. When the rms velocity of the dynamical ejecta is known, the maximum velocity can be calculated.
Neutrino-driven wind ejecta. After the remnant disk forms, copious thermal neutrinos due to the high temperature in the disk would drive an ejecta (i.e. a neutrino-driven wind) which in general tend to concentrate towards the polar regions. Due to neutrino illumination which prevents full r-process nucleosynthesis, such an ejecta would be lanthanide-free with a high electron fraction Y e 0.4. Since the neutrino luminosity decreases rapidly, only ∼ 1% of the disk material can contribute to the neutrino-driven ejecta (Fernández & Metzger 2013; Just et al. 2015) . Therefore, we set the mass of the neutrino-driven ejecta as M n = ξ n M disk , where ξ n = 0.01 is a constant and M disk is the disk mass. The simulation results from Martin et al. (2015) revealed that the neutrino-driven ejecta is a polar emission which has a rather uniform mass distribution, i.e., F n,θ (θ) ≈ const 4 , for θ ≤ θ n , where θ n is the half opening angle of the neutrino-driven ejecta in the latitudinal direction. We define θ n ≈ π/6 according to the fitting results for the lightcurves of AT 2017gfo from Perego et al. (2017) . Wollaeger et al. (2018) presented an analytical expansion profile for a spherically symmetric radiation-dominated outflow with an adiabatic index Γ = 4/3 as F (v) = (1 − (v/v max ) 2 ) 3 . We then simply define the neutrino-driven ejecta mass velocity profile along the radial direction as being similar to this velocity profile, i.e., F n,v (v) = (1 − (v/v max,n ) 2 ) 3 , where we set v max,n ≈ 0.2 c based on the simulation results from Martin et al. (2015) .
Viscosity-driven wind ejecta. As the neutrino luminosity decreases, viscous heating and angular momentum transport would play important roles in driving mass loss of the disk. The viscosity-driven ejecta ususally has an intermediate electron fraction, mainly lying in the range of Y e ∼ 0.15 − 0.25 (Fernández & Metzger 2013; Just et al. 2015; Siegel & Metzger 2017; Metzger 2019) 5 . Simulation results indicated that the fraction mass of the generated viscosity-driven ejecta in the remnant disk mass depends on the spin of the remnant BH (Just et al. 2015; Fernández et al. 2015; Fujibayashi et al. 2020 ). More specifically, it can range from ∼ 5% for a low-spin BH to ∼ 30% for a high-spin BH with the dimensionless spin parameter χ BH,f 0.95. Considering that the spin distribution of the remnant BH is almost in the range of χ BH,f 0.6 − 0.9 for a BH-NS system (e.g., Pannarale 2013) whose orbital angular momentum is aligned with the spin of the primary BH, we simply set the mass of the viscosity-driven ejecta as a constant fraction of the disk mass, i.e., Perego et al. (2017) showed that the angular profile of the viscosity-driven ejecta can be well described by F v,θ (θ) = sin 2 θ based on simulations (Wu et al. 2016; Lippuner et al. 2017; Siegel & Metzger 2017) . For velocity profiles, we define a profile similar to the neutrino-driven ejecta, i.e.,
For ejecta with different origin, we assign the opacity according to the value of Y e . Tanaka et al. (2020) presented that the average opacities for the mixture of r-process reactions are nearly a constant for different Y e values, i.e., κ ∼ 20 − 30 cm g −1 for electron fraction Y e ≤ 0.20, κ ∼ 3 − 5 cm 2 g −1 for Y e = 0.25 − 0.35, and κ ∼ 1 cm 2 g −1 for Y e = 0.4. The condition of these assumptions requires temperatures between T = 5 − 10 3 K. A steep decrease of the opacities would significantly occur at the temperature T 1500 K, so that the assumptions can be satisfied during most of the evolution times for a kilonova emission event. Based on these assumptions, we simply set the gray opacity of dynamical ejecta, neutrino-driven ejecta, and viscosity-driven ejecta as κ d = 20 cm 2 g −1 , κ n = 1 cm 2 g −1 , and κ v = 5 cm g −1 , respectively.
The energy that powers a kilonova comes from the radioactive decays of r-process. The specific energy injection rate is approximated by a power law decay, derived from detailed nucleosynthesis calculations (Korobkin et al. 2012; Tanaka & Hotokezaka 2013; Wanajo et al. 2014) :
with˙ 0 ≈ 1.58 × 10 10 erg g −1 s −1 and s ≈ 1.3, where˙ 0 is the thermalization efficiency introduced in Metzger et al. (2010) and Ye is an electron-fraction-dependent term to take into account extremely neutron-rich ejecta with a decay half-life time of a few hours (Perego et al. 2017) , i.e.,
On condition that we have known the profile for each ejecta component and specific energy injection rate, the evolution of temperature and luminosity in the rest frame of the source can be obtained. For an observer, one would see the photosphere defined by the optical depth in a particular direction which is viewing-angle-dependent. Therefore, the observed flux density in principle can be expressed as an integral of the photosphere emission in the frame of an observer:
where ν = (1 + z)ν, h is the Planck constant, k B is the Boltzmann constant. dS proj is the infinitesimal projected photosphere area along the light of sight. We define D and T as the Doppler factor and temperature of the point of this infinitesimal projected photosphere area. We note that the photosphere that observer sees is defined at the "retarded" times for the same observational time t obs , i.e., emission comes from the so-called equal-arrive-time surface. This light propagation effect is fully considered in our calculations (more details see Appendix B in Paper I).
Via the Blandford−Znajek mechanism (Blandford & Znajek 1977) , the remnant BH formed after a BH-NS merger would accrete from the surrounding remnant disk if the NS undergoes tidal disruption and could launch a pair of collimated relativistic jets possibly toward to the polar direction. The relativistic jets would collide with the surrounding interstellar medium to generate a luminous afterglow with emission ranging from radio to X-rays (Rees & Meszaros 1992; Meszaros & Rees 1993; Paczynski & Rhoads 1993; Mészáros & Rees 1997; Sari et al. 1998 ). Barbieri et al. (2019) noted that the kinetic energy of the jet generated after the BH-NS merger could be estimated as
where = 0.015, Ω H is the dimensionless angular frequency at the horizon, which is determined by the final spin of the remnant BH χ BH,f , i.e.,
and f (Ω H ) = 1 + 1.38Ω 2 H − 9.2Ω 4 H is a correction factor for high-spin values (Tchekhovskoy et al. 2010) . The final spin of BH χ BH,f , which depends on the initial BH-NS system parameter, is calculated using Equation (11) from Pannarale (2013)(see also Deng 2020) .
The observations of late-time multiband afterglow associated with GW170817 revealed that the sGRB jet plausibly carries an angular structure (Alexander et al. 2017; D'Avanzo et al. 2018; Gill & Granot 2018; Ghirlanda et al. 2019; Lamb & Kobayashi 2017; Lazzati et al. 2018; Lyman et al. 2018; Margutti et al. 2018a; Mooley et al. 2018; Ren et al. 2020; Troja et al. 2018; Tan & Yu 2020; Xie et al. 2018 ). Here, we consider a power-law jet model (Rossi et al. 2002; Zhang & Mészáros 2002) , where the angular distributions of the kinetic energy and Lorentz factor are adopted as Salafia et al. 2019) :
where we set the core Lorentz factor Γ c = 250, core half-opening angle θ c = 3.7 • , energy power law slope s 1 = 5.5 and Lorentz factor slope s 2 = 3.5 based on the observational constraints of GW170817/GRB 170817A, and the core isotropic equivalent kinetic energy E c = E K,jet /πθ 2 c . The relativistic jet sweeping into the interstellar medium would drive a forward shock. We ignore the possible energy lost caused by the interaction between the jet and ejecta material, because the ejecta for BH-NS mergers is mainly distributed in the equatorial plane which is perpendicular to the moving direction of the jet. For a relativistic shock propagating into a cold ISM, its evolution enters a self-similar stage beyond the deceleration radius. The physical condition of the shocked plasma can be described by the shock-jump conditions (Blandford & McKee 1976) , i.e., the conservation of baryon number, energy and momentum fluxes. For a relativistic shock, due to the equal-arrival-timesurface effect, photons that reached the observer at the same time would emit from different radiation zones at different times. By integrating the equal-arrival-time surfaces, the observed flux of afterglow for a given flux frequency ν is
The surface brightness I ν depends on synchrotron emission of the relativistic electrons accelerated in the shock. The standard synchrotron emission spectra from the relativistic electrons are employed following Sari et al. (1998) . The relation between the photon arrival time t obs and the emission radius can be expressed by the aberration effect:
Photons emitted at (R, θ, ϕ) would arrive at the observer with a time delay with respect to a photon emitted at R = 0. For more details of the afterglow model we used to calculate the sGRB lightcurves along the line of sight, one can see Appendix C in Paper I (for reviews of afterglow theories see Gao et al. 2013b; ).
For a BH-NS merger, whether the NS is directly plunged into the BH or tidally disrupted by the BH can be determined by the the location of the radius at which the tidal disruption occurs R tidal and the relative positions of the radius of the innermost stable circular orbit (ISCO) R ISCO (Kyutoku et al. 2011; . The radius at which tidal disruption occurs, in the Newtonian approximation, can be estimated by balancing the tidal force and the self-gravitational force at the surface of the NS, i.e.,
The normalized ISCO radius R ISCO = c 2 R ISCO /GM BH is given by Bardeen et al. (1972) , i.e.,
with
where χ BH is the dimensionless spin parameter projected onto the orientation of orbital angular momentum for the primary BH. The NS disruption is significant when R tidal R ISCO , which indicates that mass ejection prefers to occur for BH-NS binary parameters with a high BH spin, a small mass ratio, and/or a stiffer NS EoS. Foucart (2012) and Foucart et al. (2018) presented nonlinear models showing that the mass outside the remnant BH after a BH-NS merger can be determined by the relative location between R ISCO and R tidal . By fitting 75 numerical relativity simulation BH-NS events whose NSs are all assumed to be non-spinning, the fitting model of the total remnant mass outside the remnant BH given by Foucart et al. (2018) is
where η = Q/(1 + Q) 2 , M b NS is the baryonic mass of the NS, Q = M BH /M NS is the mass ratio between the BH mass and the NS mass, and C NS = GM NS /c 2 R NS is the compactness of the NS determined by the NS EoS. This fitting formula applies to the total remnant mass for binaries in the following parameter ranges: Q ∈ [1, 7], χ BH ∈ [−0.5, 0.9], and C NS ∈ [0.13, 0.182].
In Paper I, we presented a similar fitting formula for the dynamical ejecta (see also Kawaguchi et al. 2016; Krüger & Foucart 2020 , for other fitting model) by considering 66 NR simulations collected from Brege et al. (2018); Foucart et al. (2014 Foucart et al. ( , 2017 ; Kawaguchi et al. (2015 Kawaguchi et al. ( , 2016 ; Kyutoku et al. (2015 Kyutoku et al. ( , 2018 ; Lovelace et al. (2013) :
This formula covers the ranges Q ∈ [1, 7], χ BH ∈ [0, 0.9], and C NS ∈ [0.108, 0.18].
As mentioned in Paper I, the fitting formula of the total remnant mass (M total,fit ) and the dynamical ejecta mass (M d,fit ) are fitted with independent data, which may make M d,fit approximately equal to or even larger than M total,fit . Therefore, we assume an upper limit on the total remnant mass that can become unbound dynamical ejecta, i.e., M d,max = f max M total,fit , where we set f max = 0.5 based on the numerical relativity results that show the largest ratio of the unbound dynamical ejecta to the total remnant mass being ∼ 50% (the model MS1-Q7a5 in Kyutoku et al. 2015) . The mass of the dynamical ejecta can be expressed as
while the remnant disk mass around the BH can be estimated as
Numerical relativity simulations also provided information regarding the velocity of the dynamical ejecta. The fitting model for the rms velocity of the dynamcial ejecta, given from Paper I, reads v rms,d = (−0.441Q −0.224 + 0.549) c.
It can cover a wide mass ratio range from Q = 1 to Q = 7. The lightcurve of a BH-NS kilonova mainly depends on the amount of the masses of the three ejecta components. In Section 2.2.1, we mentioned that the mass of the neutrino-driven wind ejecta and the viscosity-driven wind ejecta could simply be considered as a constant fraction of the remnant disk mass. We have modeled the dynamical ejecta mass in Equation (18) and the remnant disk mass in Equation (19) which are associated with the BH mass M BH , the BH dimensionless spin parameter χ BH , the NS mass M NS , and the NS EoS. So the parameter distributions of the BH-NS binary systems would determine detection results of the kilonovae. Besides, the kinetic energy of the jet we modeled, affected on the afterglow lightcurve, is dependent on the disk mass M disk and the final spin of the BH χ BH,f . Therefore, the afterglow lightcurve is also determined by BH-NS system parameters. Furthermore, one should also consider the distribution of the viewing angles since both afterglow and kilonova emissions are highly viewing-angle-dependent.
By far, no BH-NS system has been discovered in our galaxy and no confirmed BH-NS GW signal has been reported by LIGO 6 so that we cannot consider BH-NS system parameters based on observations. Here, we assume the following Figure 2 . Distribution of the BH mass from BH-NS systems. The gray bar chart represents the BH mass distribution of population-synthesis simulations α3 with metallicity 0.008 ≤ Z ≤ 0.02 presented in . The blue line is our best fit probability density function of the BH mass distribution. independent distributions of several properties for BH-NS systems 7 based on some theories and simulation results as described below.
1. BH mass -The probability density function of BH mass distribution depends on the progenitor's metallicity Z, the common-envelope parameter α and the natal kick. However, population-synthesis simulations from ; showed that most BH masses in merging BH-NS systems are < 10 M and are mainly distributed in the range from ∼ 5 M to ∼ 8 M for different conditions of the BH-NS binaries progenitors. So we can define a specific probability density function of BH mass distribution. Observationally, the metallicity of sGRB host galaxy spans in the range of 0.6 Z − 1.6 Z where Z = 0.012 (Berger 2009 (Berger , 2014 D'Avanzo et al. 2009 ). Here, in order to obtain our probability density function of BH mass distribution, we assume a formula with an exponential rise and two broken exponential decays (see Fryer & Kalogera 2001; Özel et al. 2010; Steiner et al. 2012 , for an theoretical exponential decay) to fit the α3 model with metallicity 0.008 ≤ Z ≤ 0.02 presented in as our BH mass distribution. The best fitting results of the probability density function showed in Figure 2 is
where M TOV is the NS TOV mass (non-rotating NS maximum mass) and the best fit values of each parameter are a 1 = 0.002845, b 1 = 1.686, a 2 = 1.04 × 10 11 , b 2 = 2.1489, a 3 = 799.1, and b 3 = 0.2904. According to the rapid core-collapse SN model by Fryer et al. (2012) , set a minimum possible mass of a BH as ∼ 5 M , which is the reason of a sharp cut-off in Figure 2 . However, the existence of low-mass BHs cannot be theoretically ruled out (e.g., Drozda et al. 2020) . We extend the minimum mass limit of our model to M TOV , although the region of M TOV < M BH < 5M only occupies a very small portion.
2. BH spin -At present, there is no clear theory to predict the distributions of the magnitude and direction for the BH spin in BH-NS systems. The direction of the BH spin could retain some memory of its birth orientation and be affected by torquing of the accretion disk . It could be plausibly presumed that the orientation of the BH spin is nearly aligned with the orbit. So we assume that the spin of the primary BH is almost aligned with the orbit. The measurement of GWs from binary black holes (BBHs) by LVC showed that Note-The columns from left to right represent the NS EoSs we selected, the NS TOV mass MTOV, the NS radius for 1.4M , and the best fit values of A1 and A2 for each EoS in Equation (23). All of the values of each parameter are cited from Gao et al. (2020) the majority of the observational BBH systems have a low effective spin (Abbott et al. 2019a) 8 , i.e., χ eff ∼ 0, which indicated that the spin of each BH in BBH systems could also be low. The similar result is also obtained from population-synthesis simulations by Bogomazov et al. (2018) Miller et al. 2011 ). Since part of X-ray binaries would evolve into BH-NS binaries, the BHs in BH-NS systems could have high-spin also. The uncertainties in theory encourage us to define two simple cases for the BH spin distribution, i.e., χ BH ∼ N (0, 0.15 2 ) and χ BH ∼ N (0.8, 0.15 2 ), for the low-spin case and the high-spin case, respectively. We assume that χ BH is in the range of χ BH ∈ [0, 0.99]. Table 1 , we select three representative EoSs. From soft to stiff, they are AP4 (Akmal & Pandharipande 1997), DD2 (Typel et al. 2010) , and Ms1 (Müller & Serot 1996) . Gao et al. (2020) presented a relation between NS mass M NS and baryonic mass M b NS which reads
where M b NS and M NS in this formula are in units of M . The best fit values of A 1 and A 2 for each EoS are collected in Table 1 . The baryonic mass as a function of compactness given by Coughlin et al. (2017) is
. Therefore, for each EoS and each NS mass, the compactness of the NS can be expressed as
During the inspiral phase of a BH-NS merger, the matter effect coming from the tidal deformation could enhance GW emission and accelerate the coalescence (Flanagan & Hinderer 2008; Damour et al. 2012) . Therefore, the amount of tidal deformation, which is described by the dimensionless tidal deformalibity of the NS, i.e., Λ = (2/3)k 2 C −5 NS with k 2 is the Love number, depends on the NS mass. So the EoS could affect the GW signal. 4. NS mass -We assume that the distribution of the NS mass in BH-NS systems is approximately the observationally derived normal distribution of galactic BNS systems, i.e., M NS /M ∼ N (1.33, 0.11 2 ) (Lattimer 2012; Kiziltan et al. 2013 ). In addition, the population-synthesis simulation results from and predicted that the masses of most of NSs in BH-NS systems are possibly larger than the NSs in BNS systems, with the median value as ∼ 1.6 M . In order to consider the effect of different NS mass distribution, we assume another normal distribution M NS /M ∼ N (1.6, 0.11 2 ) for specific DD2 NS EoS. The range of NS mass adopted here is M NS /M ∈ [1.1, M TOV ], where the lower boundary 1.1 M is the minimum mass of a NS formed from a core-collapse supernova based on the rapid supernova model of Fryer et al. (2012) .
5. Viewing angles -The distributions of latitudinal viewing angle sin θ view and longitudinal viewing angle ϕ view are adopted to be uniform, so we set sin θ view ∼ U(0, 1) and ϕ view ∼ U(0, π).
Except the redshift distribution, the properties of the observed sGRB afterglow emission are mainly affected by the the isotropic equivalent kinetic energy E c , the power-law index of the electron distribution p, the density of circumburst environments n, the shock microphysics parameters e and B , and the latitudinal viewing angle θ view whose distribution has been obtained. For E c , we assume that it would be dependent on the BH-NS binary system parameters, which has been introduced in Section 2.2.2. We fix p = 2.15 and e = 0.1, since these parameters have relatively narrow distributions from the GRB afterglow modeling (e.g., Santana et al. 2014; Wang et al. 2015) . Therefore, only n and B remain free parameters that can significantly effect the brightness of the observed afterglow emission.
The ranges of n and B distributions are much wider in the literature. By systematically analysing a sample of 38 well-observed broadband sGRB afterglow lightcurves, Fong et al. (2015) found that sGRBs usually occur in low-density environments, with a median density of n ≈ (3 − 15) × 10 −3 cm −3 , and 16% of sGRBs in their sample took place at densities 1 cm −3 . Moreover, O'Connor et al. (2020) found that the medium densities of 16% sGRBs have n 10 −4 cm −3 by utilizing upper limits on the deceleration time and lower limits on the peak X-ray flux with a sample of 52 sGRB X-ray afterglow observations. For the limit of B , Santana et al. (2014) performed a systematic study of 60 sGRBs X-ray afterglows by constraining the observed flux at the end of the steep decay phase. Setting a constant density n = 1 cm −3 , they showed that B has a wide range of ∼ 10 −8 − 10 −3 and is centered at ∼ few × 10 −5 . Similar results were obtained by Wang et al. (2015) . Despite the constraints of B are usually given in the case of a given n, n and B are actually independent physical parameters. So, based on these observational constraints, we simply assume that they have independent log-normal distributions as follow:
1. Density of circumburst environments -We assume its distribution as log(n/cm −3 ) ∼ N (−2.0, 1.1 2 ).
2. Fraction of shock energy carried by magnetic field -We assume its distribution as log B ∼ N (−3.0, 1.0 2 ).
Up to present, except for GW170817/AT 2017gfo that had relatively complete observations, the scarce and ambiguous observational data of presently observed kilonova candidates are not adequate to allow us to model the kilonova luminosity function. Recently, ) used the kilonova flux upper limits of 13 GW events involving at least one NS during the LVC/O3 to constrain the kilonova luminosity function. On one hand, with the continuous upgrades and developments of GW detectors, more BNS and BH-NS events with accurate localization will be detected, which will be helpful for GW follow-up observations to search for kilonovae (e.g., Scolnic et al. 2018; Setzer et al. 2019 ). On the other hand, prospective large wide-field survey projects, e.g. LSST, can discover more kilonovae through blind searches (e.g., Cowperthwaite et al. 2019; Setzer et al. 2019) . A sufficient quantity of kilonovae observed in the future make it possible to constrain the kilonova luminosity function. In this section, we use the assumed parameter distributions and the kilonova model presented in Section 2 to predict theoretical luminosity function of kilonova emission from BH-NS mergers.
Kilonovae are thought to be fast-evolving with lightcurve peak times distributed from 0.2 day to 1 day (e.g., Korobkin et al. 2020) . Early fast evolution of kilonovae makes it difficult to observationally catch the lightcurve peak times and reconstruct the bolometric lightcurve near the peak. With the consideration of actual kilonova searches, we adopt the absolute magnitude at 0.5 day post-merger, e.g. the first detection times of GW170817/AT 2017gfo (Abbott et al. 
However, the luminosity function Φ(M ) cannot show the relative fraction of the case where the NS directly plunges into the BH without a kilonova out of all BH-NS mergers by considering the distributions of different BH-NS system parameters, so we look for other alternative parameters to display Φ(M ). We define the local event rate density per unit absolute magnitude for kilonova emission from BH-NS mergers aṡ
where M represent the absolute magnitude at 0.5 day after the merger. Therefore, the number density in the absolute magnitude interval from M to (M + dM ) can be obtained by Equation (1), i.e.,
Define the number per unit time is ∆Ṅ in a finite absolute magnitude bin from M to (M + ∆M ). We then havė
The local event rate density is thereforė
where ξ is introduced to consider the fraction of BH-NS mergers that can produce kilonovae, M min and M max are the minimum and maximum values of the absolute magnitude, respectively. With the definitions in Equation (25) and (29), the kilonova luminosity function can be derived by displaying the local event rate density per unit absolute magnitudeρ 0,M /ρ 0 ∝ Φ(M ) as a function of M . In view of the relation between absolute magnitude M and bolometric luminosity L, i.e., ∆M = 2.5∆ log L, one can also mapρ 0,log L /ρ 0 as a function of L to replaceρ 0,M /ρ 0 . Similar to the concerned absolute magnitude, L should be the kilonova bolometric luminosity at 0.5 days after a BH-NS merger.
In Section 2.3, we presented the variable parameter distributions for the NS EoS, NS mass, and BH spin in BH-NS systems. Hereafter, we consider 6 groups of parameter distributions containing a low-mass NS, 3 kinds of NS EoSs (AP4, DD2, and Ms1), and 2 kinds of BH spin distributions (Low χ BH and High χ BH ). In addition, for the DD2 EoS, we assume another two sets of massive NS distributions, so we make a total of 8 groups of BH-NS system parameter distributions to simulate the kilonova luminosity functions. Each group of simulation is named following the convention of 'EoS−NS mass distribution−BH spin distribution' which is shown in the label of Figure 3 .
Following Equation (28), we use the Monte Carlo method to simulate 1 × 10 6 BH-NS mergers in the universe and map the relevantρ 0,M /ρ 0 as a function of M (ρ 0,log L /ρ 0 as a function of L) shown in Figure 3 for these eight different distributions of system parameters. Altogether, the results indicate that the kilonova absolute magnitude M at 0.5 days postmerger for BH-NS mergers is mainly distributed in the range of ∼ −10 to ∼ −15.5 (or ∼ 3×10 39 erg s −1 to ∼ 5 × 10 41 erg s −1 for bolometric luminosity L). Since 10 41 erg s −1 is roughly 10 3 times of the nova luminosity, and hence, the terminology 'kilonova' was invented (Metzger et al. 2010) , our results actually suggest that these events can be 'hecto-' (100 times) or even (multi-) 'deka-' (10 times) novae depending on parameters.
The diversity of different system parameter distributions could be directly reflected in the distribution of the absolute magnitudes. Obviously, among these three variable parameters, the BH spin distribution has the largest impact on the kilonova luminosity distribution. As mentioned in Section 2.3.1, whether the NS is directly plunged into the BH or tidally disrupted by the BH is determined by the comparison between two radii: one is the radius at which the tidal disruption occurs R tidal (which is related to NS mass M NS and NS compactness C NS ), and the other is the ISCO radius The blue, yellow and orange lines represent three selected EoSs, i.e., AP4, DD2, and Ms1, with the Galactic NS mass (low MNS) distribution, while the green lines represent DD2 with high NS mass MNS distribution predicted by and .
R ISCO (which is dependent on BH mass M BH and dimensionless spin parameter χ BH ). If R tidal R ISCO , the NS would plunge entirely into the BH without kilonova emission. If R tidal R ISCO , the NS would be partially disrupted to so that neutron-rich materials can be ejected to power a bright kilonova. For a fixed NS mass and EoS, the BH with a lower spin has a larger R ISCO which cause the NS directly plunged for most cases. This is why the overall luminosity functionsρ 0,M /ρ 0 for low BH spin distribution shown in Figure 3 have lower values compared to those for high BH spin distribution. Tidal disruption would be significant if the BH has a higher spin, so that more material could be ejected to power a brighter kilonova. For fixed BH parameters, a smaller NS mass and a stiffer EoS are required to have tidal disruption. Because a less massive NS is less compact and because a stiffer EoS can support a larger NS radius, introducing a larger R tidal makes it easier to tidally disrupt the NS and release more neutron-rich material. As shown in Figure 3 , BH-NS kilonovae tend to be more luminous if the NS has a lower mass with a stiffer EoS.
One can conclude that the kilonova absolute magnitude at 0.5 days after a BH-NS merger is mainly distributed in the range of ∼ −10 to ∼ −15.5, corresponding to the bolometric luminosity L from ∼ 3 × 10 39 erg s −1 to ∼ 5 × 10 41 erg s −1 . The kilonova luminosity function is highly dependent on the BH-NS system parameter distributions. In the case of less massive NSs, stiffer EoSs, and higher BH spins in particular, NSs can be tidally disrupted significantly while BH-NS merger kilonovae have more luminous distributions.
Observers tend to model the luminosity function for a certain type of transient by its observed luminosity. In Paper I, we showed that the variation of the observed luminosity for radioactivity-powered kilonova emission from BH-NS mergers is mainly caused by the change of the projected photosphere area of the dynamical ejecta, which causes a difference of only a factor of ∼ (2 − 3) of the observed luminosity between the face-on view (θ view ∼ 0 • ) with the largest projected photosphere area and the edge-on view (θ view ∼ 90 • ) with the smallest projected photosphere area. Since kilonova emission is highly viewing-angle-dependent, its luminosity function should contain information of the viewing angle distribution.
In Figure 4 , we select two examples of theoretical luminosity function with the consideration of two different BH-NS system parameter distributions from Figure 3 . The example luminosity functions contain the potential viewing angle distribution information of anisotropic kilonova emission. To see how viewing angles contribute to the final luminosity function, we also separately mark the contributions from latitudinal viewing angles θ view < 45 • and θ view ≥ 45 • , respectively, in Figure 4 . One can see that the kilonova luminosity function is a superposition of a high-luminosity subset with smaller latitudinal viewing angles and low-luminosity subset with larger latitudinal viewing angles. The ratio between the peak luminosities of the two subsets is ∼ 1.4, corresponding to the difference in absolute magnitude by ∼ 0.5 mag.
If there is a long-lived energy injection from the remnant BH, e.g., fall-back accretion (Rosswog 2007; Kyutoku et al. 2015) or the Blandford-Payne mechanism (Blandford & Payne 1982; Ma et al. 2018) , the luminosity function of BH-NS kilonova would be more complex. It can be predicted that BNS kilonova luminosity function can be more interesting due to a larger difference of the observed luminosity caused by the viewing angle variation, different remnant compact objects (Gao et al. 2016; Kawaguchi et al. 2020) , and possible long-lived energy injection . et al. 2017a,b,d, 2018b, 2019b,d) . The follow-up observations using space-borne and ground-based telescopes led to the discovery of the broadband EM signals from GW170817. This event was first case of detecting EM counterparts of a GW event with target-of-opportunity follow-up observations of GW triggers (Abbott et al. 2017b). The LVC released the GWTC-1 GW catalog after the end of O2 (Abbott et al. 2019a). This catalog contains 10 BBH GW signals and 1 BNS GW signal with high significance but no BH-NS GW signal in O1 and O2. Although no significant BH-NS GW signal was announced, LIGO officially set an estimated upper limit on the BH-NS event rate using O1 data (Abbott et al. 2016e) 9 . After further upgrades, LVC started O3 in April 2019 and ended in March 2020 10 . At present, LVC has announced a total of 56 GW candidates during O3 11 , four of them were first confirmed to be high significance gravitational wave events (Abbott et al. 2020a,d,e,f) . After that, LVC released a list of high significance gravitational wave events in O3a, named GWTC-2 (Abbott et al. 2020g). Among these GW events, GW190425 is considered the second BNS event. A low significance BH-NS candidate GW190426 was reported in GWTC-2, but its astrophysical probability is strongly dependent on prior assumptions of the rate of BH-NS GW signals. The mass of the lighter object in GW190814 lies between the heaviest known neutron star and the lightest known black hole, which is in the so-called "mass gap". The tidal deformation parameters obtained through parameter estimation are not sufficient to confirm its composition (Abbott et al. 2020e), so LVC cannot definitively call it a BH-NS event. As of the writing of this paper, LVC has not announced any high significance BH-NS GW event. However, it is worth noting that among the 56 announced GW candidates, the point estimation parameters 12 of several candidates meet the requirements of a BH-NS merger . These events require further confirmation by LVC. Despite great efforts of follow-up observations, no EM signal was detected during O3 Antier et al. 2020a,b; Coughlin et al. 2020a,b; Gompertz et al. 2020; Kasliwal et al. 2020) . It remains an open question why no associated kilonova or afterglow was identified from these BH-NS merger candidates.
With continuous upgrade of existing detectors and increase in the number of GW detectors, the chance of detecting BH-NS GW signals in the near future is increasing. So far no one has made a detailed calculation on the detection capability of different detector networks in the near future for BH-NS signals. This section is to present a detailed calculation of the BH-NS GW detectability of the 2nd, 2.5th and 3rd generation GW detectors in the next 15 years. We will present the detector networks used in this work, the parameters of the simulated GW signals, optimal signalto-noise ratio (SNR), and the calculation methods of the GW signal duration. Together with the upgrade of GW detectors, several large EM survey projects will start operation and join follow-up observation collaborations that can significantly improve the search ability for target-of-opportunity observations of GW triggers. We also simulate the follow-up EM detectability at different epochs of GW detector networks. In addition, some implications for the past and future GW-triggered target-of-opportunity observations are given in this section. The key feature of matched filtering is to convolve the detector's strain data with the GW templates. Inaccurate templates will cause loss of SNR. When LVC estimates model parameters of the detected GW candidates using Bayesian methods (such as Markov Chain Monte Carlo and nested sampling), they need to calculate a large number of likelihood functions (Thrane & Talbot 2019) . such calculations of the likelihood function require accurate GW templates.
In this work, we need to predict the detection capabilities of different GW detector networks in the foreseeable future for BH-NS systems. Accurate GW templates are needed to make our prediction results accurate. At present, the most mature GW templates of compact binary coalescences are the BBH templates. For systems containing NSs, such as the BNS systems and the BH-NS systems, due to the tidal deformation of the NS, the corresponding GW templates are more complicated. There are currently several GW templates specially used for BNSs, such as the post-Newtonian approximation templates TaylorF2 (Sathyaprakash & Dhurandhar 1991; Poisson 1998; Mikoczi et al. 2005; Arun et al. 2009; Bohé et al. 2013 Bohé et al. , 2015 Mishra et al. 2016) , the phenomenological models IMRPhenomD NRTidal Khan et al. 2016; Dietrich et al. 2017 Dietrich et al. , 2019 , IMRPhenomPv2 NRTidal ) and its improvement IMRPhenomPv2 NRTidalv2 , as well as effective-one-body numerical relativity (EOBNR) models such as SEOBNRv4T (Bohé et al. 2017; Lackey et al. 2019; Hinderer et al. 2016; Steinhoff et al. 2016) , TEOBResumS (Damour & Nagar 2010; Vines et al. 2011; Damour et al. 2012; Bernuzzi et al. 2015; Nagar et al. 2018) , SEOBNRv4 NRTidal (Bohé et al. 2017; Dietrich et al. 2017 ) and its improvement SEOBNRv4 NRTidalv2 (Bohé et al. 2017; Dietrich et al. 2019 ). As the current numerical relativity simulations of BNS mergers are not perfect, the aforementioned templates only contain the GW waveforms during the inspiral phase (Bernuzzi et al. 2014; Dietrich et al. 2019) . GW templates specifically used for BH-NS systems are also under development, such as the recently developed IMRPhenomNSBH (Thompson et al. 2020) and SEOBNRv4 ROM NRTidalv2 NSBH (Matas et al. 2020 ). These two templates are based on the phenomenological model IMRPhenomD Khan et al. 2016 ) and the effective-one-body numerical relativity model SEOBNRv4 (Bohé et al. 2017) with added amplitude correction (Pannarale et al. 2015) and phase correction caused by the tidal deformation of the NS. They include GW of inspiral, merger and ringdown. However, these two templates were not publicly available when our work started, so these two waveform templates specially used for BH-NS mergers are not used in this work.
In this work, we use the phenomenological model IMRPhenomPv2 NRTidal. It is based on the GW template IMRPhenomPv2 (Hannam et al. 2014) , which describes GWs from quasicircular precessing binary black holes. Phase correction caused by the tidal deformation of the NS is introduced through the NRTidal model ). The NRTidal model itself is tuned with the numerical relativity waveform. The IMRPhenomPv2 NRTidal template was used in the parameter estimation of GW170817 (Abbott et al. 2019c). It is a relatively mature template. Although the NRTidal model itself is designed for BNSs, it has recently been proven that it can be also applied to BH-NS mergers . Since this article needs to simulate millions of BH-NS GW signals, the calculation time of the waveforms is also a factor of consideration. The frequency domain phenomenological model IMRPhenomPv2 NRTidal meets this requirement. Here we ignore the influence of higher-order modes (Ruiz et al. 2008) due to the large mass ratios, and only consider the GW signals during the BH-NS inspiral phase. We also ignore the precession effect. At present, ground-based GW detectors include two Advanced LIGO detectors in the United States (Harry et al. 2010; Aasi et al. 2015) , Advanced Virgo detector in Europe (Acernese et al. 2014) and KAGRA detector in Japan (Aso et al. 2013; Akutsu et al. 2018) . During O1 staring in 2015, only the two LIGO detectors in the United States have completed the upgrade from the 1st generation of GW detectors to the 2nd generation of GW detectors (Smith et al. 2009 ). So only LIGO Hanford (H1) and LIGO Livingston (L1) made observations (Martynov et al. 2016) , even though both detectors at the time did not reach the design sensitivity of the 2nd generation GW detectors (Martynov et al. 2016) . In O2, H1 and L1 have further improved the sensitivity over O1, which is closer to the design sensitivity of Advanced LIGO. At the end of O2, the Advanced Virgo (V1) finished upgrading from the 1st generation to the 2nd generation, formed the 2nd generation GW detector network for the first time. Due to the increase in the number of detectors during O2 and the increase in detectors' sensitivity, the number of GW events observed during O2 is far more than that in O1 (Abbott et al. 2019a). After a long time upgrade (McCuller et al. 2020; Acernese et al. 2019) , the Advanced LIGO and Advanced Virgo began their third observation run (O3). At the end of O3, the Japanese detector KAGRA (K1) and the German detector GEO600 (Dooley et al. 2016 ) conducted a 2-week joint observation. Up to now, the four 2nd generation GW detectors have not achieved their respective design sensitivities yet, but it is expected that in the fourth observation run (O4), LIGO, Virgo and KAGRA will reach their respective design sensitivity (Abbott et al. 2018a ). In the following, we refer to the time era during which these four detectors reach the design sensitivity as the "HLVK era", which is expected to begin in 2021. The detector sensitivities during this era, expressed as the amplitude spectral densities (ASD), are presented in the top left panel of Figure 5 , which are adopted from the public data 13 . We also plot GW190814's sensitivity curves 14 as representatives for O3 observations. As the location of the third LIGO detector LIGO-India (Unnikrishnan 2013) has not been determined yet, this work does not consider LIGO-India. As the 2nd generation detectors are about to reach the design sensitivity, people have begun to study how to further improve the sensitivity on the basis of the existing 2nd generation detectors. These upgrades are called 2.5th generation detectors. For example, the subsequent upgrade of Advanced LIGO is called Advanced LIGO Plus (A+) , and Advanced Virgo also has Advanced Virgo Plus (AdV+) upgrade plan (Abbott et al. 2018a). The current KAGRA detector is called baseline KAGRA (bKAGRA) (Aso et al. 2013 ), but within 5 to 10 years, it will also be upgraded to KAGRA+ (Michimura et al. 2020 ). KAGRA+ has multiple upgrade plans. For convenience, we select the sensitivity of the long-term upgrade plan "KAGRAplusCombined" (Michimura et al. 2014 ) as the sensitivity of KAGRA+. The detector network composed of these 2.5th generation detectors is referred to "PlusNetwork" in the following. The sensitivity curves are shown in the top right panel of Figure 5 , which are from public data 15 . After 2030, the 3rd generation GW detectors are expected to replace the 2nd generation and the 2.5th generation detectors. The 3rd generation detector plans currently proposed include the Einstein Telescope (ET) in Europe (Punturo et al. 2010a,b; Figure 5 . It can be seen from these sensitivity curves that the sensitivity of LIGO Voyager is much lower than that of ET and CE, so we will discuss the situation of LIGO Voyager separately. Since the location of ET and CE has not yet been determined, according to the convention (Vitale & Evans 2017; Vitale & Whittle 2018), we place ET at the current Virgo detector position and CE at the current H1 position. Since LIGO Voyager is directly using the existing LIGO site, we place two LIGO Voyagers at the sites of H1 and L1 respectively. The sensitivity curves of ET and CE used in this paper come from the official websites 1617 , noting that we have selected the sensitivity of "ET-D" and "CE1". LIGO Voyager's sensitivity curve comes from the public data 18 . We set the low cut-off frequency (f min ) of all the 2nd and 2.5th generation detectors to 10 Hz , 20 Hz for O3. For the 3rd generation gravitational wave detectors, we set f min of LIGO Voyager to 10 Hz (Adhikari et al. 2020), f min of CE to 5 Hz (Reitze et al. 2019 ), and f min of ET to 1 Hz (Punturo et al. 2010a) . From the bottom panel of Figure 5 , we can see the f min value of each detector.
We simulate multiple sets of GW signal parameters according to different degrees of NS tidal deformation, BH spin and NS mass according to the method presented in Section 2.3.1. According to these parameters, we can obtain the redshifted masses in the geocentric coordinate system. By using the IMRPhenomPv2 NRTidal template, we can simulate the GW waveform in the corresponding geocentric coordinate system, and then project it to different detectors to obtain the strain signal detected by each detector. The optimal SNR is calculated as
where S n (f ) is the one-sided power spectral density (PSD) of the GW detector, which is the square of ASD. Here f max is set to 2048 Hz, andh(f ) is the strain signal in the frequency domain. We can use the optimal SNR to approximate the matched filtering SNR of the GW signal detected by each detector, and then we can calculate the network SNR of the entire detector network, i.e., the square sum of the SNR of all detectors. When the network SNR is greater than the threshold 8, we expect that the corresponding GW signal is detected. In order to consider the best and worst situations in different eras, we consider "all detectors in the corresponding era have reached the designed sensitivity and work normally" as the optimal situation. In the mean time, we consider "only the detector with the worst sensitivity in the corresponding era works" as the worst situation. For the 2nd generation detector network, the best case is that H1, L1, V1 and K1 all work normally, we abbreviate it as "HLVK", while the worst case is that only K1 works normally.
In the same way, the best case for the 2.5th generation detector network is that A+, AdV+ and KAGRA+ all work normally, which we abbreviate as "PlusNetwork", and the worst case is that only AdV+ works. We discuss LIGO Voyager separately. "ET+CE" represents the best case of the 3rd generation era, and a single ET is considered as the worst case. The multi-messenger astronomy brought by GW170817 emphasizes the joint observations of GW and EM signals. Since the GW from inspiral of compact binaries will reach the earth earlier, early warning through GW to the EM community could be a feasible way of detecting more EM counterparts. There have been preliminary researches of early warning Sachdev et al. 2020; Kapadia et al. 2020; Singh et al. 2020 ). We will calculate the duration of BH-NS GW signals in different detector networks. This duration can be regarded as the upper limit of the early warning time 19 . Here we use the 3.5 post-Newtonian order expression for the chirp time T chirp (Allen et al. 2012) 
with v low = GM c 3 πf min
where M and η are the total mass M = m 1 + m 2 and the symmetric mass ratio η = m 1 m 2 /M 2 , m 1 and m 2 are the component masses of the BH-NS system, and γ E ≈ 0.577216 is the Euler constant. In the single detector case, f min is the lower cut-off frequency of the GW detector. In the GW detector network case, f min is set to the smallest f min among all the detectors. Note that we ignore the spin effect and EoS of NS in the calculation of the chirp time T chirp . The results of our simulations are shown in Section 4.1.3.
In this section, we show the simulated detection results of GW signals for 8 sets of BH-NS system parameter distributions. For each set of parameter distributions, we simulate 1 × 10 6 BH-NS GW signals. Here, by settinġ ρ 0 = 35 Gpc −3 yr −1 , we take the two sets of simulation results of DD2−Low M NS as an example to illustrate the detectability of BH-NS signals with GW detectors in different eras. The top panel of Figure 6 shows the detectability by the detector networks for the DD2−Low M NS −Low χ BH BH-NS signals. One can see that for the 2nd generation GW detector network, the worst case is that only bKAGRA is running so that the observed BH-NS detection rate 19 The matched filtering SNR will increase with the square root of signal duration, and the actual detection requires the matched filtering SNR to be higher than a certain threshold, so it must be shorter than the total duration of the GW signal. Figure 6 . Detectability of BH-NS mergers by various detector networks for the DD2−Low MNS−Low χBH (top panels) and DD2−Low MNS−High χBH (bottom panels) BH-NS signals in different GW detection eras. The left panels show the signals that can be detected by different detectors and detector networks (optimal SNR greater than 8). The crimson and red dots represent the detectable signals for K1 and AdV+, while the orange, yellow, green, blue, and violet dots are the GW signals detectable at HLVK, PlusNetwork, LIGO Voyager, and ET+CE networks eras, respectively. The gray dots in the background represent undetectable signals. The small panels in the upper right corner are enlarged images of the low redshift area. The middle panels show the distributions of all simulated signals on the "optimal SNR-redshift" plane, with a detection threshold. The right side of the threshold are the GW signals that can be detected. The upper right corner is the converted annual detection rates and the total event rate where we define the local event rate of BH-NS mergers asρ0 = 35 Gpc −3 yr −1 . The right panels show the distributions of BH-NS detection rates with redshift, and the insets are zoom-in pictures in the low redshift region.
is ∼ 5 yr −1 . The best case is that the entire HLVK network is observing so that the detection rate is ∼ 80 yr −1 . For the 2.5th generation detector network, the worst case is that only AdV+ is in operation, with a detection rate of ∼ 30 yr −1 . This is better than the result of a single 2nd generation GW detector, but not as good as the HLVK detector network. The optimal situation is that PlusNetwork is fully operating, with a detection rate of ∼ 650 yr −1 . For the 3rd generation detector network, LIGO Voyager can increase the detection rate to ∼ 4, 800 yr −1 . This result is better than the 2nd generation and 2.5th generation detector networks, but is far lower than the newly designed 3rd generation detectors. For a single ET, the detection rate is ∼ 2.6 × 10 4 yr −1 , while for the ET+CE network, the detection rate is ∼ 2.84 × 10 4 yr −1 . These results are very close to the total BH-NS event rate (∼ 2.87 × 10 4 per year for our models), indicating that ET and ET+CE can basically detect almost all BH-NS GW events in the entire Universe. As shown in right panels of Figure 6 , as the sensitivity of the GW detector increases, BH-NS events at higher redshifts gradually dominate the detected event rate. The event rate of ET and ET+CE is mainly dominated by the BH-NS mergers at z ∼ 1. For comparison, when we set BH spin to the high-spin prior, the simulation results are shown in the bottom of Figure 6 . Comparing with the top panel, we find that high spin of the black holes can increase the detection rate of BH-NS mergers, but the increase is less than an order of magnitude. The results of the two sets of simulations are basically similar. The reason for this result is that when other parameters are fixed, the increase of spin will make the GW waveform longer in the time domain, and the SNR of the matched filtering is proportional to the square root of the GW signal duration, so the SNR of GW signal becomes larger, and the signal itself becomes easier to detect. The simulation results of the other 6 groups are in the Appendix B. Through comparison, we can find that the simulation results of each group are the same in order of magnitude, but the effects of each variable on the detection rate of BH-NS mergers are different. The improvement brought by mass of NS is larger than by black hole spin, because the amplitude of the GW signal is proportional to chirp mass (Allen et al. 2012) , and the NS EoS has almost no effect on the detection rate of BH-NS mergers. All of our simulation results for BH-NS merger GW detection rates are summarized in Table 1 . Besides, the median detectable distances (redshift) with a 90% detectable interval, and the maximum detectable distances (redshift) for our simulated GW siginals are given in Table 2 . Our results for "HLVK era" are generally consist with LVC's prospect results in O4 (Abbott et al. 2020b).
We calculated the T chirp of DD2−Low M NS and DD2−High M NS BH-NS simulations, as shown in Figure 7 . All the BH-NS signals that can be detected by each detector or detector network are drawn in the figure. In top panels of Figure 7 , we can see the general trend that the chirp time decreases as the geocentric frame total mass increases. This is because the GW frequency of the ISCO of the binary star system is inversely proportional to the total mass (Wu et al. 2020 ). The larger the total mass, the closer the ISCO frequency is to f min . For the 2nd and 2.5th generations of GW detectors, although f min is 10 Hz, due to the increase in detector's sensitivity, more short GW signals can have sufficiently large SNR, that is, large mass BH-NS is easier to detect, which is why the median value of T chirp on the right side of each panel in Figure 7 decreases as the sensitivity of the detector increases. In the bottom panels of Figure 7 , the T chirp of ET and ET+CE are shown. For ET and ET+CE, the f min is reduced from 10 Hz to 1 Hz, which greatly increases the signal duration. All T chirp values are shown in Table 3 . From the table, we can see that for the 2nd and 2.5th generation detectors and LIGO Voyager, T chirp does not exceed 10 minutes, which means that the time for early warning is very limited. For ET and ET+CE, T chirp has reached the order of hours or even a few days. There will be enough preparation time for EM observations in these cases.
Target-of-opportunity observations after GW triggers would be an optimal strategy to search for associated EM signals, especially for kilonova, in the future. The BH-NS GW detectabilities for the networks of the 2nd, 2.5th, and 3rd generation GW detectors have been studied in detail in Section 4.1.3. Based on these results, we now discuss the kilonova detection probabilities for GW-triggered target-of-opportunity observations under ideal follow-up conditions in this section. Gravitational Wave High-energy Electromagnetic Counterpart All-sky Monitor (GECAM), which will be launched in later 2020 can cover all-sky to detect GRBs (e.g., Song et al. 2019) . GW-associated GRBs would essentially not be missed and on-axis jet afterglows and kilonovae of these events will likely be detected. In the following, we do not focus on these joint GW/sGRB detections but focus on GW triggers without an sGRB association. In the absence of an sGRB, as-yet undiscovered off-axis jet orphan afterglow (e.g., van Eerten et al. 2010 ) could establish an EM link to a GW trigger. Therefore, the detectability of orphan afterglow from BH-NS mergers for target-of-opportunity observations is also predicted together with kilonovae. In fact, the observed optical EM flux should be a superposition of the emissions from both the kilonova and the afterglow, so we also give the detection probabilities of the total emission in this section.
As Advanced LIGO and Virgo will upgrade while additional interferometers (KAGRA and LIGO-India) come online in the future, the detection horizon would be several hundred Mpc for BNS and BH-NS mergers with the localization accuracy reaching to ∼ 10 deg 2 by a network of these 2nd GW detectors (Abbott et al. 2018a ). In contrast, higher localization accuracy would be obtained by networking of 3rd generation GW detectors. For example, a network of three or four 3rd generation detectors can improve sky localization for most of BNS or BH-NS mergers to be smaller than ∼ 1 deg 2 up to z ∼ 0.5 or ∼ 10 deg 2 up to z ∼ 2 (Zhao & Wen 2018; Hall & Evans 2019; . In these regimes, present and future wide-field-of-view survey projects can be able to cover sky localization given by GW detections in few pointings and achieve deep detection depths with relatively short exposure integration times. The detection of GW170817/AT 2017gfo was made in an i-band filter and confirmed by larger telescopes (Coulter et al. 2017) . Similar to this case, we assume that one can use only one filter to make follow-up searches after GW triggers. With quick response times of GW events follow-up observations, we also assume that the survey telescopes would not miss the peaks of kilonovae in each band. Therefore, by assuming that all GW detectors in each era can operate normally, for a given X-band limiting magnitude m X,limit , the detection rates for target-of-opportunity observations can be expressed aṡ
where z max is the maximum detectable distance andρ <m X,limit (z) represents the event rate densities of the BH-NS mergers that can be triggered by GW detectors and be detected by telescopes in X band. Here, the astronomical magnitude system we adopted is the AB magnitude system, i.e., m ν = −2.5 log(F ν /3631 Jy).
We consider specific BH-NS system parameter distributions for the DD2 − Low M NS − High χ BH model as an example. Figure 8 shows the simulated detection fractions of kilonova, orphan afterglow, and total emission as a function of search limiting magnitude in three example observed filter bands, i.e., riJ bands. If the contribution from the afterglow is ignored, in the same search limiting magnitude for each filter, one can discover more kilonovae using redder bands since BH-NS merger kilonovae are optically dim but possibly infrared bright (Paper I). More specifically, kilonovae discovered in J band are almost more than double of those discovered in r band. However, in view that most survey projects always have a higher sensitivity in the optical band and that BH-NS merger kilonovae lightcurves in the optical band evolves faster and reaches the peak with a shorter time, it seems that there is no significant difference between searching in optical and in infrared. Because the kilonova reaches the peak earlier in optical, a good strategy would be to follow up early in blue bands to learn the early-stage emission of the kilonova.
In Paper I, we indicated that the AB absolute magnitudes of BH-NS merger kilonova are fainter than ∼ −15 mag and ∼ −16 mag in the optical and infrared bands, respectively. As shown in Figure 10 , since the search limiting magnitudes of present and planned wide-field survey projects are m limit 26 mag, almost all detectable BH-NS merger kilonovae would be discovered at a distance of z 0.3 which can cover the GW detection space in the HLVK era. For the GW events that are detected during the future eras except the HLVK era, one can hardly find more kilonovae by blind-search survey projects. As a result, there appears to be little difference in the amount of detectable kilonovae at different GW detection eras as shown in Figure 8 . Even so, for those BH-NS mergers whose jets beam toward us or slightly deviate from us, their sGRBs and afterglows could still be detected. For the wide-field survey projects whose search limiting magntidues are relatively shallow, Figure 8 reveals that it has a low probability of discovering BH-NS merger kilonovae after GW triggers. However, these survey projects can find bright on-axis jet afterglows by following up detections of sGRBs or orphan afterglows with a relatively higher probability than detecting kilonovae. Their relatively larger fields of view can make them more efficiently to localize the source, facilitating continued follow up observations using large telescopes. In addition, the detection fractions of orphan afterglow and total emission Figure 8 . Detection fractions as functions of 5σ limiting magnitude and search filter for target-of-opportunity observations during four GW detection eras, i.e., HLVK, PlusNetwork, LIGO Voyager, and ET+CE eras. The solid, dashed-doted, and dashed lines denote the detection fractions of kilonova, orphan afterglow, and total emission, respectively. The blue, green, and red lines represent the detection fractions in r, i, and J bands, respectively. The right axis of each panel represent the detection rates by considering the local BH-NS merger event rate asρ0 = 35 Gpc −3 yr −1 . The example BH-NS system parameter distributions are adopted for the DD2−Low MNS−High χBH model. would rise rapidly with the increasing search limiting magnitude and the continuous upgrading of the GW detectors. Although afterglow emission is beamed, it is the most probable detected EM signal in association of BH-NS GW signals in the future.
Since the detection fractions of kilonova, afterglow, and total emission are essentially not influenced by the search filter, we only show the r-band detection fractions for different BH-NS system parameter distributions in Figure 9 . Similar to the discussion for the kilonova luminosity function in Section 3.1, the detection fractions are also highly dependent on the BH-NS system parameter distributions, especially for the BH spin parameter. The BH-NS mergers can usually generate bright kilonovae and afterglows if the spins of the primary BHs are high. This would be beneficial for target-of-opportunity follow-up observations. However, for those GW-triggered BH-NS mergers that have a massive NS component with soft EoS and a low-spin primary BH, it is barely possible to detect their EM counterparts. Our simulated detection fractions and detection rates of the total emission for eight BH-NS system parameter distributions are shown in Table 5 (see Appendix C).
Despite great efforts, follow-up observations to O3 LVC events revealed no creditable BH-NS merger kilonova. So, why did this happen? One possibility is that there was actually no EM counterpart produced from most BH-NS Figure 9 . r-band EM detection fractions of GW-triggered events as a function of r-band 5σ limiting magnitude. The panels from left to right represent the detection fractions of kilonova, orphan afterglow, and total emission, respectively. The panels from top to bottom denote GW detection at four eras, i.e., HLVK, PlusNetwork, LIVO Voyager, and ET+CE eras. The right axis of each panel represent the detection rates by considering the local BH-NS merger event rate density asρ0 = 35 Gpc −3 yr −1 . The line convention is the same as that in Figure 3 . Figure 10 . Left panel: the maximum search depth as a function of 5σ limiting magnitude for different peak absolute magnitudes of BH-NS merger kilonovae (blue lines). The gray lines represent the r-band 5σ limiting magnitude of present and planned survey projects, including ZTF, Mephisto, PS1, LSST, and CSST, with different exposure times. Right panel: GW detection rates per unit redshift depending on distance for four GW detection eras. The example BH-NS system parameter distributions are adopted for the DD2−Low MNS−High χBH model. Figure 11 . Comparison between the detectabilities for the network HLV (O3, striped histograms), for one single K1 (carmine histograms), and for the network of 2nd generation GW detectors (orange histograms). Here we assume the example BH-NS system parameter distributions for the DD2−Low MNS−High χBH model. mergers. If the BH-NS system distributions tend to have low-spin primary BHs and massive NSs with a soft EoS, for a large portion of BH-NS mergers the NSs would directly plunge into the BHs without being tidally disrupted, so that no sGRB or kilonova is expected from these systems 20 . Another possibility is that we could not detect the EM counterparts because the search limiting magnitude for follow-up survey observations is too shallow.
Our simulation results show that the median detectable distance of targeted BH-NS mergers GW events for a single and a network of 2nd generation GW detectors are ∼ 300 Mpc and ∼ 700 Mpc, respectively (see Table 2 ). For comparison, Figure 11 shows that the detection rate and detectable distance for HLV (O3) are approximately the a Area survey rate for each survey is calculated with the assumptions of an average 6 hr observing time per night and two exposure bands for each visit. Since Mephisto can achieve simultaneously imaging in three bands, we calculate its area survey rate by assuming single exposure for each visit. same as those for the case that only a bKAGRA is running. This is basically consistent with the detection rate and the distance distribution of BH-NS merger candidates detected during LVC O3 (e.g., Anand et al. 2020; Antier et al. 2020b; Coughlin et al. 2020b; Gompertz et al. 2020; Kasliwal et al. 2020 ). In Section 3, we present that the kilonova absolute magnitude at 0.5 days after a BH-NS merger is mainly distributed in the range of ∼ −10 to ∼ −15.5. In view that the limiting magnitude of the follow-up wide-field survey projects is almost 21 mag (e.g. Antier et al. 2020b; Gompertz et al. 2020; Coughlin et al. 2020b; Kasliwal et al. 2020; Wyatt et al. 2020) , the maximum detectable distance for BH-NS kilonovae would be 200 Mpc which can hardly cover the horizon of GW-triggered BH-NS merger events that O3 found (as shown in Figure 10 ). However, although BH-NS merger kilonovae can hardly be detected for the present search depths, Figure 8 and Figure 9 reveal that there are great opportunities to discover on-axis afterglows associated with sGRBs or orphan afterglows if the BH components have a high-spin distribution. In order to cover the distance range for searching BH-NS kilonovae for the network of 2nd generation GW detectors as complete as possible, a search limiting magnitude m limit 23 − 24 is required as shown in Figure 10 . Present survey projects could reach this search limiting magnitude by increasing exposure times and number of simultaneous exposures. However, the GW candidates during O3 had very large localization areas with an average of thousands of square degrees (Antier et al. 2020b) . Increasing exposure times makes it hard for the present survey projects to cover such large localization areas. Therefore, during the HLVK era, we recommend that survey projects may search for jet afterglows after GW triggers with a relatively shallow search limiting magnitude. If BH-NS mergers have a high location precision, reaching a limiting magnitude to m limit 23 − 24 can be made, which gives a higher probability to discover associated kilonovae. At the late period of the HLVK era, LSST (LSST Science Collaboration et al. 2009) with its 8-meter-class telescope will come online by 2021 and begin full science operation by 2023. It can detect kilonoa emission of almost all BH-NS merger GW signals to be triggered by the 2nd generation detector network (see Figure 10 ). Due to its large field of view (9.6 deg 2 ) and deep search limiting magnitude (e.g., m r,limit = 24.7 mag with 30 s exposure time), LSST can be used to achieve follow-up detections and identifications of kilonovae with a high efficiency (Cowperthwaite et al. 2019 ). If the localization of BH-NS merger GW signals is only a few hundred square degrees, Mephisto (Er et al. 2020 ), which will see the first light by 2021 and is expected to operate in 2022, can allow a deep enough volumetric coverage of BH-NS merger events within a short period of time. Its capability of simultaneous imaging in three bands allows to obtain color information that is conducive to rapidly identifying kilonovae and afterglows. In addition, the CSST (Zhan 2011) , which also has a wide and deep survey capability, is planned to be launched in 2024. For those GW signals with sky positions to a few square degrees (Abbott et al. 2018a), CSST with a 1.1 deg 2 effective field-of-view can also efficiently detect the EM counterparts of the GW events. Its good performance in ultraviolet band can provide essential early observations that allow one to study the initial emission properties of kilonovae. Figure 12 . A schematic diagram of observed afterglow and kilonova lightcurves from the BH-NS merger exploded at low density environment, for observers with a distance of DL = 400 Mpc. The green, blue, and orange lines represent r-band afterglow emission, r-band kilonova emission, and total emission, respectively. For a given limiting magnitude of telescope, we define the timescales that the kilonova emission and total emission larger than the magnitude are ∆tr,K and ∆t r,A&K , respectively.
Even though GW-triggered follow-up observations is the best way to achieve an optimal searching for kilonovae, the time in survey projects that can be assigned to target-of-opportunity follow-up observations of GW triggers is limited. For example, only a time investment of ∼ 1% − 2% of the total survey time for LSST is designated to follow-up observations (Margutti et al. 2018b) . Therefore, it is also important to discuss how many EM counterparts of BH-NS mergers can be discovered from time-domain survey observations. In this section, we show the simulation results of the detection rates of kilonova, orphan afterglow, and total emission from BH-NS mergers for time survey observations as functions of search filter, limiting magnitude, and cadence time. In practice, if one detects a fast-evolving transient, in order to confirm whether it is a kilonova, one needs to use its recorded color evolution to perform analyses. Therefore, for specific survey projects we discuss, including PS1, ZTF, Mephisto, LSST, WFIRST, and CSST, we assume that these survey projects would search each point in two filters that can be beneficial for identifications. Values of some technical parameters for these survey telescopes are given in Table 4 . Hereafter, ρ 0 = 35 Gpc −3 yr −1 is adopted as the local event rate density of BH-NS mergers.
For the time-domain survey observations, the probability that a source is detected could be considered as the ratio of survey area within the time that the kilonova brightness is above limiting magnitude to the area of the celestial sphere. By defining ω as the areal survey rate (i.e., sky area covered per unit time) and ∆t X,K is the time during which kilonova brightness is above limiting magnitude in X band m X,limit (see Figure 12 for a schematic diagram), one can obtain the survey area as ω∆t X,K . However, high-cadence observations would restrict the survey area to ωt cad , where t cad is the cadence time, so we express the survey area as ω∆t, where ∆t is the minimum value between t cad and ∆t X,K . The detection rate for time-domain observations can be written aṡ
where ∆t = min (t cad , ∆t X,K ). If one considers both emission from kilonova and afterglow, ∆t X,K in Equation (34) can be substituted by ∆t X,A&K . We note that ω∆t is required to be smaller than 4π.
As listed in Table 6 (see Appendix C), we attach the detection rates of total emission per unit search rate (i.e., N X /ω × 10 −3 ) with different BH-NS system parameter distributions, cadence times, search filters, and search limiting Figure 13 . Detection rates of kilonovae (left panels), orphan afterglows (middle panels), and total emissions (right panels) as functions of 5σ limiting magnitude, cadence time, and search filter for time-domain survey observation. The solid, dashed, dashed-doted lines represent the detection rates with no cadence, 1 day cadence, and 3 days cadence, respectively. The panels from left to right represent the detection rates of kilonova, orphan afterglow, and total emission, respectively. We show two example BH-NS system parameters distributions, i.e., DD2−Low MNS−Low χBH for top panels, and DD2−Low MNS−High χBH for bottom panels.
magnitudes. The detection rates for time-domain survey observations tend to be higher if BH-NS systems have a less massive NS component with a stiff EoS and a high-spin primary BH. Particularly, compared with other system parameter distributions, the BH spin distribution has a more significant effect on the detection rates. Such effects of BH-NS system parameters have been discussed in detail in Section 3.1 and Section 4.2.2.
By giving two example BH-NS system parameter distributions, i.e., DD2 − Low M NS − Low χ BH and DD2 − Low M NS − High χ BH , we show the detection rates of BH-NS merger kilonova, orphan afterglow, and total emission for time-domain observations with the consideration of different limiting magnitudes, search filters and cadence times in Figure 13 . In Paper I, we revealed that BH-NS merger kilonovae could be dim in optical but bright in infrared. Due to the large amount of lanthanide-rich ejecta formed after BH-NS mergers, its longer diffusion timescale makes the infrared kilonova emission evolve slower. Therefore, the time duration when the kilonova brightness is above the limiting magnitude would be longer in the infrared band. It is obvious that searching in redder filter bands can discover more sources for time-domain survey observations. Figure 13 shows that the detection rate of kilonova in J band is ∼ 8 − 10 times more than in r band. In contrast, the detection rate ratios between J-band and r-band for orphan afterglow and total emission are close to 1.5 − 2. Figure 13 also shows that high-cadence time-domain observations, e.g., t cad = 1 day, would lose detection rates significantly. This is because the shorter the search cadence, the smaller the sky area that telescopes can cover at the same time while the detection rates would be constrained. However, it makes kilonovae or afterglows detectable in a few epochs by survey projects, which is of great advantage for identification. Using a three-day cadence with limiting magnitude of m X,limit 23 to search for orphan afterglows and total emissions would obtain the same conclusion, Figure 14 . Detection rates with different detectable epoches for kilonovae (left panels), orphan afterglows (middle panels), and total emissions (right panels). The panels from top to bottom represent EM counterparts to be detected in r band, i band, and J band, respectively. The observations can be detected in at least one, two, three, and four epoches are denoted as solid lines, dashed lines, dotted lines, and dashed-dotted lines. The red, green, and light blue lines represent the detection rates by adopting one-day cadence searching strategy, while yellow, mazarine, and purple lines represent the detection rates by adopting three-day cadence searching strategy. The example BH-NS system parameters distribution is DD2 − Low MNS − High χBH.
because afterglows, especially off-axis ones, are usually bright and last longer. If one searches kilonovae or afterglows using small telescopes with a relatively longer cadence time, e.g., t cad = 3 day, the detection rates would not decrease significantly since the cadence time is almost larger than the time when the brightness is above limiting magnitude. However, this would cause that most of sources can be only detected in one or two epochs, which makes it hard to identify them. Figure 15 . Detection rates of kilonvoa, orphan afterglow, and total emission for specific survey projects, including PS1, ZTF, Mephisto, LSST, WFIRST, and CSST. The green, yellow, and red circles represent the detection rates by adopting a one-day cadence strategy, three-day cadence strategy, and no cadence searching, respectively. The example BH-NS system parameter distributions are those of the DD2 − Low MNS − High χBH model.
In view that the kilonova and afterglow emission from BH-NS mergers whose BH component has low spin is difficult to detect for time-domain searches, we only show the detection rates of kilonova, orphan afterglow, and total emission by setting the BH-NS system parameter distributions as those of the DD2 − Low M NS − High χ BH model in Figure  14 . A three-day cadence, especially for searches in r band and i band, can barely detect kilonovae in only one epoch. However, if one adopts high-cadence searches, kilonovae are expected to be detected many times. Figure 14 shows that there are ∼ 25% detectable kilonovae that can be probed in at least 4 epochs by using a one-day cadence probed in J band. Therefore, we suggest that the time-domain survey observations for kilonovae should adopt a one-day cadence search strategy in redder bands. A similar strategy is recommended for orphan afterglow and total emission searches with a relatively shallow limiting magnitude. As also shown in Figure 14 , the difference of orphan afterglow and total emission detection rates between high-cadence and low-cadence searches is small for deep survey projects.
For time-domain survey observations, we simulate the detection rates of kilonova, orphan afterglow, and total emission based on the survey information for PS1, ZTF, Mephisto, LSST, WFIRST, and CSST, collected in Table 4 . The predicted values of detections for each survey by considering an example BH-NS system parameter distributions, i.e., DD2−Low M NS −High χ BH , are exhibited in Figure 15 . The possibility for detecting EM counterparts from BH-NS mergers is extremely low if the BH component has a low spin distribution.
Overall, for time-domain survey observations, the primary EM counterpart targets for BH-NS mergers would be on-axis jet afterglows and orphan afterglows. Compared with other survey projects, the relatively deep depth and wide field of view of LSST with a 8-meter-class telescope would detect most EM counterparts from BH-NS mergers. ZTF has the shallowest search depth, but could still discover many BH-NS merger afterglows due to its ultra wide field of view. Mephisto should have similar results to ZTF. For other survey projects, despite their deep search limiting magnitudes, the possibilities for searching EM counterparts from BH-NS mergers would be limited due to their relatively small field of view.
In this paper, based on our model proposed in the companion paper (Paper I), we have presented simulated BH-NS merger kilonova luminosity function, detectability of GWs for different generations of GW detectors, joint-search GW signals and optical EM counterparts, as well as time-domain survey blind-search detectability for BH-NS EM signals. Our simulated results provide important implications for future time-domain survey searches and targeted-ofopportunity observations of GW triggers.
Effect of BH-NS system parameter distributions -We discuss the BH-NS system parameter distributions in detail and discuss their effect on the kilonova luminosity function, GW signal detections, GW-triggered target-of-opportunity observations, and time-domain survey blind searches for BH-NS merger EM signals. Different system parameter distributions show only small effects on the GW detections, detectable distances, and chirp times of BH-NS merger events. On the contrary, different distributions can significantly influence kilonova luminosity function and EM detections. In the case of primary BH having a high-spin distribution and NS component being less massive with a stiff EoS, the NS can be disrupted by the BH for almost all cases to power a bright kilonova and afterglow (e.g., Barbieri et al. 2019 Barbieri et al. , 2020 Bhattacharya et al. 2019 ). Our simulated results show that the kilonova luminosity function and EM detectability might rely more on the BH spin distribution. For those BH-NS mergers whose primary BH has low-spin, the NS would generally plunge into the BHs, or the NS would undergo tidal disruption but only a small amount of material can remain outside the post-merger BH remnant. If the BH spin in BH-NS systems is universally low, the brightness distribution of BH-NS merger kilonovae could be dim, and the detectable event rate would be low. Therefore, in such a case, the EM counterparts from BH-NS mergers could be difficult to discover (see also Zappa et al. 2019; Drozda et al. 2020) . At present, the BHs predicted by several studies (e.g., Fuller & Ma 2019) and detected from GW observations (e.g., Abbott et al. 2019a) usually rotate very slowly, while the BHs in some observed high-mass X-ray binaries (which could be progenitors of BH-NS systems) are born with a high spin (e.g., McClintock et al. 2011; Miller et al. 2011) . In view of the little present knowledge about BH spin in BH-NS systems, the constraints by future GW detections could improve our understanding of BH-NS mergers and their EM signals.
Luminosity Function -By adopting the viewing-angle-dependent model of radioactivity-powered BH-NS merger kilonova model presented in Paper I, we for the first time present a detailed calculation of BH-NS merger kilonova luminosity function with the consideration of different BH-NS system parameter distributions. The predicted kilonova absolute magnitude at 0.5 days after a BH-NS merger is mainly distributed in the range of ∼ −10 to ∼ −15.5, corresponding to the bolometric luminosity from ∼ 3 × 10 39 erg s −1 to ∼ 5 × 10 41 erg s −1 . The simulated luminosity functions potentially contain the viewing angle distribution information of the anisotropic kilonova emission. Besides, the possible additional energy injection from the remnant BH, e.g., due to fall-back accretion (Rosswog 2007; Kyutoku et al. 2015) or the Blandford-Payne mechanism (Blandford & Payne 1982; Ma et al. 2018) , would make the kilonova luminosity function more complex. Compared with BH-NS mergers, BNS mergers could form more polar-dominated lanthanide-free neutrino-driven wind ejecta (e.g., Grossman et al. 2014; Perego et al. 2014; Kasen et al. 2015; Martin et al. 2015 ) whose mass could be dependent on the survival times of remnant NSs (e.g., Kasen et al. 2015; Fujibayashi et al. 2017; Gill et al. 2019) , so that lightcurves of BNS merger kilonovae would be diversified and highly viewing-angledependent. Furthermore, the long lasting X-ray afterglows in sGRB events (Gompertz et al. 2013 Rowlinson et al. 2014; Lü et al. 2015 ) indicated continuous energy ejection from a long-lived remnant magnetar after some BNS mergers. The spin-down of magnetars could also power X-ray tranisents (Zhang 2013; Sun et al. 2017 Sun et al. , 2019 Xue et al. 2019; Xiao et al. 2019) , give rise to a strong broadband afterglow-like emission upon interaction between the merger ejecta and the ambient medium (Gao et al. 2013a; Liu et al. 2020 ) and heat the ejecta to enhance the brightness of the kilonovae . Such magnetar-powered kilonovae could be a few tens or hundreds times brighter than the radioactive-powered kilonovae. Therefore, it can be predicted that the kilonova luminosity function would be more complex for BNS mergers. Future extensive researches of diverse kilonovae from different GW events, constraints on the event rates of compact binary mergers, NS EoS, and NS mass distribution can lead to better constraints on the kilonova luminosity function. For future detected fast-evolving transients, one can make use of of their peak luminosity and inferred event rate to judge whether these events could be kilonovae.
GW Detections -We for the first time make detailed calculations on the detection capabilities of the 2nd, 2.5th, and 3rd generation detector networks in the near future for BH-NS signals. By setting the local event rate of BH-NS mergers asρ 0 = 35 Gpc −3 yr −1 , the GW detection rates of the worst and best cases for the 2nd generation GW detector network are ∼ 6 yr −1 and ∼ 90 yr −1 , respectively. The values of the worst and best detection rates for 2.5th generation GW detector network would increase to ∼ 30 yr −1 and ∼ 700 yr −1 , respectively. For the 3rd generation detector network, LIGO Voyager can detect ∼ 6, 000 BH-NS merger events per year. A single ET and the ET+CE network are expected to detect all BH-NS merger events in the entire Universe, with the detection rates ∼ 2.6 × 10 4 yr −1 and 2.8 × 10 4 yr −1 , respectively. As the sensitivity of the GW detector increases, BH-NS events at high redshifts gradually dominate of the detected events. For ET and ET+CE, the detection rate is mainly dominated by the BH-NS mergers at z ∼ 1. We also show the chirp time of different generation detector networks. For the 2nd generation, 2.5th generation detectors, and LIGO voyager, the median chirp time is distributed in the range of ∼ 2 − 4 min, and does not exceed 10 min. Since the minimum frequency reduces from 10 Hz to 1 Hz for ET and ET+CE, the chirp time would reach the order of hours or even days with a median value close to ∼ 10 hr, which can provide an enough preparation time for EM observations. The detailed GW detection rates, median and maximum detectable distances, chirp time and total mass in geocentric frame of detectable GW signals can be consulted from Table 1, Table 2, and Table 3 . We also verify that the detection rate and detectable distance for HLV (O3) are approximately similar to those assuming only one bKAGRA running, which are basically consistent with the detection rate and distance distribution of BH-NS merger candidates observed during LVC O3. We predict the detection rates of BH-NS mergers for different generation detector networks. The predicted median detectable distances and chirp times can provide important implications for future EM follow-up observations of GW triggers.
GW-triggered target-of-opportunity observations -In view that the EM detectability is highly dependent on the BH-NS system parameter distributions, one can hardly give explicit predictions of the detection rates for GW-triggered target-of-opportunity observations. Here, we simulate the detectability of kilonovae, orphan afterglows, and total emissions with the consideration of different parameter distributions and provide some implications for the past and future EM follow-up observations. We discuss why no kilonova from BH-NS merger candidates has been detected during LVC O3. One possibility is that there is actually no EM signal produced from BH-NS mergers or the EM signals are too dim for us to detect, because the bulk of BH-NS systems possibly have low-spin primary BHs and massive NSs with a not very stiff EoS. Alternatively, it may be that the previous follow-up observation searches have been too shallow to detect BH-NS merger kilonovae. In order to cover the median BH-NS merger detectable distance for the 2nd generation GW detector networks, our simulated results show that a limiting magnitude m limit ∼ 23 − 24 mag is required. At the late period of the HLVK era, LSST (LSST Science Collaboration et al. 2009 ) with its 8-meter-class telescope will join the operation and would detect almost all BH-NS merger kilonovae with GW triggers.For those survey projects whose limiting magnitudes are m limit 23 − 24 mag, we suggest that they could search for afterglow signals after GW triggers. For the 2.5th and 3rd generation detector networks, the kilonovae are even fainter for GW-triggered BH-NS mergers. Nonetheless, follow-up observations may detect the associated sGRBs or afterglows and look for possible kilonova signatures in the lightcurves.
Time-domain blind survey observations -Since survey telescopes only invest a relatively small percentage of time in GW trigger follow-up observations, we also discuss the detectability of BH-NS merger EM counterparts for timedomain blind survey observations. In order to find more BH-NS merger kilonovae with multiple observations, a high-cadence search in redder filters is recommended. Compared with kilonova emission, on-axis jet afterglows and orphan afterglows are more likely to be discovered by time-domain blind survey observations. The detectability for partial survey projects, including ZTF, PS1, WFIRST, Mephisto, LSST, and CSST, are discussed in detail. We indicate that LSST would provide a unique combination of deep searches and wide field-of-view survey observations, which will be well suited to detect more EM counterparts from BH-NS mergers. The dimensionless redshift distribution factor f (z) describes the redshift-dependent merger event rate, i.e.,
The merger event rate of BH-NS binary systemρ(z) at redshift z can be connected to the cosmological star formation rate by accounting for the merger delay time from the formation due of the BH-NS binary system and the orbital decay time scale through gravitational wave radiation. Therefore, the merger event rateρ(z) can be expressed as a convolution of the cosmological star formation rateρ * (z) and the probability density function of delay time P (τ ), i.e.,
where τ = t(z) − t(z ) is the delay time, τ min is the minimum delay time, t(z ) is the time when the binaries are formed, and t(z) is the delay time for binaries to merge, which is related to the source redshift as Sun et al. (2015) collected three major types of merger delay time models, including the Gaussian delay model (Virgili et al. 2011) , log-normal delay model (Wanderman & Piran 2015) , and power-law delay model (Virgili et al. 2011) . They simulated the binary systems following the distribution of star formation rateρ * (z) in the redshift space, and randomly generated the merger delay timescales of all the systems based on different delay models. Then, they derived the lookback time of each system and converted the lookback time to redshift. By these steps, Sun et al. (2015) suggested some empirical formulae of the dimensionless redshift distribution f (z) for the three merger delay models. However, due to the different independent variables ofρ * (z) and P (τ ), one needs to use unified variables to calculate event rate of BH-NS mergers. Here, we use the method of convolution with redshift z by Equation (A2) to obtain the redshift distributions, which are normalized to unity at the local universe (z = 0) for the three merger delay models collected by Sun et al. (2015) . We adopt the analytical model of Yüksel et al. (2008) as the star formation rate which readsρ * (z) ∝ (1 + z) 3.4η + 1 + z 5000 −0.3η
where η = −10.
The updated empirical formulae of f (z) shown in Figure 16 are given as follows. All of these three empirical formulae can well fit the simulated data from z = 0 to z = 8. For the Gaussian delay model (Virgili et al. 2011 ) 
where η = −7.553. Figure 16 . Redshift distribution derived from the method of convolution considering three merger delay models: Gaussian delay model (red), log-normal delay model (green), and power-law delay model (blue). The corresponding colored curves represent empirical fits for each merger delay model. The gray curve is the star formation rate modeled by Yüksel et al. (2008) For the log-normal delay model (Wanderman & Piran 2015 )
where t d,LN = 2.9 Gyr and σ t,LN = 0.2, the empirical dimensionless redshift distribution f (z) is 
where η = −5.51. For the power-law delay model (Wanderman & Piran 2015) P PL (τ ) = τ −αt ,
where α t = 0.81, one has f PL (z) = (1 + z) 1.895η + 1 + z 5.722 Note-Values represent EM detection rates for GW-triggered target-of-opportunity observations are calculated by setting the BH-NS local event rate asρ0 = 35 Gpc −3 yr −1 at different GW detection eras. While the values in the brackets are the EM detection fractions in detectable GW signals. Note-Values times 10 −3 represent BH-NS total emission detection rates per unit search rate (ṄX /ω × 10 −3 ) for time-domain survey observations are calculated by setting the BH-NS local event rate asρ0 = 35 Gpc −3 yr −1 .