Abstract Quantification of the significance of a candidate multi-messenger detection of cosmic events is an emerging need in the astrophysics and astronomy communities. In this paper we show that a model-independent optimal search does not exist, and we present a general Bayesian method for the optimal model-dependent search, which is scalable to any number and any kind of messengers, and applicable to any model. In the end, we demonstrate it through an example for a joint gravitational wave, high-energy neutrino, short gamma-ray burst event search; which has not been examined heretofore.