A high-resolution self-consistent whole sky foreground model

TL;DR

This paper introduces a high-resolution, self-consistent sky foreground model combining diffuse emissions and point sources, crucial for simulating and analyzing 21cm cosmological experiments amidst strong foreground contamination.

Contribution

It presents a novel method to construct a high-resolution, self-consistent sky model at low frequencies, integrating diffuse foregrounds and point sources using actual survey data and statistical methods.

Findings

The model effectively combines diffuse and point source emissions.

Adding point sources significantly impacts the angular power spectrum.

The generated sky maps and catalogs are publicly available.

Abstract

The neutral hydrogen 21cm line is potentially a very powerful probe of the observable universe, and a number of on-going experiments are trying to detect it at cosmological distances. However, the presence of strong foreground radiations such as the galactic synchrotron radiation, galactic free-free emission and extragalactic radio sources make it a very challenging task. For the design of 21cm experiments and analysis of their data, simulation is an essential tool, and good sky foreground model is needed. With existing data the whole sky maps are available only in low angular resolutions or for limited patches of sky, which is inadequate in the simulation of these new 21cm experiments. In this paper, we present the method of constructing a high resolution self-consistent sky model at low frequencies, which incorporates both diffuse foreground and point sources. Our diffuse map is constructed by generating physical foreground components including the galactic synchrotron emission and galactic free-free emission. The point source sample is generated using the actual data from the NRAO VLA Sky Survey (NVSS) and the Sydney University Molonglo Sky Survey (SUMSS) where they are available and complete in flux limit, and mock point sources according to statistical distributions. The entire model is made self-consistent by removing the integrated flux of the point sources from the diffuse map so that this part of radiation is not double counted. We show that with the point sources added, a significant angular power is introduced in the mock sky map, which may be important for foreground subtraction simulations. Our sky maps and point source catalogues are available to download.

Tables

Table 1. Table 1: Sky maps which are used to model the diffuse foreground.

Frequency	FWHM	Reference
10 MHz	${2.6}^{\circ }\times {1.9}^{\circ }$	[34]
22 MHz	${1.1}^{\circ }\times {1.7}^{\circ }$	[35]
45 MHz	$5^{\circ }$	[36]
85 MHz	${3.8}^{\circ }\times {3.5}^{\circ }$	[37]
150 MHz	${2.2}^{\circ }$	[37]
408 MHz	${0.85}^{\circ }$	[22]
1420 MHz	${0.6}^{\circ }$	[38]
2300 MHz	${2.3}^{\circ }\times {1.9}^{\circ }$	[39]
23, 33, 41, 61, 94 GHz (WMAP)
	$\sim {0.5}^{\circ }$	[40]

Table 2. Table 2: The fitted β i subscript 𝛽 𝑖 \beta_{i} parameters.

parameter	value and error
${\beta }_0$	$2.55\pm 0.035$
${\beta }_1$	$0.06\pm 0.002$
${\beta }_2$	$\left(3.33\pm 0.069\right)\times {10}^{-2}$
${\beta }_3$	$\left(2.57\pm 0.093\right)\times {10}^{-3}$
${\beta }_4$	$\left(-1.7\pm 0.028\right)\times {10}^{-3}$
${\beta }_5$	$\left(-9.14\pm 0.162\right)\times {10}^{-5}$
${\beta }_6$	$\left(2.86\pm 0.038\right)\times {10}^{-5}$

Equations

$$I_{\Theta }(\vec{n},\nu )=I_{\Theta }^f(\vec{n},\nu )+\sum _{p\in \text{M}(\Theta )}{I}_{\Theta }^p(\vec{n},\nu )$$

$$I_{\Theta }(\vec{n},\nu )=I_{{\Theta }_0}^f(\vec{n},\nu )+\sum _{p\in \text{M}(\Theta )}{I}_{\Theta }^p(\vec{n},\nu )$$

$$n(S)=k{S}^{-\gamma }$$

$$S_c={\left[\frac{Q^{2}k\pi {\theta }_0^2}{(3-\gamma )(\gamma -1)4\ln 2}\right]}^{\frac{1}{\gamma -1}}$$

$$T_{\mathrm{diffuse}}(\nu )=T_{\mathrm{diffuse}}({\nu }_{\ast }){\left(\frac{\nu }{{\nu }_{\ast }}\right)}^{-\beta }$$

$$\ln \left[\frac{T_{\mathrm{diffuse}}(\vec{n},\nu )}{T_{\mathrm{diffuse}}(\vec{n},{\nu }_{\ast })}\right]={\beta }_0(\vec{n})\left[\ln \frac{\nu }{{\nu }_{\ast }}\right]+{\beta }_1(\vec{n}){\left[\ln \frac{\nu }{{\nu }_{\ast }}\right]}^2+\cdots$$

$$y=\text{A}x+n$$

$$\text{A}=\left[\begin{array}{cccc}\ln \frac{{\nu }_1}{{\nu }_{\ast }}& {\left[\ln \frac{{\nu }_1}{{\nu }_{\ast }}\right]}^2& \cdots & {\left[\ln \frac{{\nu }_1}{{\nu }_{\ast }}\right]}^{{\mathrm{N}}_{\beta }}\\ \ln \frac{{\nu }_2}{{\nu }_{\ast }}& {\left[\ln \frac{{\nu }_2}{{\nu }_{\ast }}\right]}^2& \cdots & {\left[\ln \frac{{\nu }_2}{{\nu }_{\ast }}\right]}^{{\mathrm{N}}_{\beta }}\\ ⋮& & \ddots & \\ \ln \frac{{\nu }_{{\mathrm{N}}_{\nu }}}{{\nu }_{\ast }}& {\left[\ln \frac{{\nu }_{{\mathrm{N}}_{\nu }}}{{\nu }_{\ast }}\right]}^2& \cdots & {\left[\ln \frac{{\nu }_{{\mathrm{N}}_{\nu }}}{{\nu }_{\ast }}\right]}^{{\mathrm{N}}_{\beta }}\end{array}\right]$$

$$\hat{x}=[{\text{A}}^{\mathrm{T}}{\text{N}}^{-1}\text{A}{]}^{-1}{\text{A}}^{\mathrm{T}}{\text{N}}^{-1}y$$

$$n(S)\approx 1300S^{-1.77}$$

$$\displaystyle P(>\theta)=\left\{\begin{array}[]{ll}\left(\frac{\theta}{\theta_{0}}\right)^{-\gamma}\,,&\quad\theta\geq\theta_{0}\\ \quad 1\,,&\quad\theta<\theta_{0}\end{array}\right.$$

$$B_{\mathrm{circ}}(\theta ,\phi )=\exp \left\{-4\ln 2\frac{{\theta }^2}{{\theta }_0^2}\right\},{\theta }_0=1.029\frac{\lambda }{D}$$

$$B_{\mathrm{ell}}(\theta ,\phi )=\exp \left\{-4\ln 2{\theta }^2\left[{\left(\frac{\cos \phi }{{\theta }_1}\right)}^2+{\left(\frac{\sin \phi }{{\theta }_2}\right)}^2\right]\right\}$$

$$R(x)={\int }_{{\Omega }_b}\frac{n\left(\frac{x}{B(\theta ,\phi )}\right)}{B(\theta ,\phi )}\mathrm{d}\Omega$$

$$R(x)=k{x}^{-\gamma }\frac{\pi {\theta }_0^2}{(\gamma -1)4\ln 2}$$

$${\sigma }_c^2={\int }_0^{S_c}{x}^{2}R(x)\mathrm{d}x$$

$$S_c={\left[\frac{Q^{2}k\pi {\theta }_0^2}{(3-\gamma )(\gamma -1)4\ln 2}\right]}^{\frac{1}{\gamma -1}}$$

$$S^{5/2}n(S)=k_n{S}^{5/2-{\gamma }_{{}_n}},n=1,2,3,4,5,\cdots$$

$$\displaystyle\begin{array}[]{lrl}\gamma_{1}=1.57,\quad k_{1}=29500&\quad{\rm for}&\quad\qquad\quad\,\,S\leq 10^{-4}\,\,{\rm Jy}\\ \gamma_{2}=2.23,\quad k_{2}=60&\quad{\rm for}&\quad 10^{-4}<S\leq 10^{-3}\,\,{\rm Jy}\\ \gamma_{3}=1.77,\quad k_{3}=1300&\quad{\rm for}&\quad 10^{-3}<S\leq 10^{-1}\,\,{\rm Jy}\\ \gamma_{4}=2.20,\quad k_{4}=440&\quad{\rm for}&\quad 10^{-1}<S\leq 1\,\,\,\,\quad{\rm Jy}\\ \gamma_{5}=2.77,\quad k_{5}=390&\quad{\rm for}&\quad\qquad\quad\,\,S>1\,\,\,\,\quad{\rm Jy}\end{array}$$

$$\displaystyle S_{c}=\left\{\begin{array}[]{rl}0.086\,{\rm mJy}\cdot Q^{3.51}\left(\frac{\theta_{0}}{\rm arcmin}\right)^{3.51}\left(\frac{\nu}{{\rm GHz}}\right)^{-2\alpha}\,,&\quad\qquad\quad\,\,\theta_{0}\leq 0.25^{\prime}\\ 0.12\,{\rm mJy}\cdot Q^{1.53}\left(\frac{\theta_{0}}{\rm arcmin}\right)^{1.53}\left(\frac{\nu}{{\rm GHz}}\right)^{-2\alpha}\,,&\quad 0.25^{\prime}<\theta_{0}\leq 1.13^{\prime}\\ 0.017\,{\rm mJy}\cdot Q^{2.56}\left(\frac{\theta_{0}}{\rm arcmin}\right)^{2.56}\left(\frac{\nu}{{\rm GHz}}\right)^{-2\alpha}\,,&\quad 1.13^{\prime}<\theta_{0}\leq 7.18^{\prime}\\ 0.51\,{\rm mJy}\cdot Q^{1.63}\left(\frac{\theta_{0}}{\rm arcmin}\right)^{1.63}\left(\frac{\nu}{{\rm GHz}}\right)^{-2\alpha}\,,&\quad 7.18^{\prime}<\theta_{0}\leq 30^{\prime}\\ 5\,{\rm mJy}\cdot Q^{1.18}\left(\frac{\theta_{0}}{\rm arcmin}\right)^{1.18}\left(\frac{\nu}{{\rm GHz}}\right)^{-2\alpha}\,,&\quad\qquad\quad\,\,\theta_{0}>30^{\prime}\end{array}\right.$$

$${\rho }_m(L)=\sum _{i=1}^N{\left(\frac{1}{V_m}\right)}_i$$

$$\sigma =\sqrt{\sum _{i=1}^N{\left(\frac{1}{V_m}\right)}_i^2}$$

$${\rho }_m(L)=C{\left[\frac{L}{L_{\ast }}\right]}^{1-a}\exp \left\{-\frac{1}{2}{\left[\frac{{\log }_{10}\left(1+\frac{L}{L_{\ast }}\right)}{\sigma }\right]}^2\right\}$$

$$C$$

$$L_{\ast }$$

$$a$$

$$\sigma$$

$${\log }_{10}({\rho }_m)=Y-\frac{3}{2}{\log }_{10}(L)-\sqrt{B^2+{\left[\frac{{\log }_{10}(L)-X}{W}\right]}^2}$$

$$Y$$

$$B$$

$$X$$

$$W$$

$$S=\frac{L}{4\pi D_L^2(1+z{)}^{1+\alpha }}=\frac{L}{A(1+z{)}^{1+\alpha }}$$

Full text

\PACS

95.80.+p, 95.85.Bh, 98.70.Dk, 95.75.Mn, 95.75.Pq

1 Introduction

Tomographic observation of the redshifted 21cm signal from neutral hydrogen will allow us to study the Epoch of Reionization (EoR) [1, 2] and the Cosmic Dawn [3], to trace the evolution of the large scale structure of our universe and constrain the cosmological parameters including the equation of state of dark energy [4, 5], and even to probe the cosmic dark age [6]. A number of telescopes have been developed to detect the cosmological 21cm signal, including the EoR experiments such as LOFAR [7], MWA [8], PAPER [9], and HERA [10]; and the dark energy experiments such as Tianlai [11] and CHIME [12]; and the next generation telescopes such as SKA [13] which are expected to provide considerably more accurate power spectrum and directly map the large scale structures.

However, detecting the cosmological 21cm signal is very difficult due to the presence of foregrounds which are several orders of magnitude stronger. At low frequencies the main astrophysical foregrounds include the galactic synchrotron emission which arises from the relativistic electrons moving in the galactic magnetic field, the galactic free-free emission which is produced by the scattering of free electrons with diffuse warm ionized gas, strong radio sources such as supernova remnants (SNRs), and extragalactic radio sources such as radio-loud galaxies, quasars and BL Lac objects [14, 15]. In order to study the foreground removal methods and evaluate its impacts on the 21cm experiments by numerical simulation, mock sky images with foregrounds are needed.

Various studies on the foregrounds have been carried out for the 21cm experiments [16, 17, 18, 19, 20, 21]. To investigate how the different foreground subtraction techniques work, or to make forecast of the experiments, simulations with foreground sky models are needed. The Global Sky Model (GSM; [17]) is widely used for this purpose, with a recent updated version [19]. We denote the original and updated GSM models as GSM2008 and GSM2016 respectively. The foreground map in the GSM model is based on the method of Principal Component Analysis (PCA). GSM uses several whole sky or large area sky maps, especially the Haslam map [22] at 408 MHz as input data, to extract the principal components in the pixel-frequency space, to generate maps at other frequencies. GSM2008 applies the PCA once on the whole pixel-frequency data set, while GSM2016 applies the PCA iteratively on the outputs of each iterations. However, owing to the limitation of the available whole sky surveys (e.g. the Haslam map has an angular resolution of $0.85^{\circ}$), the angular resolutions of these models are limited at low frequencies.

A high resolution whole sky model is provided in the 21cm simulation software package CORA [20, 21]. It uses the 2003 version of the Haslam map 111https://lambda.gsfc.nasa.gov/product/foreground/haslam_408.cfm as the base map, with the spectrum index given by GSM, plus a randomly generated variation in the spectrum index. To include structures on the small scales, CORA added bright point sources drawn from the matched ones in the VLSS [23] and NVSS [24] surveys, plus a mock sample of randomly generated point sources according to the distribution given in Ref.[25], and finally, based on Ref.[26], very faint point sources are taken into account as a whole in the small scale angular power spectrum. There are also other high resolution models, e.g. the Tiered Radio Extragalactic Continuum Simulation (T-RECS; [27]), but only a limited patches of the sky are covered.

In this study, we construct a high-resolution whole sky foreground model by modeling the sky as the sum of a diffuse component and a point source component. For the diffuse component, our model differs from the GSM model in that instead of using the purely mathematical approach of PCA in GSM (2008 and 2016), we fit the spectral index, varying both with frequency and angular direction using the multi-wavelength data. Such model is less general than the PCA model, but may have some advantages in certain applications, e.g. employing data from other wave bands to improve the foreground model, or to make more extreme extrapolation. We also use the updated Haslam map which is source-subtracted and de-striped [28]222http://www.jb.man.ac.uk/research/cosmos/haslam_map as our basis. We achieve a higher angular resolution by adding point sources to the low resolution maps. This procedure is possible because point sources, i.e. sources unresolved up to the angular resolution of the survey, have been detected in surveys of higher angular resolutions, further, their statistical distribution is known or can be inferred. However, these point sources are not detected individually in surveys with low angular resolutions, but are merged into the diffuse foreground. Developed independently, our procedure differs from that of CORA in detail, for example in the choice of statistical model, and the placement of the points–CORA assumed a uniformly random distribution of points, whereas we consider angular clustering of the radio sources and use Rayleigh-Lévy random walk to place the mock point sources. However, we found that our results are similar to that obtained using CORA.

The remainder of this paper is organized as follows. In Section 2, we outline the construction of our foreground model. In Sections 3 and 4, we present our method of constructing the diffuse map and point source contribution respectively. Finally, the resulting total sky map is presented in Section 5 and discussed.

2 Foreground Model

In an astronomical observation, a source of radiation can be classified either as an extended source or a point source, where the latter is unresolved up to the limit of angular resolution. This classification depends on the angular resolution of the observation: at higher resolution some point sources in low resolution surveys will turn out to be extended ones.

In our sky model, we consider two types of contributions: a whole sky diffuse foreground which include all extended sources, and the point sources. The sky intensity for an angular resolution $\Theta$ is then given by

[TABLE]

where $\vec{n}$ denotes the direction of the sky, $\nu$ denotes the frequency, $I^{f}_{\Theta}$ denotes the whole sky diffuse foreground, $\textbf{M}_{\Theta}$ denotes the set of point sources at resolution $\Theta$, and $I^{p}_{\Theta}$ denotes the radiation from point source $p$ convolved with a beam of angular resolution $\Theta$. The observed diffuse foreground map and catalog of point sources are available only at some frequencies, and in most cases only a part of the whole sky. However, we may construct statistical models with parameters fitted with observations, then use the resulting model to generate mock diffuse foreground maps and point source catalogs.

However, in practice, diffuse maps are usually obtained not at the simulated resolution $\Theta$, but at a lower angular resolution $\Theta_{0}$, as the future experiments being studied are usually designed to have higher resolution than the existing ones. To address this issue, we shall adopt a simple approach: the sky is modeled as a diffuse foreground derived from the fixed angular resolution, plus point sources, i.e.

[TABLE]

For example, for the low frequency sky map we use the Haslam map as a basis, which has $\Theta_{0}=0.85^{\circ}$. The construction of our diffuse foreground model is described in section 3.

Compared with Eq. (1), the model given by Eq. (2) essentially assumes that the small scale angular power originates from point sources, whose statistical distribution is assumed to be known, and that there is no significant contribution of the extended sources at the scales between $\Theta$ and $\Theta_{0}$. Of course, if the presence of certain extended sources are known, it may also be possible to generate mock models of these. In the present paper we simply ignore such extended sources. Even so, many tests of foreground subtraction techniques or forecast on sky surveys can be conducted using such a model.

The radiation of a source can be detected above the noise if its flux is higher than the sensitivity limit of the observation, but whether this radiation can be detected as from a distinct source also depends on the number density of sources. Only above a certain threshold, i.e. the confusion limit [29, 30, 31], can the source be distinguished from other sources, below the threshold there are likely multiple faint sources of comparable flux within the same beam that they merge into each other and become indistinguishable from the diffuse foreground. Thus, in the simulation, we consider point sources above the confusion limit, while taking the ones below the threshold as part of the diffuse foreground.

Assuming that the sources are randomly distributed over the sky with a differential source count approximated as a power-law,

[TABLE]

For a circular Gaussian beam, the confusion limit flux is [30, 31]

[TABLE]

where $\theta_{0}=1.029\lambda/D$, $Q$ is a factor of 3–5 [32]. Note that the confusion limit depends on the angular resolution of the observation and the distribution of sources, but not on the sensitivity of the telescope. A detailed computation of the confusion limit is given in Appendix A.

As noted above, we may generate the radiation from the set of point sources M, which includes all sources with flux above the confusion limit $S_{c}$. However, in practice such data may not be available. Available surveys may only cover a small part of the sky, and even in the covered regions, the survey may not be complete up to the confusion limit. However, if a statistical distribution of the point sources can be inferred from the observations, we may be able to generate a mock sample of point sources according to the statistical distribution of the sources, such that it covers the whole sky up to the confusion limit. The aforementioned approach is adopted here: we use the combination of the NVSS and SUMSS surveys. For those regions of sky covered by these surveys, especially for the brighter sources which these surveys are essentially complete, we will take the point sources from the actual catalog with measured position and flux. At the same time, we also derive the statistical distribution of the sources and extrapolate it to faint fluxes. For the regions which are not covered, or for fainter sources which the two surveys are incomplete, we will generate mock point sources according to the statistical distribution.

Note that when we constructed the whole sky diffuse maps, except for a few very bright ones, the majority of the point sources are not individually detected in the input data, so the diffuse map actually included the radiation from these sources. If we are to include the point sources, to avoid double counting, we should remove their radiation from the diffuse map to ensure that the model is self-consistent. We do this by computing the integrated flux of all the point sources for each pixel of the whole sky map and subtract them out.

Below we present the detailed procedures for constructing the diffuse map and the point source contributions.

3 Diffuse Emission

For the whole sky, we pixelize the sky map using the HEALPix scheme [33]. In HEALPix the whole sky is divided into $12n_{\rm side}^{2}$ pixels, the pixel size of the map is $\theta_{\rm pix}\approx\sqrt{4\pi/(12n_{\rm side}^{2})}$. The pixel size should not be larger than the FWHM resolution of the sky map, i.e. $\theta_{\rm pix}\leq\theta_{\rm FWHM}$. In the current model we use $n_{\rm side}=512$, which can support resolutions of up to $\sim 0.1^{\circ}$. If higher resolutions are needed, we can pixelize the sky with larger $n_{\rm side}$.

First we construct the whole sky diffuse foreground map. To model the diffuse emission, the GSM uses the observation data of whole sky surveys or large area surveys at frequency bands ranging from 10 MHz to 94 GHz (10 MHz to 5 THz for improved GSM), all with different angular resolutions at different frequencies [17, 19]. The dataset is listed in Table 1 (only up to 94 GHz for our purpose), and the maps of different frequencies are shown in Fig. 1. The principal components are obtained by solving the eigenvalue or singular value problem of the image correlation matrix, and the sky map at arbitrary frequencies are obtained by extrapolating these components, with the Haslam map at 408 MHz as the base map. As such, the principal components are purely mathematical and do not have obvious physical origins.

However, we consider a simple physical model of the foreground radiation, which includes the two major physical components of diffuse foregrounds at low frequencies: the galactic synchrotron emission and the galactic free-free emission [15]. The galactic synchrotron emission arises from the acceleration of the free high-energy (cosmic ray) electrons in the galactic magnetic field. The energy spectrum of the cosmic ray electrons is approximately a power-law $E(\gamma)\propto{\gamma}^{-p}$ [41, 42, 43], the brightness temperature of the synchrotron emission is also power-law. The galactic free-free (thermal) emission is produced by the scattering of free electrons with diffuse warm ionized gas. Although it is a small contribution compared with the synchrotron at low frequency [15], its brightness and fluctuation are still much larger than the 21cm signal. At high galactic latitude, H$\alpha$ and the diffuse warm ionized gas are optically thin and in local thermodynamic equilibrium, they are both proportional to the emission measure (EM), so we can use the galactic H$\alpha$ as a tracer of the galactic free-free emission [44].

The brightness temperature of the diffuse emission can be described as a power law

[TABLE]

For the spectrum index $\beta$, at the high galactic latitude, [45] obtained $\beta=2.81\pm 0.16$ in 1–10 GHz, they then combined it with the Haslam 408 MHz map [22] and Reich & Reich 1420 MHz map [46, 47] to obtain $\beta=2.76\pm 0.11$ in 0.4–7.5 GHz. At low frequencies, such as 150 MHz, in the same sky region as above, $\beta\approx 2.55$ and increases as frequency increases. At about 1 GHz, it steepens to $\beta=$ 2.8–3 [48, 49, 50, 51, 52, 46, 53, 45].

The spectrum index $\beta$ varies not only with frequency but also with different directions $\vec{n}$. We have several observations at mid- and low-frequencies from 10 MHz to 2.3 GHz, a feasible method to obtain $\beta(\nu,\vec{n})$ is to use the polynomial fitting [17].

Expanding $\beta(\nu,\vec{n})$ in the logarithm of the frequency,

[TABLE]

rewriting in matrix form,

[TABLE]

For each line of sight, $y$ is a ${\rm N}_{\nu}\times 1$ vector of the data, and $y=\left[\ln\frac{T_{\rm diffuse}(\vec{n},\nu_{1})}{T_{\rm diffuse}(\vec{n},\nu_{*})},\ln\frac{T_{\rm diffuse}(\vec{n},\nu_{2})}{T_{\rm diffuse}(\vec{n},\nu_{*})},\cdots\right]^{\rm T}$, $x$ is a ${\rm N}_{\beta}\times 1$ vector of the estimator, and $x=[\beta_{0}(\vec{n}),\beta_{1}(\vec{n}),\cdots]^{\rm T}$, A is a ${\rm N}_{\nu}\times{\rm N}_{\beta}$ matrix which is constructed as

[TABLE]

where ${\rm N}_{\nu}$ and ${\rm N}_{\beta}$ are the number of frequency bins and the number of parameters required to describe the spectrum index of diffuse emission respectively. The extra term $n$ denotes the noise with $\langle n\rangle=0$. The estimator $\hat{x}$ is the minimum variance solution of Eq. (7)

[TABLE]

To compute Eq. (9), we use the observation maps listed in Table 1, the procedure for obtaining the diffuse foreground map is as follows:

Pre-process the observation maps. As the 23 GHz, 33 GHz, 41 GHz, 61 GHz, 94 GHz maps are the WMAP-9yr diffuse radiation map products, we use them directly. The 2300 MHz map used here has already subtracted out the CMB temperature. For the other maps, we subtract the average CMB temperature 2.7255 K. The CMB anisotropy [54] is too small to affect any result. 2. 2.

Subtract the contribution of radio point sources (described in next section) from the maps. Note that in the 1420 MHz map some of the strongest sources such as Cassiopeia A and Cygnus A have already been removed. For each map we generate the point sources as described in section 4, first, generate a point source map at the same frequency as the sky map, with a resolution of $45^{\prime\prime}$; second, we use a circular Gaussian function (Eq. 14) or a elliptic Gaussian function (Eq. 15) to convolve the point source map if the sky map has a circular beam $\theta_{0}$ or a non-circular beam $\theta_{1}\times\theta_{2}$, we then obtained a point source map with the same resolution as the sky map; finally we subtract the point source map from the sky map. 3. 3.

Maps produced by following the above two steps (pre-processed maps) represent the diffuse emission. In order to compute and fit the spectrum indices, as the input maps are in different angular resolutions, all of them are smoothed to a $5^{\circ}$ FWHM resolution, which is the lowest resolution (45 MHz map) among all maps. But for the diffuse map simulation, we use the pre-processed 408 MHz map in resolution $0.85^{\circ}$ as the base.

We find that for the expansion in Eq. (6), an order of 6 is a good choice, the fitting curve is shown in Fig. 2, and the $\beta_{i}$ s are given in Table 2.

With the spectrum index of diffuse emission fitted above, we can simulate the diffuse emission map by scaling the pre-processed 408 MHz map obtained in item 3. As an example, Fig. 3 is the simulated diffuse map at 1.4 GHz.

4 Point Source

Most radio point sources are extragalactic. It has been noted that for flux above mJy level, most radio sources are radio-loud galaxies (radio-loud AGNs), quasars and BL Lac objects; while at sub-mJy, the star-forming (late-type) galaxies (mostly spiral galaxies) dominates [55, 56]. Recent studies suggest that there is another population, radio-quiet AGNs, which are star-forming galaxies but hosting an active nucleus [57, 58]. However, here we are mainly concerned with the flux density distribution???regardless the nature of the sources.

We derive the statistical distribution of the radio point sources from two surveys, namely the NRAO VLA Sky Survey (NVSS; [24]) and the Sydney University Molonglo Sky Survey (SUMSS; [59, 60]). The combination of these two surveys covers nearly the entire sky, and their angular resolutions and flux limits are comparable. It is more difficult to incorporate other survey data, e.g. the 3C [61], 5C [62], FIRST [63, 64] surveys, etc., because their sky coverages, resolutions and sensitivities are significantly different from each other.

The NVSS covers the sky of $\delta\geq-40^{\circ}$, which is approximately 82% of the celestial sphere, its observation frequency is 1.4 GHz, with a resolution $\theta_{\rm FWHM}=45^{\prime\prime}$ and nearly uniform sensitivity. It produces a catalog of almost $2\times 10^{6}$ discrete sources stronger than $S\approx 2.5$ mJy. SUMSS covers the sky of $\delta\leq-30^{\circ}$, its observation frequency is 843 MHz, and the angular resolution (FWHM) is $45^{\prime\prime}\times 45^{\prime\prime}{\rm cosec}|\delta|$ where $\delta$ is the declination, it produces a catalog of more than $2\times 10^{5}$ discrete sources stronger than $S\approx 6$ mJy.

In Fig. 4 we show the source count distributions $n(S)$ of the two surveys, and we see that above a certain threshold ($S=2.7$ mJy for NVSS and $S=12$ mJy for SUMSS), the $n(S)$ follows a power law distribution. The sample of the survey is incomplete when the flux drops below the threshold (sensitivity limit) and the density drops steeply at this point. Therefore, in the following analysis we shall only use the sample above the threshold.

The survey regions of NVSS and SUMSS overlap for $-40^{\circ}<\delta<-30^{\circ}$, and the resolutions of NVSS and SUMSS in the overlap region are similar ($\sim 45^{\prime\prime}$). The sources of the two catalogs in the overlap region are cross-identified using the following criterion: if the coordinate difference is within $30^{\prime\prime}$ and there is only one match within the surrounding region of $90^{\prime\prime}$, the two sources are identified as the same one. In the overlap region, 7515 matched sources are found, in agreement with the SUMSS survey team [60]. We compute the spectrum index $\alpha$ of these sources from the observations at 1.4 GHz (NVSS) and 843 MHz (SUMSS). The distribution of the spectrum index is well fitted by a Gaussian function, as shown in Fig. 5, with $\alpha=0.8157\pm 0.3885$.

With the spectrum index of flux density $\alpha$ fitted above (the index of brightness temperature is $\beta=2+\alpha$), we may select a conservative sample limit above which the survey is complete with high confidence. For NVSS we set this flux limit as $S=15$ mJy, and by scaling the corresponding SUMSS limit is $S=22$ mJy, i.e. $S_{1.4\,{\rm GHz}}=S_{843{\rm MHz}}\left(\frac{1.4\,{\rm GHz}}{843\,{\rm MHz}}\right)^{-\alpha}$ where $\alpha=0.8157$, the surface densities of sources in the two samples are then similar. We assume the NVSS and SUMSS surveys obtained complete samples of low frequency radio sources above this limit, then by combining the two surveys, we obtain a complete catalog of radio point sources above this limit for most of the sky. Thus, if the point source catalog is to be limited above 15 mJy at 1.4 GHz, we can use this actual catalogue to generate the point sources.

However, the confusion limit for our input all-sky map is well below this value, which is only 0.5 mJy (see Appendix A). We need to take into account of the point sources in the flux range of 0.5–15 mJy. For these sources the complete catalog is not available, but we can generate mock point sources according to their statistical distribution. The $n(S)$ distribution can be derived from the existing data, as detailed in Appendix B. For a sky region whose solid angle is $\Omega$, the total number of sources is $N_{\rm tot}=\int_{S}n(S)\Omega\,{\rm d}S$. There are about $2.5\times 10^{8}$ such sources in the whole sky, its distribution can be written approximately as

[TABLE]

The simplest way for generating the mock point sources would be to place them with a uniform random distribution on the celestial sphere. However, it has been observed that there is angular clustering among the radio sources, as would be expected for general matter distribution [65, 16]. To generate point sources with such angular clustering, we use the Rayleigh-Lévy random walk method, in which the distribution of the angular distance is [66, 67, 68, 14, 69]

[TABLE]

The parameters $\theta_{0}$ and $\gamma$ can be obtained from the two-point correlation function between NVSS and FIRST. Here we follow the result from [65] to use $\theta_{0}=6^{\prime}$ and $\gamma=0.8$.

We execute the following steps to generate the point sources of the Rayleigh-Lévy random walk:

Start from an arbitrary position, place the first source; 2. 2.

Choose a random direction $\vec{n}$, a random angular distance $\theta$ and an uniform random number $u$. If $u>P(>\theta)$, then place the next source at the coordinate $(\vec{n},\theta)$; 3. 3.

Repeat the Step 2, until all sources are placed; 4. 4.

Specify the flux of each source at 1.4 GHz. For each flux bin centered at $S_{0}$ with width $\Delta S_{0}$, the number of sources over the whole sky is $N_{0}=4\pi\int_{S_{0}-\Delta S_{0}/2}^{S_{0}+\Delta S_{0}/2}n(S)\,{\rm d}S$. Therefore, we randomly select $N_{0}$ sources from the mock map, give flux densities as uniform random values from $S_{0}-\Delta S_{0}/2$ to $S_{0}+\Delta S_{0}/2$, with the spectrum index given by a Gaussian random number with $N(\mu=0.8157,\sigma^{2}=0.1509)$.

Indeed, even for flux density above 15 mJy at 1.4 GHz, there are a few gaps on the sky not covered by the combined NVSS+SUMSS catalog. The same procedure can be used to generate sources which fill these gaps.

5 Results and Discussions

Now we can make the total sky map by summing up the diffuse foreground map (which is itself the sum of Galactic synchrotron and galactic free-free emission map) and the point sources. In the sky regions covered by the NVSS and SUMSS surveys, the bright sources are from the actual catalog. For sky regions not covered by the surveys, or for the fainter sources, point sources are generated according to statistical distribution. Using the spectrum indices we obtained for the diffuse maps and the point sources, the flux at the desirable frequencies are computed, and then convolved with a beam of assumed resolution.

We produce the total sky maps at 150 MHz, 408 MHz, 750 MHz, 1.4 GHz, and 2.3 GHz with FWHM resolutions $0.1^{\circ}$ and $0.015^{\circ}$ (shown in Fig. 6), as well as a point source catalogue above 0.5 mJy at 1.4 GHz. The sky appears brighter at lower frequencies, but in the frequency range discussed here the large features are similar. Our sky maps and point source catalogs are available for download333http://tianlai.bao.ac.cn/~huangqizhi/SSM/.

In the sky regions covered by the NVSS and SUMSS surveys, the brighter sources are from the actual catalogue. For sky regions not covered by the surveys, or for the fainter sources, point sources are generated according to statistical distribution. Using the spectrum indices we obtained for the diffuse maps and the point sources, the flux at the desirable frequencies are computed, and then convolved with a beam of assumed resolution.

In Fig. 7 we show a zoom up at the mid-declination for the maps at the different resolutions of $0.7^{\circ}$, $0.1^{\circ}$, and $0.015^{\circ}$. We can see that as the resolution improved, more radio sources can be seen clearly in the map, some sources which are merged in the low resolution maps appear distinct at higher resolutions, while the overall temperature of the map remains the same.

In Fig. 8 we plot the relative differences of the GSM (2008 and 2016) and our model with respect to the 2014 version of the Haslam map at 408 MHz in a resolution of \sim$$1^{\circ}. Compared with the original 1982 version, in the 2014 version of the Haslam map [28] strong baseline stripings are removed, however the stronger radio sources are also removed. So we can see that the GSM 2008 model which is based on the original 1982 version of Haslam map show differences in the form of stripes and points. The GSM 2016 model is itself based on the 2014 version of Haslam map, so there is very small difference. Our own model also has small differences with respect to the 2014 version of Haslam map, but as we added back the point sources, we can see some differences appearing as points.

With the addition of the point sources, we expect significant changes in the angular power spectrum. In Fig. 9 we plot the angular power spectrum for both our diffuse map and total map, at $1^{\circ}$ and $0.1^{\circ}$ resolutions. For comparison, we also use the CORA code to generate mock maps at these resolutions and plot the resulting angular power spectrum in this figure. For the diffuse map, the angular power is basically flat until it begins to drop at the scale corresponding to the map resolution. With the addition of the point sources, there are significant powers on the small scales. The angular power spectrum actually raises before reaching the map resolution. The general shape of our angular power spectrum is similar to that of the CORA map, though there are slight differences in the amplitude. As the two maps are both generated with random distribution, and there are differences in the details of both the assumed distribution and the method of generation, such difference is not unexpected. Compared with the low resolution GSM model, such angular power will have substantial impacts on the foreground removal results.

Note however our model consists a low resolution diffuse sky model and a point source catalog model, but no high resolution extended sources, so the angular power spectrum may not be completely accurate. Extended sources do exist, e.g. some galactic supernova remnants. Many extragalactic sources, e.g. radio galaxies, radio halo of galaxy clusters and groups, etc., are also extended source under higher resolutions. Such extended sources could be modeled in a manner similar to our current modeling of point sources but with more parameters each: in addition to the center location coordinates, total flux and spectral index, also the parameters describing its shapes (modeled e.g. with one or more elliptic Gaussian profiles) [27].

There are some limitations in the model presented here. In terms of applicable frequency range, here we considered only below 2.3 GHz. At this mid-to-low frequency range, the main components of the foreground radiation can be modeled using the galactic synchrotron emission, galactic free-free emission and the extragalactic sources which are also primarily radiating by these mechanism. We have not considered the propagation effects (absorption, refraction and scattering) in this treatment. At the very low frequencies these effects may have significant impacts. For example, the radiowaves at a few MHz may suffer significant absorption because of the interstellar medium electrons. These shall be studied in the future.

We have considered only unpolarized emission in this study. In fact, the galactic synchrotron emission is highly polarized, and some point sources are also polarized [70]. The polarized emission is also very important in the foreground subtraction analysis, because the Faraday rotation induce an oscillation along the frequency which is similar to the 21cm signal [71, 72, 73, 74, 75, 76].

Despite these limitations, the inclusion of the point source in this model increased the usable range of angular resolution, making it a useful tool for the study of foreground subtraction in 21cm experiment. Our model can be further improved in the future by implementing the extended sources, polarization and propagation effects.

*We thank Réza Ansari for discussions. The computations of this paper is performed on the Tianhe-2 supercomputer at the Guangzhou supercomputing center with the support of NSFC-Guangdong Joint grant U1501501. This research is supported by the Ministry of Science and Technology grant 2016YFE0100300, NSFC grants 11473044, 11761141012, 11633004, 11653003, and CAS grant QYZDJ-SSW-SLH017. F. Q. Wu also acknowledge the support by the CSC Cai Yuanpei grant. *

\InterestConflict

The authors declare that they have no conflict of interest.

Appendix A Computation of the Confusion Limit

The beam of a telescope can be modeled as a circular Gaussian function

[TABLE]

where $\theta_{0}$ is the FWHM resolution, $\lambda$ is the observation wavelength, $D$ is the aperture size of the telescope, or an elliptic Gaussian function

[TABLE]

For a telescope with a beam response $B(\theta,\varphi)$, a source with flux $S$ will induce an output $x=S\cdot B(\theta,\varphi)$, the average number of response $R(x)$ with an amplitude between $x$ and $x+{\rm d}x$ in a solid angle $\Omega_{b}$ is then [29, 30, 31]

[TABLE]

For a circular Gaussian beam, the integration of Eq. (16) is given by

[TABLE]

Let $S_{c}=Q\sigma_{c}$ where Q is a factor of 3–5 [32], the rms confusion $\sigma_{c}$ is [30, 31]

[TABLE]

If the differential count of radio sources $n(S)$ is a single power-law as given in Eq. (3), Eq. (18) becomes

[TABLE]

For broader range of frequencies, we may model $n(S)$ using a piece-wise power law [24, 60],

[TABLE]

The best-fit parameters are

[TABLE]

With these parameters we obtain the confusion limit, we find

[TABLE]

where $\alpha$ is the spectrum index of flux for the sources, typically $\alpha=$ 0.7–0.8 (Fig. 5), For example, for the NVSS+SUMSS catalogue, $\theta_{0}=45^{\prime\prime}=0.75^{\prime}$, if choose $Q=5$, the confusion limit $S_{c}=0.5$ mJy and the rms $\sigma_{c}=0.1$ mJy. So point source whose flux density above 0.5 mJy at 1.4 GHz can be distinguished.

Appendix B Estimation of Differential Surface Density of Sources

In order to calculate $n(S)$, we first compute the radio luminosity function. Two methods are usually used to calculate the local luminosity function, the maximum volume $V_{m}$ method [55, 77] and the STY maximum-likelihood method [78]. Here we use the maximum volume $V_{m}$ method, in which the target source is moved through the redshifts to find the maximum volume it could have inhabited and still been observed:

[TABLE]

where $N$ is the number of sources. Assuming the distribution of sources is random, the root mean square (r.m.s.) error is

[TABLE]

We use the data from [79], with the cosmological parameters from [80].

For start-forming galaxies, the luminosity function is of the form[81]

[TABLE]

The best-fit parameters of Eq. (30) are

[TABLE]

For radio-loud AGNs, the luminosity function is of the form[82]

[TABLE]

The best-fit parameters of Eq. (32) are

[TABLE]

Figure 10 shows the fitted local radio luminosity function.

For an isotropic source with spectrum index $\alpha\equiv-{\rm d}[\ln(S)]/{\rm d}[\ln(\nu)]$ at redshift $z$, its flux is given by

[TABLE]

where $L$ is the luminosity, and the luminosity distance of this source is $D_{L}=\frac{c(1+z)}{H_{0}}\int_{0}^{z}\frac{{\rm d}z^{{}^{\prime}}}{\sqrt{\Omega_{\Lambda}+\Omega_{m}(1+z^{{}^{\prime}})^{3}}}$, and we adopt the $\Lambda$CDM cosmological model with the parameters from Plank 2015 [80]. Then we have

[TABLE]

In a comoving distance $r$ to $r+{\rm d}r$, the number of sources with luminosities $L$ to $L+{\rm d}L$ is $\rho_{m}(L,z)A{\rm d}L{\rm d}r$. This number can also be expressed by the sources with flux densities $S$ to $S+{\rm d}S$ in redshift range $z$ to $z+{\rm d}z$, that is $\eta(S,z){\rm d}S{\rm d}z$. The comoving distance is given by $dr=cdz/H(z)$. Let $\eta(S,z){\rm d}S{\rm d}z=\rho_{m}(L,z)A{\rm d}L{\rm d}r$, we obtain

[TABLE]

The weighted differential count is

[TABLE]

where $n(S)$ is the number of sources per steradian per flux density.

We model the redshift evolution of the radio luminosity function as a combination of “pure luminosity evolution” and “pure density evolution” [55],

[TABLE]

where

[TABLE]

We use the full-sky complete sample from NVSS and SUMSS to fit the evolution of radio luminosity function, the weighted different count of radio sources (Eq. 37) is shown in Fig. 11, and the best-fit parameters are

[TABLE]

If approximate as a single power law, we dind $n(S)\approx 1300\,S^{-1.77}$.