Statistical estimation of full-sky radio maps from 21cm array visibility data using Gaussian constrained realizations

An important application of next-generation wide-field radio interferometers is making high dynamic range maps of radio emission. Traditional deconvolution methods like CLEAN can give poor recovery of diffuse structure, prompting the development of wide-field alternatives like Direct Optimal Mapping and $m$-mode analysis. In this paper, we propose an alternative Bayesian method to infer the coefficients of a full-sky spherical harmonic basis for a drift-scan telescope with potentially thousands of baselines, that can precisely encode the uncertainties and correlations between the parameters used to build the recovered image. We use Gaussian constrained realizations (GCR) to efficiently draw samples of the spherical harmonic coefficients, despite the very large parameter space and extensive sky-regions of missing data. Each GCR solution provides a complete, statistically-consistent gap-free realization of a full-sky map conditioned on the available data, even when the interferometer's field of view is small. Many realizations can be generated and used for further analysis and robust propagation of statistical uncertainties. In this paper, we present the mathematical formalism of the spherical harmonic GCR-method for radio interferometers. We focus on the recovery of diffuse emission as a use case, along with validation of the method against simulations with a known diffuse emission component.