Davide Tornotti, Matteo Fossati, Michele Fumagalli, Davide Gerosa, Lorenzo Pizzuti, Fabrizio Arrigoni Battaia

We present a hierarchical Bayesian framework designed to infer the luminosity function of any class of objects by jointly modeling 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 could be generalized to any extensive quantity, such as mass. We validate the model using mock catalogs, recovering how deep data reaching $\gtrsim 1.5$ dex below a characteristic luminosity $\tilde{L}^\star$ are essential to reduce biases at the faint end ($\lesssim 0.1$ dex), while wide-area data help constrain the bright end. As proof of concept, we consider a combined sample of 1176 Ly$\alpha$ emitters at redshift $3 < z < 5$ drawn from several MUSE surveys, ranging from ultra-deep ($\gtrsim 90$ hr) and narrow ($\lesssim 1$ arcmin$^2$) fields to shallow ($\lesssim 5$ hr) and wide ($\gtrsim 20$ arcmin$^2$) 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 studies based on gravitational lensing that reach $\log(L/\mathrm{erg\,s^{-1}}) \approx 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 Ly$\alpha$ emitters [$\log(L/\mathrm{erg\,s^{-1}}) \gtrsim 43$], where statistical uncertainties still dominate.