Gaussian Process Foreground Subtraction and Power Spectrum Estimation for 21 cm Cosmology

One of the primary challenges in enabling the scientific potential of 21 cm intensity mapping at the Epoch of Reionization (EoR) is the separation of astrophysical foreground contamination. Recent works have claimed that Gaussian process regression (GPR) can robustly perform this separation, particularly at low Fourier $k$ wavenumbers where the signal reaches its peak signal-to-noise ratio. We revisit this topic by casting GPR foreground subtraction (GPR-FS) into the quadratic estimator formalism, thereby putting its statistical properties on stronger theoretical footing. We find that GPR-FS can distort the window functions at these low k modes, which, without proper decorrelation, make it difficult to probe the EoR power spectrum. Incidentally, we also show that GPR-FS is in fact closely related to the widely studied optimal quadratic estimator. As a case study, we look at recent power spectrum upper limits from the Low Frequency Array (LOFAR) that utilized GPR-FS. We pay close attention to their normalization scheme, showing that it is particularly sensitive to signal loss when the EoR covariance is misestimated. This implies possible ramifications for recent astrophysical interpretations of the LOFAR limits, because many of the EoR models ruled out do not fall within the bounds of the covariance models explored by LOFAR. Being more robust to this bias (although not entirely free of it), we conclude that the quadratic estimator is a more natural framework for implementing GPR-FS and computing the 21 cm power spectrum.