2025A&A...704A.201T

We present a hierarchical Bayesian framework designed to infer the luminosity function of any class of object by jointly modelling data from multiple surveys with varying depth, completeness, and sky coverage. Our method explicitly accounts for selection effects and measurement uncertainties (e.g. in luminosity) and can be generalized to any extensive quantity, such as mass. We validated the model using mock catalogues; from this we determined that deep data reaching ≳1.5 dex below a characteristic luminosity ($ \tilde{L}^\star $) are essential to reducing biases at the faint end (≲0.1 dex) and that wide-area data help constrain the bright end. As a proof of concept, we considered a combined sample of 1176 Lyman α emitters at redshift 3 < z < 5 drawn from several MUSE surveys, ranging from ultra-deep (≳90 h) and narrow (≲1 arcmin²) fields to shallow (≲5 h) and wide (≳20 arcmin²) fields. With this complete sample, we constrain the luminosity function parameters $ \log(\Phi^\star/\mathrm{Mpc^{-3}}) = -2.86^{+0.15}_{-0.17} $, $ \log(L^\star/\mathrm{erg\,s^{-1}}) = 42.72^{+0.10}_{-0.09} $, and $ \alpha = -1.81^{+0.09}_{-0.09} $, where the uncertainties represent the 90% credible intervals. These values are in agreement with the results of studies based on gravitational lensing that reach log(L/erg s^–1)≈41, although differences in the faint-end slope underscore how systematic errors are starting to dominate. In contrast, wide-area surveys represent the natural extension needed to constrain the brightest Lyman α emitters [log(L/erg s^–1)≳43], where statistical uncertainties still dominate.