Joint Bayesian calibration and map-making for intensity mapping experiments

Abstract

Line-intensity mapping (LIM) is an emerging cosmological technique that traces large-scale structure through the integrated spectral-line emission of unresolved sources. Reconstructing unbiased sky maps requires careful joint treatment of instrumental calibration and map-making, a task made challenging by time-varying receiver gains, thermal drifts, and correlated $1/f$ noise intrinsic to single-dish radio telescopes. We present a Bayesian framework for joint calibration and map-making using Gibbs sampling, giving access to the full joint posterior of calibration and sky map parameters. Our data model is grounded in the radiometer equation, capturing the coupling between noise level and system temperature without assuming a fixed noise amplitude. Gain and system temperature are estimated via an iterative generalised least squares (GLS) scheme, while absolute flux calibration is achieved either with external calibrators or via known signal injections such as noise diodes. We further introduce a $1/f$ noise model that avoids spurious periodic correlations arising from the common assumption of a diagonally structured noise covariance in the frequency domain. The workflow is implemented in an efficient software package using the Levinson algorithm and a polynomial emulator to reduce computational cost. Demonstrated on simulations representative of MeerKLASS single-dish observations, the framework generalises to other single-dish surveys and to cross-correlation and interferometric data.

Figures (19)

• Figure 1: Comparison between the conventional $1/f$ noise model [Uncorrelated DFT modes (Case 1); see Appendix \ref{['Appendix: PSD']}] and our analytical approach (Case 2; see Section \ref{['sec: flicker model']}). Left: Time-domain correlation function. Middle: Power Spectral Density (PSD; defined with DFT), normalised by that of white noise. The vertical line marks the knee frequency ($1$ mHz), where the flicker noise power matches that of white noise (shown by the horizontal line). Right: Time-time covariance matrix. This figure highlights the consequences of periodicity assumptions in DFT-based $1/f$ noise models. In Case 1, the diagonal covariance in DFT space leads to artificial periodic and symmetric correlations in the time domain. In contrast, our model (Case 2) yields a covariance consistent with the underlying Gaussian statistics and statistical homogeneity of the noise. Both models share the same PSD power law and knee frequency. For Case 2, the parameters are $f_c = 1.099 \times 10^{-3}$ (rad/s) and $f_0 = 1.335 \times 10^{-5}$ (rad/s).
• Figure 2: Flowchart of the Gibbs sampling procedure. Note that there is no strict preference for the order of the different Gibbs sampling steps; in practice, however, the first step is determined by the setup of initial conditions: Parameters without specified initial values should be sampled first. Also, step (b) can be incorporated into step (d) as part of a large joint linear sampling task. In this task, the local system temperature parameters from all TOD sets are sampled alongside the sky parameters.
• Figure 3: Execution time comparison of three computational methods for evaluating the log-likelihood with real symmetric Toeplitz noise covariance matrices. We compare three implementations, listed in order of decreasing computational complexity: (1) using NumPy's slogdet and solve; (2) using NumPy's slogdet combined with SciPy's solve_toeplitz; (3) using our own comat implementation. As shown, comat outperforms the other methods in terms of execution time. This benchmark analysis was performed on a Mac Studio equipped with an Apple M3 Ultra chip.
• Figure 4: Scan pattern, covered sky and integrated beam intensity for the "$1\times$TOD" configuration.
• Figure 5: Scan pattern, covered sky and integrated beam intensity for the "$2\times$TOD" configuration.