  Kilonovae produced by mergers of binary neutron stars (BNSs) are important
transient events to be detected by time domain surveys with the alerts from the
ground-based gravitational wave detectors. The observational properties of
these kilonovae depend on the physical processes involved in the merging
processes and the equation of state (EOS) of neutron stars (NSs). In this
paper, we investigate the dependence of kilonova luminosities on the parameters
of BNS mergers, and estimate the distribution functions of kilonova peak
luminosities (KLFs) at the u, g, r, i, y, and z bands as well as its dependence
on the NS EOS, by adopting a comprehensive semi-analytical model for kilonovae
(calibrated by the observations of GW170817), a population synthesis model for
the cosmic BNSs, and the ejecta properties of BNS mergers predicted by
numerical simulations. We find that the kilonova light curves depend on both
the BNS properties and the NS EOS, and the KLFs at the considered bands are
bimodal with the bright components mostly contributed by BNS mergers with total
mass $\lesssim 3.2M_\odot$/$2.8M_\odot$ and fainter components mostly
contributed by BNS mergers with total mass $\gtrsim 3.2M_\odot$/$2.8M_\odot$ by
assuming a stiff/soft (DD2/SLy) EOS. The emission of the kilonovae in the KLF
bright components is mostly due to the radiation from the wind ejecta by the
remnant discs of BNS mergers, while the emission of the kilonovae in the KLF
faint components is mostly due to the radiation from the dynamical ejecta by
the BNS mergers.
