Correlation-based Beam Calibration of 21cm Intensity Mapping (2024)

Generally, a realistic beam function $B_{\nu}(\mathbf{u})$ typically has a complex functional form, showing variations in both $\mathbf{u}$ and $\nu$. As demonstrated in the previous section, the efficacy of foreground removal depends on the beam’s relative response. This prompts the decomposition of the beam into two distinct components:

where $\bar{B}$ denotes the frequency-independent part, aligned with the beam at the reference frequency, i.e., $\bar{B}=B_{\nu_{\rm ref}}$. On the other hand, the frequency-dependent relative beam $B^{\rm rel}$ is defined as $B^{\rm rel}_{\nu}=B_{\nu}/\bar{B}$. In the following, the primary objective of our method is to determine the relative beam $B^{\rm rel}_{\nu}$ with existing data.

To achieve this, we aim to develop a filtering function designed to uniformize the beam response across all frequencies

Consequently, applying this filter $f_{\nu_{\rm ref}\nu}(\mathbf{u})$ to the observed map $\widetilde{T}_{\nu}(\mathbf{u})$

results in a map $\bar{B}(\mathbf{u})\,T_{\nu}(\mathbf{u})$ that is convolved with the same reference beam $B_{\nu_{\rm ref}}$ at every frequencies.

Therefore, the core idea of our algorithm revolves around obtaining the relative beam $B^{\rm rel}_{\nu}$ through internal calibration.First, let us re-express the temperature map at a given frequency $\nu$ and $\mathbf{u}$ cell as the product of the frequency-dependent relative beam $B^{\rm rel}_{\nu}(\mathbf{u})$ and the remaining component

Here we have ignored the noise term $n_{\nu}(\mathbf{u})$ and absorbed the frequency-independent beam in the temperature map and define

To proceed, we aim to construct a filter denoted as $F_{\nu_{\rm ref}\nu}(\mathbf{u})$. Applying this filter to observed temperature $\widetilde{T}_{\nu}(\mathbf{u})$ at each frequency channel effectively inverts the relative beam, yielding an estimator for $\bar{T}_{\nu}$

Consequently, our goal is to design the filter $F_{\nu_{\rm ref}\nu}(\mathbf{u})$ so that it functions similarly to the relative filter $f_{\nu_{\rm ref}\nu}(\mathbf{u})$ as described in Eq.(10).This could be achieved via the minimal variance estimator, derived by examining the variance

where $\langle\cdots\rangle$ presents the averaging within specific grid cell in $\mathbf{u}$. Minimizing $r^{2}(\mathbf{u})$ leads to the well-known Wiener filter (WF) ¹⁰¹⁰10hence this explains the meaning of the superscript ${\rm WF}$ in Eq. (14)

Here, we have introduced the correlation function $C^{\widetilde{T}}_{\nu_{1}\nu_{2}}$, defined as

Substituting Eq.(3) of the observed temperature map $\widetilde{T}$ and temporarily neglecting the noise matrix $N_{\nu_{1}\nu_{2}}=\langle n_{\nu_{1}}n_{\nu_{2}}\rangle$, we obtain

It is apparent that

where the additional intrinsic temperature factor $F^{T}_{\nu_{\rm ref}\nu}(\mathbf{u})$ is defined as

It is apparent that if we ignore the frequency variation of the temperature map, causing the factor $F^{T}_{\nu_{\rm ref}\nu}(\mathbf{u})$ to become unity, the constructed filter simply represents the average value of the relative beam across each given $\mathbf{u}$ cell. Therefore, $F_{\nu_{\rm ref}\nu}(\mathbf{u})$ is indeed the filter that meets our requirements.However, in reality, the image of the foreground varies with frequency, thus introducing an additional undetermined function on top of the relative beam. Various strategies can be employed to mitigate its influence, including using a foreground model as prior information. As we will discuss in detail in Section 5, the impact of this term is not significant.

Furthermore, observational data contains instrumental noise (Equation 3), making the filter described in Equation (19) usually sub-optimal for analyzing real data.Various methods can be adopted to effectively isolate the beam variations in such situation. We will explore these approaches in more detail and discuss their applications to real data analysis in a separate paper.In this work, however, it is worth mentioning that a straightforward solution involves considering the full correlation matrix $C^{\widetilde{T}}_{\nu_{i}\nu_{j}}(\mathbf{u})$ across all frequencies

Here the observed temperature map $\widetilde{T}_{\nu_{i}}(\mathbf{u})$ is a combination of sky $B^{\rm rel}_{\nu}(\mathbf{u})\bar{T}_{\nu}(\mathbf{u})$ and noise $n_{\nu}(\mathbf{u})$. Assuming that the amplitude of the sky is much larger than instrumental noise, we can then perform the SVD of matrix $\bm{C}^{\widetilde{T}}(\mathbf{u})$

where the diagonal matrix ${\bm{\Sigma}}={\rm diag}(\lambda_{i})$ store all eigen-values in decending order, and $\bm{U}$ (and $\bm{V}$) lists the corresponding eigen-vectors.As long as the sky, primarily comprised of the foreground, dominates the matrix $C^{\widetilde{T}}_{\nu_{i}\nu_{j}}(\mathbf{u})$, one could approximate the observed temperature with the value of the first eigen-mode.Therefore, in a similar approach to constructing the filter in Eq.(19), we can define the filter using only the first eigenvector $U^{(1)}_{\nu}(\mathbf{u})$, which can be expressed as:

The only underlying assumption here is the dominance of the sky temperature relative to the instrumental noise, which is achievable with sufficient integration time.Consequently, Eq. (23), representing the ratio between the first eigen-modes at two different frequencies, characterizes the desired filter with reduced noise influence.Therefore, this approach has the advantage of requiring minimal information about the noise covariance matrix.

For simplicity, the sky model for our simulation is generated using the publicly available $\mathtt{CRIME}$ code Alonso etal. (2014) ¹¹¹¹11http://intensitymapping.physics.ox.ac.uk/CRIME.html,which includes five distinct components:(a) HI signal,(b) Galactic synchrotron emission,(c) Galactic free-free emission,(d) Extra-galactic free-free emission,(e) Extra-galactic point sources emission.As demonstrated in Alonso etal. (2014), the HI signal here is derived from the lognormal model (Coles & Jones, 1991) of matter fluctuation $\delta(\mathbf{x})$, which is then transformed into a map of 21cm brightness temperature fluctuations

where $\delta_{\rm HI}$ represents the HI overdensity, smoothed over the volume element defined by $\delta\nu$ and $\delta\Omega$. Here the HI fraction $x_{\rm HI}(z)=0.008(1+z)$ is assumed to be fixed. Although a more complex modeling approach might be beneficial, it is not critical for the purposes of our study. In this paper, we consider the frequency range between $700\!-\!900\,\rm MHz$ corresponding to the redshift range $0.58\!<\!z\!<\!1.03$ in our choice of cosmology.

In terms of foreground modeling, $\mathtt{CRIME}$ models the large-scale Galactic synchrotron radiation (b) by extrapolating the $408\,\rm MHz$ Haslam map(Haslam etal., 1982). This extrapolation utilizes a direction-dependent spectral index derived from the Planck Sky model(Delabrouille etal., 2013), formulated as:

where $\beta(\mathbf{\hat{n}})$ is the synchrotron spectral index, varying with the line-of-sight (LoS) direction $\mathbf{\hat{n}}$. To capture structures on scales smaller than Haslam resolution, $\mathtt{CRIME}$ incorporates a constrained Gaussian realization with parameters from (Santos etal., 2005; Alonso etal., 2014). Similarly, for the Galactic free-free (c) and Extra-galactic point sources (e) emissions, $\mathtt{CRIME}$ also uses Gaussian realizations based on specific angular power spectra (Alonso etal., 2014).

In this subsection, we introduce the noisy beam model utilized in this study.For simplicity, we consider a single dish telescope with a Gaussian model as the unperturbed beam

with $\sigma=\theta_{\rm\tiny FWHM}/(2\sqrt{2\ln 2})$. The full-width at half-maximum (FWHM) $\theta_{\rm\tiny FWHM}$ is determined by

where $D$ is the diameter of the telescope dish. In many instances, this model offers a sufficiently accurate approximation. For example, as demonstrated in Liao etal. (2016), it provides a robust first-order approximation for both the intensity and polarization of the Green Bank Telescope.

As demonstrated in section 3, our method is implemented in Fourier space. In this domain, the Gaussian beam is represented as

In the following, we focus on a single dish telescope resembling the Green Bank Telescope (GBT) with a $100$-meter diameter, corresponding to the angular resolution $\theta_{\rm\tiny FWHM}$ that varies from $14$ to $18~{}{\rm arcmin}$ across the frequency range $\nu=700\!-\!900~{}{\rm MHz}$.

In this paper we consider two simple models of beam fluctuation. The first model (Model-I) treats the beam fluctuation as a three-dimensional random field in the $\mathbf{u}-\nu$ space,

where $\delta B_{\nu}(\mathbf{u})$ is the beam noise, and $\epsilon_{\nu}(\mathbf{u})$ is a zero mean random field in $\mathbf{u}-\nu$ space, defined as

Here, $\varepsilon$ is the amplitude coefficient of the noise, and ${\rm rand}$ is a uniformly distributed random variable within the range of $-1$ to $1$. Given the entirely random nature of the noises across frequencies, this model represents a pessimistic scenario regarding potential beam noise.

In our second model of beam noise (Model-II), we consider a more realistic scenario where the noise originates within the telescope’s spatial configuration. This noise manifests as a two-dimensional random field

where $\mathbf{d}$ represents the distance within the telescope’s frame. This noise field is subsequently rescaled and projected onto the $\mathbf{u}$ space at various frequencies using the relation $\mathbf{u}={\mathbf{d}}/\lambda$

After convolving our simulated datacube with the noisy beam model, we derive the mock observation data $\widetilde{T}_{\nu}(\mathbf{u})$ in Fourier space

as well as the corresponding map in real space $\widetilde{T}_{\nu}(\mathbf{\hat{n}})$.

After estimating and applying the CBC filter to our mock data, we now proceed with foreground removal and statistical analysis. Specifically, we adopt the blind Principal Component Analysis (PCA) cleaning technique following the methods detailed in Masui etal. (2013) and Switzer etal. (2013).We begin by representing the temperature map as a matrix $\bm{T}$, with dimensions $N_{\nu}\times N_{\theta}$. The PCA modes are then obtained through the decomposition of $\bm{C}=\bm{TT}^{\rm T}/N_{\theta}=\bm{U\Lambda U}^{\rm T}$, where $\bm{\Lambda}$ is a diagonal matrix arranged in descending order. The ‘cleaned’ map, obtained by subtracting the first $\rm m$ PCA modes, is calculated using Eq. (5).We continue by calculating the auto-power spectrum of the cleaned map

Additionally, we also compute the cross-power spectrum with the injected true HI signal. Specifically, our interest lies in the cross-correlation ratio between the two fields, defined as

Here $T^{\rm HI}$ represents the injected HI signal and $P^{\rm cl\times HI}(k)$ is the cross-power spectrum between the signal and the cleaned map.