On using the counting method to constrain the anisotropy of kilonova radiation

A large number of binary neutron star (BNS) mergers are expected to be detected by gravitational wave (GW) detectors and the electromagnetic (EM) counterparts (e.g., kilonovae) of a fraction of these mergers may be detected in multi-bands by large area survey telescopes. For a given number of BNS mergers detected by their GW signals, the expected numbers of their EM counterparts that can be detected by a survey with given selection criteria depend on the kilonova properties, including the anisotropy. In this paper, we investigate whether the anisotropy of kilonova radiation and the kilonova model can be constrained statistically by the counting method, i.e., using the numbers of BNS mergers detected via GW and multi-band EM signals. Adopting simple models for the BNS mergers, afterglows, and a simple two (blue and red)-component model for kilonovae, we generate mock samples for GW detected BNS mergers, their associated kilonovae and afterglows detected in multi-bands. By assuming some criteria for searching the EM counterparts, we simulate the observations of these EM counterparts and obtain the EM observed samples in different bands. With the numbers of BNS mergers detected by GW detectors and EM survey telescopes in different bands, we show that the anisotropy of kilonova radiation and the kilonova model can be well constrained by using the Bayesian analysis. Our results suggest that the anisotropy of kilonova radiation may be demographically and globally constrained by simply using the detection numbers of BNS mergers by GW detectors and EM survey telescopes in multi-bands.


INTRODUCTION
Mergers of binary neutron stars (BNSs) or neutron star-black hole binaries (NSBHs) are predicted to eject neutron-rich material and produce "kilonovae" phenomena (Li & Paczyński 1998;Metzger et al. 2010;Kasen et al. 2017) by the heating from the radiative decay of heavy elements formed via the rapid neutron capture process (r-process) in the ejecta (e.g., Lattimer & Schramm 1974).This prediction has confirmed by the joint detection of the gravitational wave (GW) signal from the merger of a BNS (GW170817) by the advanced Laser Interferometer Gravitational wave Observatory (LIGO) and Virgo (Abbott et al. 2017a) and its associated electromagnetic (EM) counterparts, i.e., the short gamma-ray burst (sGRB), its afterglow, and kilonova (Goldstein et al. 2017;Savchenko et al. 2017;Hallinan et al. 2017;Alexander et al. 2017;Arcavi et al. 2017;Pian et al. 2017).The discovery of GW170817 marks the beginning of the multi-messenger astronomy by using both the GW and EM signals to study the physical processes in strong gravity field and extreme physical environments (Abbott et al. 2017b).
LIGO and Virgo have detected more than 90 merger events of compact binaries (The LIGO Scientific Collaboration et al. 2021), among them only two are BNS mergers, i.e., GW170817 and GW190425 (Abbott et al. 2017a(Abbott et al. , 2020b)).It is expected that many more BNS mergers will be detected by future GW detectors, including the ★ E-mail: luyj@nao.cas.cnLIGO Voyager (Adhikari et al. 2020), Einstein Telescope (ET; Punturo et al. 2010) and Cosmic Explorer (CE; Reitze et al. 2019).In the coming fourth and fifth observation runs (O4 and O5) of LIGO/Virgo/KAGRA, it is anticipated that ∼ 10 1 − 10 3 of BNS mergers will be detected per year (Petrov et al. 2022;Weizmann Kiendrebeogo et al. 2023).In the era of the third generation GW detectors, it is expected that upto ∼ 10 3 and ∼ 10 4 − 10 5 of BNS mergers per year will be detected by LIGO Voyager and CE/ET, respectively (e.g., Belgacem et al. 2019;Singh et al. 2022;Ma et al. 2023).
The EM counterparts of some BNS mergers detected by GW detectors are also expected to be detected by the time domain large sky surveys, such as the Rubin Observatory Legacy Survey of Space and Time (hereafter Rubin; LSST Science Collaboration et al. 2009;Ivezić et al. 2019), the Zwicky Transient Facility (ZTF; Bellm et al. 2019;Graham et al. 2019), the Panoramic Survey Telescope and Rapid Response System (Pan-STARRS; Chambers et al. 2016), and the SiTian Projects (SiTian; Liu et al. 2021).Several works have investigated the detection rates of kilonovae and the afterglows of short Gamma Ray Bursts (sGRBs) associated with BNS mergers, and found that about 20, 2, 0.1, or 7 (800/50/1/2000) kilonovae (afterglows) can be discovered by Rubin, ZTF, Pan-STARRS, or SiTian (Zhu et al. 2021b;Scolnic et al. 2018a;Setzer et al. 2019a;Andreoni et al. 2022a;Yu et al. 2021;Roberts et al. 2011).Currently only GW170817 has both the GW and the EM (afterglow/kilonova) observations.There are indeed some sGRBs that have afterglow observations with kilonova signatures (Jin et al. 2016(Jin et al. , 2015(Jin et al. , 2020(Jin et al. , 2013;;Berger et al. 2013;Troja et al. 2018aTroja et al. , 2019;;Lamb et al. 2019;Kasliwal et al. 2017), of which the GW signals could not be detected because no GW detector worked at the merger time of these sGRBs.One may thus expect to detect a large number of kilonovae and sGRBs associated with BNS merger GW events in the near future.Some of the detected kilonovae may have detailed multi-epoch multiband observations as GW170817, which enables detailed studies of each one.Some of them may not have detailed observations because of their low luminosity or faint magnitude, or the long cadence time (defined as the interval between consecutive observations of the same sky area chosen from different detection strategy) (Andreoni et al. 2022a,c).The samples of BNS mergers, kilonovae, and sGRBs as a whole may also offer demographic studies on the intrinsic properties of these systems and relevant underlying physics, as well as the individual ones with detailed observations.Kilonova radiation is anisotropic because of different properties of the ejecta at different directions.The neutron rich material ejected at directions around the equator has larger opacity and thus emits relatively redder radiation (red component), while that at directions close to the polar has smaller opacity and emits relatively bluer radiation (blue component) (e.g., Metzger 2019;Villar et al. 2017).This difference is the main reason for the anisotropy of the kilonova radiation.Many efforts have been made to study such kilonova anisotropy in the literature.Metzger & Fernández (2014) first performed an axisymmetric hydrodynamic simulation to generate the light curves of two-component kilonova.Such a two-component model is required for successfully fitting the light curves of AT2017gfo (Chornock et al. 2017;Cowperthwaite et al. 2017;Coulter et al. 2017;Kasen et al. 2017;Villar et al. 2017).Furthermore, the light curves of AT2017gfo may be better fitted by a three-component model with a purple component adding into the two-component model (Villar et al. 2017;Perego et al. 2017).Recently, many useful simulations and tools were proposed and developed to study the anisotropic for kilonovae, such as the Los Alamos National Laboratory (LANL) grid of radiative transfer kilonova simulations and the time-dependent multidimensional Monte Carlo radiative transfer code POlarization Spectral Synthesis In Supernovae (POSSIS) (e.g., Wollaeger et al. 2021Wollaeger et al. , 2018;;Bulla 2019;Dietrich & Ujevic 2017;Bulla 2023;Heinzel et al. 2021).According to these studies of kilonova anisotropy, the observational properties of a kilonova observed at different viewing angle are expected to be different (e.g., Darbha & Kasen 2020;Zhao et al. 2023).Therefore, it is possible to study the anisotropy of kilonova radiation by using the statistical distribution of kilonova properties obtained from observations if a large number of kilonovae could be detected.
A bayesian analysis to infer the width of sGRB jets from BNS mergers was proposed by Farah et al. (2020) and the method is denoted as the "counting method".Based on the number of sGRB and GW detections, they showed that the posterior distribution of the effective angular width for the sGRB jets can be well constrained by assuming a binomial distribution of the likelihood.Recently, Chen & Zhang (2023) also adopted a similar method to constrain the geometry of fast radio bursts (FRB).The reason why sGRB and FRB jet can be estimated by this method is that the emission of them is highly anisotropic.Compared with the sGRB and FRB, the emission from a kilonova is relatively less anisotropic, however, the anisotropy may be still significant.It is possible to adopt the "counting method" to obtain constraint on the anisotropy of kilonova emission from a large sample of kilonovae observations.The counting method requires little observational information, and it is expected to use for The blue component ejecta is around the polar direction with an opening angle of  b , and the red component ejecta is around the merger plane with an opening angle of  r .The viewing angle  is defined as the angle between the line of sight (LOS) and the normal of the BNS orbit.a quick analysis of the large amount of data obtained by the gradually upgraded GW and optical detectors in the future.
In this paper, we investigate whether the counting method can give some constraints on the structure parameters of kilonova according to simple kilonova and afterglow models of BNS mergers.The paper is organized as follows.In Section 2, we introduce the structured twocomponent kilonova model and the afterglow model.In Section 3 and Section 4, we describe the counting method and perform simulations to demonstrate that the counting method can put substantial constraints on the aniostropy of the kilonova radiation.Conclusions and Discussions are given in Section 5.

KILONOVA AND SGRB MODELS
In this section, we introduce simple models for kilonovae and sGRB afterglows, respectively.With these models, we can estimate the light curves for any given kilonova and sGRB.

Kilonova
We adopt a simple two-component model for kilonovae, of which the diagram is schematically shown in Figure 1.In this model, a kilonovae is axi-symmetric and can be characterized by two components, i.e., a blue component around the polar angle with opening angle of  b and a red component around the equator with opening angle of  r .Note that three-component model or even more complicated models which are calculated by simulations such as POSSIS and LANL were proposed in the literature, but the fittings to the kilonova light curves by using different models only show small differences (see Villar et al. 2017;Perego et al. 2017;Darbha & Kasen 2020;Wollaeger et al. 2021Wollaeger et al. , 2018;;Bulla 2019).For the purpose of demonstration and simplification in this paper we adopt the two-component semi-analytic model.We calculate the anisotropic radiation from any kilonova according to the two-component model as those in the previous studies (Darbha & Kasen 2020;Villar et al. 2017;Kasen et al. 2017;Metzger et al. 2010;Korobkin et al. 2012).
A small amount of neutron material can be ejected out from the merger of a BNS, in which the r-process nucleosynthesis leads to the formation of radioactive heavy elements.These heavy elements are unstable and quickly decay to release a large amount of energy, leading to the kilonova phenomenon (Metzger et al. 2010).The temporal evolution of the radioactive heating rate for each component can be approximated as (Korobkin et al. 2012) where  0 = 1.3 s, and  = 0.11,  eff ej is the effective mass of the r-process nucleosynthesis ejecta, which equals the ejecta mass of the blue or red component divided by a solid angle factor of (1 − cos  b ) or sin  r .
A fraction of the total radioactive decay power  eff in in the ejecta may lead to the thermalization of the ejected material and be transferred into the kilonova radiation (Metzger et al. 2010), i.e.,  th (t) = 0.36 where , , and  are fitting constants.These constants are functions of the mass and velocity of the ejecta as seen in Barnes et al. (2016, see their Table 1).
After considering the effective radioactive power  eff in and the thermalized fraction  th , we take the form of the common analytic solution derived by Arnett (1982), and then obtain the output effective bolometric luminosity as Here  d = (2 eff ej / ej ) 1/2 , with  = 13.4 denoting a shape constant,  ej representing the mass and velocity of r-process for red or blue ejecta,  is the speed of light, and  denoting the opacity for the red or blue component.
Due to the large opacity of the ejecta, the radiation from kilonovae can be approximated as blackbody radiation.In order to simulate the multi-band light curve, each component is assumed to have a blackbody photosphere with constant radius expansion (Arnett 1982).At any given time, the temperature of each component is defined by its bolometric luminosity and radius, using the Stefan-Boltzmann law.After the ejecta cool to a critical temperature ( c ) at which the photosphere begins to recede into the ejecta, then the temperature remains constant.The temperature and radius of the photosphere for different component are given by and (5) It is worth to mention that the  phot and  phot do not change with the  eff bol scaled to the radiation inside the cone.Therefore, we can easily get the radiation from the core for each component by these parameters.

Anisotropic flux
The actual area of radiation received by the detector is equivalent to the projection of the radiation region of red and blue component on the plane which normal to the LOS.To calculate the anisotropic flux, we divide the angle of the radiation part into small solid angles as Ω.We define unit vector  as the normal to each Ω and  as the unit vector pointing to the observer.Since we can only receive the part of emission which facing to us, the condition  •  > 0 is required.The anisotropic radiation for each component can be written as where  is the frequency, () is the intensity of blackbody emission which is calculated from the  phot described in Equation ( 4),  L is the luminosity distance and the subscript "cp" denotes different component of the kilonova.We obtain the total flux  (, ) = cp  cp (, ) by summing those from different components and convert it to the AB magnitude as where () is the transmission for different filters and  is the redshift.Figure 2 shows some example light curves for kilonovae in the u-, r-, and y-bands viewing at different .As seen from this Figure, the light curves increases at the beginning, reach their peaks around 0.5-2 days after the BNS merger, and then decrease with elapsing time .After the peak, there is a bump caused by the red component, evolving slower than the blue component.Since the blue component distributes around the poles and the red component distributes near the equator, the blue component dominates the kilonova flux received by an observer viewing it at small  and the bump after the peak, due to the red component, is less obvious; the red component is prominent when  is large, and the late bump is more obvious.The light curve in each band depends on  because of the anisotropic nature of kilonovae, and this dependence in the u-band is more stronger than that in the y-band.If  b = 60 • or 30 • but  r fixed at 60 • , the peak magnitude of the light curve in the u-/y-band is expected to change by upto 0.78/0.55 or 1.33/0.90magnitudes with  changing from 0 • to 90 • .

Afterglow of GRB
GRBs can have afterglows in a wide range of bands due to the interaction of their jets with the surrounding interstellar medium (ISM), which creates internal and external counter-shocks and thus accelerates electrons to produce synchrotron radiation (Mészáros & Rees 1997).Afterglow emission depends on the structure of GRB jets and the environment of GRBs (e.g., Rossi et al. 2002;Zhang & Mészáros 2002).Lazzati et al. (2018) found that a Gaussian structure jet model can well match the observations of the GRB (GRB170817A) and afterglows associated with GW170817.Here we adopt the Gaussian structure jet model, for which the viewing-angle-dependent energy can be described as where  0 is the on-axis equivalent isotropic energy,  c represents the typical opening angle of the jets.We adopt the python package afterglowpy (Ryan et al. 2020) to calculate the light curves of afterglows.We adopt the model parameters for afterglows given by the best fit to GRB170817A observations (Troja et al. 2018b), as listed in Table 1. Figure 3 shows example light curves in the u band obtained by observers viewing the afterglows of a sGRB with redshift 0.005 at different viewing angles.As seen  Notes:  0 is the on-axis equivalent isotropic energy,  c is the opening angle of the jets,   is a truncation angle,  0 is the number density of the interstellar medium,  is the power-law index of accelerated shock,   is the magnetic field energy fraction and   is the accelerated electron energy fraction.We adopt the best fit values for these model parameters given by Troja et al. (2018b).from this Figure, the peaks of the light curves decrease sharply with increasing  especially at small viewing angles.The light curves of afterglows in other bands have similar evolution pattern and only u band is shown in Figure 3.

METHOD
In this section, we introduce the counting method and make some simulations to limit the parameters of kilonova.Simulation results for four different scenarios are also obtained in this section.
We adopt a Bayesian method to statistically constrain the parameters of kilonovae in this paper, which is similar to the one by Farah et al. (2020) for inferring the width of sGRB jets from BNS mergers.The model of the kilonova is more complicated compared with that for the sGRB afterglow.The expected anisotropy of kilonova emission could be considerable, but is not as significant as that for sGRBs.The kilonova light curves at different bands depend on the viewing angle differently, thus multi-band observations encode important information of the anisotropy.Therefore, it is possible to constrain the anisotropy of kilonova radiation by using multi-band observations.For some individual kilonovae with detailed light curve observations, one may fit the multi-band light curves by using sophisticated kilonova model to reveal the anisotropic nature of kilonova radiation, as that done for GW170817 (e.g., Perego et al. 2017;Villar et al. 2017;Drout et al. 2017;Zhu et al. 2021a;Zhao et al. 2023), with which one could obtain strong constraints on related physical parameters.
For those kilonovae without detailed or even no light curve observations, however, one may still use the information of multi-band "observed" and "unobserved" kilonovae to constrain their properties via statistical methods, such as the counting method.Several parameters that are mostly relevant to the counting method are described as follows.
• The total number of BNS merger events being detected by GWs and EM waves (either afterglows or kilonovae) in any given band with given S/N thresholds are denoted as  GW and  EM , respectively.The number of EM counterparts detected in different bands can be different, and we use the vector  EM to denote the different numbers of detected EM counterparts in different bands.In this paper, for illustration purpose, we adopt the six bands of Rubin for our calculation, i.e., u, g, r, i, z, and y bands, therefore  EM is written as  EM = { EM,u ,  EM,g ,  EM,r ,  EM,i ,  EM,z ,  EM,y }.
• The viewing angle  is defined as the angle between the LOS and the normal of the BNS orbit.
• The total mass of kilonova ejecta is denoted as  ej .We set the mass ratio (denoted as ) of the blue component to the red component at a fixed value of 0.46 (Villar et al. 2017) in our calculations.In principle, this parameter can have a wide distribution among different kilonovae, of which the effect will be discussed in Section 5.
• The fraction of the GW events that have detected EM counterparts is denoted as   =   EM / GW for one band in {, , , , , }.The goal of this paper is to investigate whether the anisotropy of kilonava radiation can be constrained by using the counting method, i.e., constraining  b ,  r , and  ej via  GW and  EM based on some mock observations.According to the Bayesian theorem, one can obtain the posterior probability distribution where the superscript '' represents any given band, and  indicate the total number of bands being considered, and it equals to six for Rubin in our calculation.We assume that the detection in each band is independent, ignoring that the observations in different bands may not be independent from each other with a certain observation strategy.With this assumption, the information from multi-bands may be cumprod in the likelihood.In the model, there is naturally a correlation among detection numbers of the events observed in different bands due to the fact that, we use different bands to observe the same mock sources which are calculated by a given set of model parameters.Therefore, the correlation between the detections in different bands is already included in the calculation of ( GW ,   EM | ej ,  b ,  r ).However, one should note that the searches of the EM counterparts in different bands may be not independent and the correlations among different bands may be determined by the searching strategy.In this paper, we do not consider this in the likelihood for simplicity.We also assume uniform priors for parameters  ej ,  b , and  r , therefore, the posterior is proportional to the likelihood.Since  can be calculated from these three parameters, the likelihood for one band can be rewritten as ( GW ,   EM |  ( ej ,  b ,  r )).We assume that the GW sources are independent from each other, then the likelihood should be a binomial distribution which can be written as where  is a bool indicator, with  EM, = 1 and  GW = 1 for the detection of EM and GW, respectively,  EM = 0 for the non-detection of EM.The detection rate is also related to redshift  and inclination angle .For GW, it has additional parameters Ψ = {, , }. and  denote the source's ecliptic colatitude and longitude,  is the polarization angle defined by Apostolatos et al. (1994) denoting the orientation of the source relative to the detector.Equation ( 11) can be then marginalized as We define a kilonova as observable, i.e., (  EM,i = 1|, ), if its peak magnitude is brighter than the threshold of the detection magnitude where Θ is a step function,   is the peak magnitude of a kilonova at a given band, which is related to the total mass of the ejecta  ej , the opening angles of the blue and red components  b ,  r . thr, is the threshold for the 5 limiting magnitude of Rubin in different band.
where (  GW |, , Ψ) = Θ((, , Ψ) −  thr ), with the signal to noise ratio (SNR) of the GW signal given by (Finn & Chernoff 1993) Here with the antenna pattern functions These angles are independent with each other, therefore we can randomly generate points for cos , /, /,  and cos  to decided the SNR for each GW.
In Equation ( 15),  0 represents the characteristic distance parameter for a GW detector and we set  0 = 1527 Mpc for ET, M  is the observed chirp mass (i.e., M  = M 0 (1 + )), Z(  max ) is a dimensionless function, reflecting the overlap between the GW signal and the effective bandwidth of a GW detector, which increases monotonically as a function of the maximum orbital frequency  max .Taylor et al. (2012) estimated that Z(  max ) is close to unity.Therefore, we set Z(  max ) = 1 for our following calculations.When LIGO and Virgo is used for calculation, a large number of GWs cannot be detected within the redshift of 0.5.This result leads to a very high value of , which affects the estimation of parameters using our method.Meanwhile, the strong detection ability of ET make the value of  become very small, because many GW can be detected at high redshift, while the corresponding optical signals are basically not detected for now.To balance this two aspects, we choose ET as the GW detector and raise the GW detection threshold to  thr = 21 in our calculation.
Substituting Equations ( 13) and ( 14) into Equation ( 12), we have According to Equation ( 19), we generate  = 10 5 mock events with different , , , ,  via Monte-Carlo (MC) method.The values of cos , cos , / and / are assumed to be randomly distributed in the range of [−1, 1].We set the maximum value for redshift as 0.5 for convenience, as kilonovae at lower redshifts are easier to be detected.For the distribution of , we consider that the merger rate of BNS evolves with redshifts.The number density per unit time for BNS mergers is given by (Chen et al. 2022;Zhao & Lu 2021) where () is the BNS merger rate density and  c is the comoving volume.We adopt the results of BNS merger rate density by the model "10.kb0.9" in Chu et al. (2022).By adopting a Lambda Cold Dark Matter (ΛCDM) universe with  0 = 70 kms −1 Mpc −1 , and Ω m = 0.3, Ω Λ = 0.7, one can obtain the comoving volume element The MC approximate shows as The masses of the ejecta for different kilonovae may be different, which depend on the total mass and mass ratio of each BNS merger.
The distribution of such masses can be given by numerical simulations (e.g., Radice et al. 2018;Dietrich & Ujevic 2017;Dietrich et al. 2017;Hotokezaka et al. 2013) or from the fittings to the observations of kilonovae (such as GW170817, see Villar et al. 2017).For simplicity, in the present paper, we assume a simple distribution for the mass of the kilonova ejecta, i.e., either an uniform distribution in the range from log  1 to log  2 or a log-normal distribution with a mean and a standard deviation of  and , respectively.Besides the distributions of the kilonova mass ejecta,  b and  r , there are some other kilonova parameters (e.g., the ejecta mass ratio , and  ej , ,  c for each component).For demonstration purposes, we fix the values of these parameters as those obtained by Villar et al. (2017) in our calculations and give a further discussion on the influence that may be caused by the distribution of these parameters in Section 5. Given the distribution of the kilonova mass ejecta,  b ,  r , and other kilonova model parameters, we can generate mock samples and obtain the   value [see Eq. ( 22) in different bands, rewritten as   ( 1 ,  2 ,  b ,  r ) or   (, ,  b ,  r )].To generate mock observations, we firstly simulate a number of BNS mergers, and calculate the expected S/N of the GW signal from each BNS merger.For any input parameter set  = { 1 ,  2 ,  b ,  r } or {, ,  b ,  r }, we can also generate the kilonova signal for each source.Combined with the detection capabilities of EM and GW detectors, we can obtain the expected  GW and  EM according to the simulations.After having the dependence of   with four parameters and the mock observation numbers, we are expected to use Equation ( 10) to re-constrain the input parameters.One may keep in mind that the expected numbers of detectable events depend on the searching strategy of the EM counterparts.Therefore, we consider four different cases as listed below.
1) "Rubin": we consider the third generation GW detector ET for the GW detection and the Rubin with the 5 limiting magnitudes given by its wide-field sky survey for searching kilonovae.The EM counterpart model considers kilonova-only.
2) "Rubin+": same as "Rubin" listed in the above item but with the 5 limiting magnitudes one magnitude deeper than those for "Rubin" to demonstrate the effect of detector sensitivity on the constraint.
3) "KN+AG&info": same as "Rubin+" but considering the influence of observations of identified afterglows in the counting.
For all the above scenarios, we further consider the influence of the observation strategy.A kilonova may be not detectable even if the peak magnitude of the light curve is above the threshold because the detection capability is influenced by different telescope detection strategies.Following Zhu et al. (2021b), for demonstration purpose, we simply adopt a factor to describe this effect for the detection probability in our calculations as min where Ω sph = 41252.96deg 2 is the total sky area, Ω cov , Ω FoV ,  ope ,  cad , Δ  ,  exp ,  exp ,  oth are the sky coverage, field of view (FoV), average operation time per day, the optimal interval between consecutive observations of the same sky area, the length of time when a kilonova is above the limiting magnitude of telescope, the condition that the detection of the kilonova by serendipitous search is exposure on at least two filters, the optimal exposure time, time spent for each visit, respectively, for a transient survey.For Rubin, we have Ω cov = 20000 deg 2 , Ω FoV = 9.6 deg 2 ,  ope = 6 hr day −1 ,  cad = 1 day, Δ  ,  exp = 2,  exp = 30 s,  oth = 15 s.For a small part of well-localized sources, Rubin will likely be used to focus on their areas to get more details of their EM signals via the Target of Opportunity (ToO) observation strategy (Andreoni et al. 2022c;Margutti et al. 2018).In this paper, we try to demonstrate a general method which only needs the information about whether the EM counterparts are detected, for which the Wide Fast Deep (WFD) survey strategy can be adopted (Chase et al. 2022;Andreoni et al. 2022a;Sagués Carracedo et al. 2021;Setzer et al. 2019b;Andreoni et al. 2022b;Scolnic et al. 2018b).We do not intend to use the detailed information that may be obtained for those EM counterparts obtained by the ToO strategy for a small number of well-localized sources.
Moreover, to illustrate the robustness of our results, we use different distribution of  ej in mock observations and reconstruction.In mock observations we still use the log-uniform distribution which is controlled by two parameters  1 and  2 (i.e., with the input parameter set  = { 1 ,  2 ,  b ,  r }) to generate samples while in our reconstruction we assume the  ej obeys log-normal distribution which controlled by  and  (i.e., consider   (, ,  b ,  r ) in posterior probability distribution).constrain the anisotropy of kilonova radiation.Therefore, we also calculate the light curves of the afterglow multi-band radiation for each BNS merger.We consider two different cases by invoking the effect of afterglow on the detection of BNS EM counterparts.The first case is for the "KN+AG&info" case listed in Section 3, in which we adopt the KN+AG model to constrain the kilonova model parameters with the counting method by assuming the afterglow signatures are identified, and the results are listed in the 5th row as "KN+AG&info" in Table 3 and shown as yellow lines in Figure 6.The model parameters  b and  r for kilonovae can still be well reconstructed, i.e.,  b = 0.61 +0.19  −0.12 rad,  r = 0.97 +0.42 −0.49rad, which means that the anisotropy of the kilonovae can be revealed via the counting method.The second case is for the "KN+AG&noinfo" case listed in Section 3, in which we only adopt the KN model to constrain the kilonova model parameters with the counting method by ignoring the afterglow signals, though the afterglow radiations are included in the mock observations.For this situation, the obtained constraints on  b and  r for the kilonova model are slightly larger than other cases, as seen from the 6th row in Table 3 and the green lines in Figure 6.To obtain robust constraint on the aniostropy of kilonova radiation and/or the kilonova model, one needs to carefully consider the effect of the afterglow radiation on the light curves of EM counterparts though it is small.
In the above calculations, we set  GW = 3000 as an example.The of their EM counterparts detected in different bands by a survey with given limiting magnitudes and searching strategy.By generating mock samples with a simple two-component kilonova model, performing mock observations, and doing MCMC fittings to reconstruct the model parameters, we demonstrate that the kilonova model can be well constrained and the anisotropy of kilonova radiation, especially represented by the opening angle of the blue components, can be revealed by future observations via the counting method.The advantage of this counting method is that it needs only a little information and can give global constraints on the model parameters for the whole population by using both the detection and non-detection of kilonovae produced from the BNS mergers observed by GW detectors.
The constraints on the opening angles of both blue and red components obtained from the counting method may depend on the prior about the ejecta mass distribution of BNS mergers.If this prior deviates from the true one, the obtained constraints may be biased away from the true one.In principle, if there are a fraction of the detected kilonovae have multi-band multi-epoch detailed observations, one may be able to extract information about the ejecta mass distribution and then use it as the prior for the counting method.We note that 1000 2000 3000 4000 5000 6000 7000 8000 The constraining capabilities  of the counting method for  b and  r as a function of the detected number of BNS merger GW events.Here  is the width of the probability distribution (from the 84% to the 16% quantiles) for the reconstructed parameter.Using the same input parameters shown in Table 3 and consider the "Rubin+" case, we carry out 30 realizations of mock observations for each given  GW to obtain the corresponding  EM .The solid lines and filled circles represent the mean values of  and its errorbars obtained by the reconstruction for all the realizations for these events.Farah et al. (2020) and Biscoveanu et al. (2020) have investigated about the prior information that can be obtained from the two observed BNS merger events, GW170817 and GW190425.We expect that one could have some prior information about the BNS merger events via both the GW and EM signals, and thus enables tight constraints on the anisotropy of kilonova radiation from the counting method.
In our calculations, we fix ,  and  ej , and assume that these parameters are consistent with the results given by the fitting results of Villar et al. (2017) for GW170817, the only BNS merger that have both GW and EM observations.In principle, these parameters can have wide distributions among different kilonovae.Ignoring the wide-range distributions of these parameters may lead to bias in the constraining of the anisotropy of kilonova radiation and the kilonova model parameters.However, if we have prior information about the distributions of these parameters, which can be included in the kilonova modelling, then we would still obtain tight constraints on the anisotropy parameters (including  b and  r ).It is expected that there will be a number of kilonovae and afterglows associated with the BNS merger GW events being detected with detailed LCs, similar as GW170817, in the near future.Specifically, in the O4 and O5 operation runs of LIGO, Virgo, and KAGRA (Abbott et al. 2020a), and the era of the third-generation GW detectors, the detection rates of kilonovae are expected to be ∼ 1 − 5, ∼ 10 − 50, and ∼ 1000 per year, respectively (e.g., Petrov et al. 2022;Weizmann Kiendrebeogo et al. 2023;Shah et al. 2023;Colombo et al. 2022;Zhu et al. 2021b).A significant fraction of these events may have multi-band light curves as those obtained for GW170817.With these detections, the detailed EM observations should enable precise extraction of the kilonova parameters like , , and  ej , and thus give their distributions.One could take these distributions into the kilonova model when using the counting method to constrain the anisotropy of kilonova radiation and the kilonova model, and thus avoid the possible bias in those constraints obtained by fixing , , and  ej , at the values as that obtained from GW170817.For simplicity, we do not incline to consider this complication but defer it to future works.
It should be emphasized that for a given set of model parameters, the lower BNS detection rate may result in fluctuations in the number of  EM across different bands due to small number statistics.However, with the high detection rate ∼ 10 4 − 10 5 of BNS mergers per year by the third generation ground-based GW detectors (ET and CE) (Belgacem et al. 2019;Singh et al. 2022;Ma et al. 2023), one would expected that the small number statistics would be overcome.
Furthermore, the kilonova model adopted in this paper is a simplified two-component model.A kilonova model consists of three components may provide a better fit to the kilonova light curves of AT2017gfo (Cowperthwaite et al. 2017).In addition, many processes that could contribute additional energy to kilonova, such as the fallback accretion (Rosswog 2007;Rossi & Begelman 2009;Zhao et al. 2023) or the magnetic winds from a long-lived magnetar remnant (Yu et al. 2013;Gao et al. 2015;Radice et al. 2018;Metzger & Piro 2014), are not considered here.These all deserve further study.Nevertheless, it is expected that the counting method can still provide a global constraint on the kilonova model and the anisotropy of the kilonova radiation, even with a more complicated kilonova model.

Figure 1 .
Figure 1.Schematic diagram of the two-component model for kilonovae.The blue component ejecta is around the polar direction with an opening angle of  b , and the red component ejecta is around the merger plane with an opening angle of  r .The viewing angle  is defined as the angle between the line of sight (LOS) and the normal of the BNS orbit.

Figure 2 .
Figure 2. Light curves of kilonovae viewing at different viewing angles in different bands.Top and bottom panels show the cases with input parameters of ( ej ,  b ,  r , ) = (0.05 ⊙ , /3, /3, 0.005) and (0.05 ⊙ , /6, /3, 0.005), respectively.Solid, dashed, and dotted dash lines represent the light curves, and the contributions to them from the red component and the blue component, respectively.Colors of the curves represent the viewing angles as indicated by the colorbar at the right side.

Figure 3 .
Figure 3.The u-band light curves from the afterglow produced by a Gaussian structured jet viewing at different angles.Solid, dashed, and dotted-dashed lines show cases with redshift  = 0.05, 0.1, and 0.5, respectively.The color of the curves represents the viewing angle of the afterglow as indicated by the right colorbar.The parameters are similar with the afterglow of GRB170817A.

Figure 4 .
Figure 4. Distributions of   ( ej ,  b ,  r ) in the  b versus  ej plane (top panels),  r versus  ej plane (middle panels), and  b versus  r plane (bottom panels) three most sensitive bands of EM detector g-, r-, and i-bands (from left to right).The values of  r for top panels are fixed at /3, the values of  b for middle panels are fixed at /6, and the values of  ej are fixed at 0.05 ⊙ for bottom panels.Different colors in this figure represent different values of   as indicated by the right colorbar, which is obtained from mock observations for BNS merger GW events with specific settings on the corresponding kilonova parameters under the "Rubin+" scenario.

Figure 6 .
Figure 6.The Probability distributions of  b and  r for five different situations which are described in Section 3. Different color histograms represent different situations.The vertical red lines are the input values to generate mock observations.
Figure7.The constraining capabilities  of the counting method for  b and  r as a function of the detected number of BNS merger GW events.Here  is the width of the probability distribution (from the 84% to the 16% quantiles) for the reconstructed parameter.Using the same input parameters shown in Table3and consider the "Rubin+" case, we carry out 30 realizations of mock observations for each given  GW to obtain the corresponding  EM .The solid lines and filled circles represent the mean values of  and its errorbars obtained by the reconstruction for all the realizations for these events.

Table 2 .
The 5 limiting magnitude ( lim ) of Rubin for six different bands.  ( ej ,  b ,  r ) in each band is the key to obtaining the likelihood function.The detection rate of EM counterparts in one band is given by