A high-resolution self-consistent whole sky foreground model
Qizhi Huang Fengquan Wu Xuelei Chen

subject




Article 

2018
October
??
?
?
000000
September 30, 2018
??, 2018





A high-resolution self-consistent whole sky foreground model





1,2]Qizhi Huang
1]Fengquan Wu
1,3,4]Xuelei Chenxuelei@cosmology.bao.ac.cn



Huang Q Z


Huang Q Z, Wu F Q & Chen X L




[1]Key Laboratory of Computational Astrophysics, National Astronomical Observatories, Chinese Academy of Sciences, Beijing 100101, China
[2]Université Paris-Sud, LAL, UMR 8607, F-91898 Orsay Cedex, France & CNRS/IN2P3, F-91405 Orsay, France
[3]School of Astronomy and Space Science, University of Chinese Academy of Sciences, Beijing 100049, China
[4]Center for High Energy Physics, Peking University, Beijing 100871, China









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. 





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

A high-resolution self-consistent whole sky foreground model
    [
    December 30, 2023
============================================================









2








§ INTRODUCTION


Tomographic observation of the redshifted 21cm signal from neutral hydrogen will allow us to study the 
Epoch of Reionization (EoR) <cit.> and the Cosmic Dawn <cit.>, 
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 <cit.>, 
and even to probe the cosmic dark age <cit.>. 
A number of telescopes have been developed to detect the cosmological 21cm signal, including the EoR experiments such as
LOFAR <cit.>, MWA <cit.>, PAPER <cit.>,  
and HERA <cit.>;  
and the dark energy experiments such as Tianlai <cit.> and CHIME <cit.>;   
and the next generation telescopes such as SKA <cit.> 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 <cit.>. 
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 
<cit.>. 
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; <cit.>) is widely used for this purpose, 
with a recent updated version <cit.>. 
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 <cit.> 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^∘), 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 <cit.>.
It uses the 2003 version of the Haslam map [<https://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 <cit.> and NVSS <cit.> surveys, plus a mock sample of randomly generated
point sources according to the distribution given in Ref.<cit.>, and finally, 
based on Ref.<cit.>,
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; <cit.>), 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 <cit.>[<http://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.













§ 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 Θ is then given by 

    I_Θ(n⃗,ν) = I^f_Θ(n⃗,ν) + ∑_p ∈M(Θ) I^p_Θ(n⃗,ν)

where n⃗ denotes the direction of the sky, ν denotes the frequency, I^f_Θ denotes
 the whole sky diffuse foreground, 
M_Θ denotes the set of point sources at resolution Θ, and I^p_Θ denotes the 
radiation from point source p convolved with a beam of angular resolution Θ. 
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 Θ, but at 
a lower angular resolution Θ_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. 

    I_Θ(n⃗,ν)=  I^f_Θ_0(n⃗,ν) + ∑_p ∈M(Θ) I^p_Θ(n⃗,ν)

For example, for the low frequency sky map we use the Haslam map as a basis, which has Θ_0=0.85^∘. 
The construction of our diffuse foreground model is described in <ref>.

Compared with Eq. (<ref>), the model given by  Eq. (<ref>) 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 Θ and Θ_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 <cit.>, 
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,

    n(S) = k  S^-γ

For a circular Gaussian beam, the confusion limit flux is <cit.>

    S_c = [ Q^2  k π θ_0^2/(3-γ)(γ-1)  4ln2]^1/γ-1
	
where θ_0 = 1.029 λ/D, Q is a factor of 3–5 <cit.>. 
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 <ref>. 

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.















§ DIFFUSE EMISSION



For the whole sky, we pixelize the sky map using the HEALPix scheme <cit.>. 
In HEALPix the whole sky is divided into 12 n_ side^2 pixels, the pixel size of the map is 
θ_ pix≈√(4π/(12 n_ side^2)). 
The pixel size should not be larger than the FWHM resolution of the sky map, i.e. θ_ pix≤θ_ FWHM. 
In the current model we use n_ side=512, which can support resolutions of up to ∼ 0.1^∘. 
If higher resolutions are needed, we can pixelize the sky with larger n_ 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 <cit.>. 
The dataset is listed in <ref> (only up to 94 GHz for our purpose), and the maps of different 
frequencies are shown in Fig. <ref>.
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 <cit.>. 
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(γ) ∝γ^-p <cit.>, 
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 <cit.>, 
its brightness and fluctuation are still much larger than the 21cm signal. 
At high galactic latitude, Hα 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α as a tracer of the
galactic free-free emission <cit.>. 


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

    T_ diffuse(ν) = T_ diffuse(ν_*) ( ν/ν_*)^-β

For the spectrum index β, at the high galactic latitude, <cit.> obtained 
β=2.81±0.16 in 1–10 GHz, they then combined it with the Haslam 408 MHz map <cit.> and 
Reich & Reich 1420 MHz map <cit.> to obtain β=2.76±0.11 in 0.4–7.5 GHz. 
At low frequencies, such as 150 MHz, in the same sky region as above, β≈ 2.55 and increases as 
frequency increases. At about 1 GHz, it steepens to β = 2.8–3 <cit.>.


The spectrum index β varies not only with frequency but also with different directions n⃗.
We have several observations at mid- and low-frequencies from 10 MHz to 2.3 GHz, 
a feasible method to obtain β(ν,n⃗) is to use the polynomial fitting <cit.>. 





Expanding β(ν,n⃗) in the logarithm of the frequency, 

    ln[T_ diffuse(n⃗,ν)/T_ diffuse(n⃗,ν_*)] = 
    			β_0(n⃗)[lnν/ν_*]
    	    + β_1(n⃗)[lnν/ν_*]^2  +⋯
	
rewriting  in matrix form,

    y = A x + n

For each line of sight, y is a N_ν× 1 vector of the data, and
y=[ lnT_ diffuse(n⃗,ν_1)/T_ diffuse(n⃗,ν_*), lnT_ diffuse(n⃗,ν_2)/T_ diffuse(n⃗,ν_*), ⋯]^ T, 
x is a N_β× 1 vector of the estimator, and
x=[β_0(n⃗), β_1(n⃗), ⋯]^ T, 
A is a N_ν× N_β matrix which is constructed as

    A = [           lnν_1/ν_*       [lnν_1/ν_*]^2                   ⋯    [lnν_1/ν_*]^ N_β;           lnν_2/ν_*       [lnν_2/ν_*]^2                   ⋯    [lnν_2/ν_*]^ N_β;                   ⋮                                       ⋱;        lnν_ N_ν/ν_*    [lnν_ N_ν/ν_*]^2                   ⋯ [lnν_ N_ν/ν_*]^ N_β ]

where N_ν and N_β 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 ⟨ n ⟩=0.
The estimator x̂ is the minimum variance solution of Eq. (<ref>)

    x̂ = [A^ TN^-1A]^-1A^ TN^-1 y





To compute Eq. (<ref>), we use the observation maps listed in <ref>, 
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 <cit.> is too small to affect any result. 


  * 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 <ref>,
first, generate a point source map at the same frequency as the sky map, with a resolution of 45”; 
second, we use a circular Gaussian function (Eq. <ref>) or a elliptic Gaussian function (Eq. <ref>)
to convolve the point source map if the sky map has a circular beam θ_0 or a non-circular beam θ_1 ×θ_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. 


  * 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^∘ 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^∘ as the base. 







We find that for the expansion in Eq. (<ref>),  an order of 6 is a good choice, the fitting curve is shown in 
Fig. <ref>, and the β_i s are given in <ref>.

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 <ref>. As an example, 
Fig. <ref> is the simulated diffuse map at 1.4 GHz. 













§ 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 <cit.>. 
Recent studies suggest that there is another population, radio-quiet AGNs, which are star-forming galaxies but hosting 
an active nucleus <cit.>. 
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; <cit.>) and the Sydney University Molonglo Sky 
Survey (SUMSS; <cit.>). 
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 <cit.>, 5C <cit.>, 
FIRST <cit.> surveys, etc., 
because their sky coverages, resolutions and sensitivities are significantly different from each other. 

The NVSS covers the sky of δ≥ -40^∘, which is approximately 82% of the celestial sphere, 
its observation frequency is 1.4 GHz, with a resolution θ_ FWHM = 45” and nearly uniform sensitivity. 
It produces a catalog of almost 2 × 10^6 discrete sources stronger than S ≈ 2.5 mJy. 
SUMSS covers the sky of  δ≤ -30^∘, its observation frequency is 843 MHz, and
the angular resolution (FWHM) is 45”× 45” cosec|δ| where δ is the declination, 
it produces a catalog of more than 2 × 10^5 discrete sources stronger than S ≈ 6 mJy. 
	





In Fig. <ref> 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^∘ < δ < -30^∘,
and the resolutions of NVSS and SUMSS in the overlap region are similar (∼ 45”). The sources of the two catalogs 
in the overlap region are cross-identified using the following criterion: if the coordinate difference
is within  30” and there is only one match within the surrounding region of 90”, 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 <cit.>.
We compute the spectrum index α 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. <ref>, with α=0.8157 ± 0.3885.


With the spectrum index of flux density α fitted above (the index of brightness temperature is β=2+α), 
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  GHz = S_843 MHz(1.4  GHz/843  MHz)^-α where α=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 <ref>). 
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 <ref>. 
For a sky region whose solid angle is Ω, the total number of sources is N_ tot=∫_S n(S) Ω  dS. 
There are about 2.5 × 10^8 such sources in the whole sky, its distribution can be written approximately as 

    n(S) ≈ 1300  S^-1.77


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 <cit.>. 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 <cit.>

    P(>θ) = {[ (θ/θ_0)^-γ ,        θ≥θ_0;         1  ,      θ < θ_0 ].

The parameters θ_0 and γ can be obtained from the two-point correlation function between NVSS and FIRST. 
Here we follow the result from <cit.> to use θ_0=6' and γ=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;


  * Choose a random direction n⃗, a random angular distance θ and an uniform random number u. 
If u > P(>θ), then place the next source at the coordinate (n⃗,θ);


  * Repeat the Step 2, until all sources are placed; 


  * Specify the flux of each source at 1.4 GHz. 
For each flux bin centered at S_0 with width Δ S_0, the number of sources over the whole sky is 
N_0 = 4π∫_S_0-Δ S_0/2^S_0+Δ S_0/2 n(S)  dS. 
Therefore, we randomly select N_0 sources from the mock map, 
give flux densities as uniform random values from S_0-Δ S_0/2 to S_0+Δ S_0/2, 
with the spectrum index given by a Gaussian random number with N(μ=0.8157, σ^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.  














§ 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^∘ and 0.015^∘ (shown in Fig. <ref>), 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 download[<http://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. <ref> we show a zoom up at the mid-declination for the maps at the different 
resolutions of 0.7^∘, 0.1^∘, and 0.015^∘. 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. <ref> 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 ∼1^∘.
Compared with the original 1982 version, in 
the 2014 version of the Haslam map <cit.> 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. <ref> 
we plot the angular power spectrum for both our diffuse map and total map, at 1^∘ and 0.1^∘ 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) <cit.>.
 



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 <cit.>. 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
<cit.>. 

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. 








The authors declare that they have no conflict of interest.













abbrv,sort



10

1997ApJ...475..429M
P. Madau, A. Meiksin, and M. J. Rees.
21 Centimeter Tomography of the Intergalactic Medium at High
  Redshift.
ApJ, 475:429–444, February 1997.

2003ApJ...596....1C
B. Ciardi and P. Madau.
Probing beyond the Epoch of Hydrogen Reionization with 21 Centimeter
  Radiation.
ApJ, 596:1–8, October 2003.

2015MNRAS.447.1806G
R. Ghara, T. R. Choudhury, and K. K. Datta.
21 cm signal from cosmic dawn: imprints of spin temperature
  fluctuations and peculiar velocities.
MNRAS, 447:1806–1825, February 2015.

2008PhRvL.100i1303C
T.-C. Chang, U.-L. Pen, J. B. Peterson, and P. McDonald.
Baryon Acoustic Oscillation Intensity Mapping of Dark Energy.
Physical Review Letters, 100(9):091303, March 2008.

2008PhRvD..78b3529M
Y. Mao, M. Tegmark, M. McQuinn, M. Zaldarriaga, and O. Zahn.
How accurately can 21cm tomography constrain cosmology?
Phys. Rev. D, 78(2):023529, July 2008.

2004PhRvL..92u1301L
A. Loeb and M. Zaldarriaga.
Measuring the Small-Scale Power Spectrum of Cosmic Density
  Fluctuations through 21cm Tomography Prior to the Epoch of Structure
  Formation.
Physical Review Letters, 92(21):211301, May 2004.

2013A   A...556A...2V
M. P. van Haarlem, M. W. Wise, A. W. Gunst, G. Heald, J. P. McKean,
  J. W. T. Hessels, A. G. de Bruyn, R. Nijboer, J. Swinbank,
  R. Fallows, M. Brentjens, A. Nelles, R. Beck, H. Falcke,
  R. Fender, J. Hörandel, L. V. E. Koopmans, G. Mann, G. Miley,
  H. Röttgering, B. W. Stappers, R. A. M. J. Wijers, S. Zaroubi,
  M. van den Akker, A. Alexov, J. Anderson, K. Anderson, A. van
  Ardenne, M. Arts, A. Asgekar, I. M. Avruch, F. Batejat,
  L. Bähren, M. E. Bell, M. R. Bell, I. van Bemmel, P. Bennema,
  M. J. Bentum, G. Bernardi, P. Best, L. Bîrzan, A. Bonafede,
  A.-J. Boonstra, R. Braun, J. Bregman, F. Breitling, R. H. van de
  Brink, J. Broderick, P. C. Broekema, W. N. Brouw, M. Brüggen,
  H. R. Butcher, W. van Cappellen, B. Ciardi, T. Coenen, J. Conway,
  A. Coolen, A. Corstanje, S. Damstra, O. Davies, A. T. Deller, R.-J.
  Dettmar, G. van Diepen, K. Dijkstra, P. Donker, A. Doorduin,
  J. Dromer, M. Drost, A. van Duin, J. Eislöffel, J. van Enst,
  C. Ferrari, W. Frieswijk, H. Gankema, M. A. Garrett, F. de
  Gasperin, M. Gerbers, E. de Geus, J.-M. Grießmeier, T. Grit,
  P. Gruppen, J. P. Hamaker, T. Hassall, M. Hoeft, H. A. Holties,
  A. Horneffer, A. van der Horst, A. van Houwelingen, A. Huijgen,
  M. Iacobelli, H. Intema, N. Jackson, V. Jelic, A. de Jong,
  E. Juette, D. Kant, A. Karastergiou, A. Koers, H. Kollen, V. I.
  Kondratiev, E. Kooistra, Y. Koopman, A. Koster, M. Kuniyoshi,
  M. Kramer, G. Kuper, P. Lambropoulos, C. Law, J. van Leeuwen,
  J. Lemaitre, M. Loose, P. Maat, G. Macario, S. Markoff,
  J. Masters, R. A. McFadden, D. McKay-Bukowski, H. Meijering,
  H. Meulman, M. Mevius, E. Middelberg, R. Millenaar, J. C. A.
  Miller-Jones, R. N. Mohan, J. D. Mol, J. Morawietz, R. Morganti,
  D. D. Mulcahy, E. Mulder, H. Munk, L. Nieuwenhuis, R. van
  Nieuwpoort, J. E. Noordam, M. Norden, A. Noutsos, A. R. Offringa,
  H. Olofsson, A. Omar, E. Orrú, R. Overeem, H. Paas,
  M. Pandey-Pommier, V. N. Pandey, R. Pizzo, A. Polatidis,
  D. Rafferty, S. Rawlings, W. Reich, J.-P. de Reijer, J. Reitsma,
  G. A. Renting, P. Riemers, E. Rol, J. W. Romein, J. Roosjen,
  M. Ruiter, A. Scaife, K. van der Schaaf, B. Scheers, P. Schellart,
  A. Schoenmakers, G. Schoonderbeek, M. Serylak, A. Shulevski,
  J. Sluman, O. Smirnov, C. Sobey, H. Spreeuw, M. Steinmetz, C. G. M.
  Sterks, H.-J. Stiepel, K. Stuurwold, M. Tagger, Y. Tang,
  C. Tasse, I. Thomas, S. Thoudam, M. C. Toribio, B. van der Tol,
  O. Usov, M. van Veelen, A.-J. van der Veen, S. ter Veen, J. P. W.
  Verbiest, R. Vermeulen, N. Vermaas, C. Vocks, C. Vogt, M. de Vos,
  E. van der Wal, R. van Weeren, H. Weggemans, P. Weltevrede,
  S. White, S. J. Wijnholds, T. Wilhelmsson, O. Wucknitz,
  S. Yatawatta, P. Zarka, A. Zensus, and J. van Zwieten.
LOFAR: The LOw-Frequency ARray.
A&A, 556:A2, August 2013.

2013PASA...30....7T
S. J. Tingay, R. Goeke, J. D. Bowman, D. Emrich, S. M. Ord, D. A.
  Mitchell, M. F. Morales, T. Booler, B. Crosse, R. B. Wayth, C. J.
  Lonsdale, S. Tremblay, D. Pallot, T. Colegate, A. Wicenec,
  N. Kudryavtseva, W. Arcus, D. Barnes, G. Bernardi, F. Briggs,
  S. Burns, J. D. Bunton, R. J. Cappallo, B. E. Corey, A. Deshpande,
  L. Desouza, B. M. Gaensler, L. J. Greenhill, P. J. Hall, B. J.
  Hazelton, D. Herne, J. N. Hewitt, M. Johnston-Hollitt, D. L.
  Kaplan, J. C. Kasper, B. B. Kincaid, R. Koenig, E. Kratzenberg,
  M. J. Lynch, B. Mckinley, S. R. Mcwhirter, E. Morgan, D. Oberoi,
  J. Pathikulangara, T. Prabu, R. A. Remillard, A. E. E. Rogers,
  A. Roshi, J. E. Salah, R. J. Sault, N. Udaya-Shankar,
  F. Schlagenhaufer, K. S. Srivani, J. Stevens, R. Subrahmanyan,
  M. Waterson, R. L. Webster, A. R. Whitney, A. Williams, C. L.
  Williams, and J. S. B. Wyithe.
The Murchison Widefield Array: The Square Kilometre Array Precursor
  at Low Radio Frequencies.
PASA, 30:7, January 2013.

2010AJ....139.1468P
A. R. Parsons, D. C. Backer, G. S. Foster, M. C. H. Wright, R. F.
  Bradley, N. E. Gugliucci, C. R. Parashare, E. E. Benoit, J. E.
  Aguirre, D. C. Jacobs, C. L. Carilli, D. Herne, M. J. Lynch, J. R.
  Manley, and D. J. Werthimer.
The Precision Array for Probing the Epoch of Re-ionization: Eight
  Station Results.
AJ, 139:1468–1480, April 2010.

2017PASP..129d5001D
D. R. DeBoer, A. R. Parsons, J. E. Aguirre, P. Alexander, Z. S. Ali,
  A. P. Beardsley, G. Bernardi, J. D. Bowman, R. F. Bradley, C. L.
  Carilli, C. Cheng, E. de Lera Acedo, J. S. Dillon, A. Ewall-Wice,
  G. Fadana, N. Fagnoni, R. Fritz, S. R. Furlanetto, B. Glendenning,
  B. Greig, J. Grobbelaar, B. J. Hazelton, J. N. Hewitt, J. Hickish,
  D. C. Jacobs, A. Julius, M. Kariseb, S. A. Kohn, T. Lekalake,
  A. Liu, A. Loots, D. MacMahon, L. Malan, C. Malgas, M. Maree,
  Z. Martinot, N. Mathison, E. Matsetela, A. Mesinger, M. F. Morales,
  A. R. Neben, N. Patra, S. Pieterse, J. C. Pober, N. Razavi-Ghods,
  J. Ringuette, J. Robnett, K. Rosie, R. Sell, C. Smith, A. Syce,
  M. Tegmark, N. Thyagarajan, P. K. G. Williams, and H. Zheng.
Hydrogen Epoch of Reionization Array (HERA).
PASP, 129(4):045001, April 2017.

2012IJMPS..12..256C
X. Chen.
The Tianlai Project: a 21CM Cosmology Experiment.
In International Journal of Modern Physics Conference Series,
  volume 12 of International Journal of Modern Physics Conference Series,
  pages 256–263, March 2012.

2014SPIE.9145E..22B
K. Bandura, G. E. Addison, M. Amiri, J. R. Bond, D. Campbell-Wilson,
  L. Connor, J.-F. Cliche, G. Davis, M. Deng, N. Denman, M. Dobbs,
  M. Fandino, K. Gibbs, A. Gilbert, M. Halpern, D. Hanna, A. D.
  Hincks, G. Hinshaw, C. Höfer, P. Klages, T. L. Landecker,
  K. Masui, J. Mena Parra, L. B. Newburgh, U.-l. Pen, J. B. Peterson,
  A. Recnik, J. R. Shaw, K. Sigurdson, M. Sitwell, G. Smecher,
  R. Smegal, K. Vanderlinde, and D. Wiebe.
Canadian Hydrogen Intensity Mapping Experiment (CHIME) pathfinder.
In Ground-based and Airborne Telescopes V, volume 9145 of procspie, page 914522, July 2014.

2013arXiv1311.4288H
M. Huynh and J. Lazio.
An Overview of the Square Kilometre Array.
ArXiv:1311.4288, November 2013.

2008MNRAS.389.1319J
V. Jelić, S. Zaroubi, P. Labropoulos, R. M. Thomas, G. Bernardi,
  M. A. Brentjens, A. G. de Bruyn, B. Ciardi, G. Harker, L. V. E.
  Koopmans, V. N. Pandey, J. Schaye, and S. Yatawatta.
Foreground simulations for the LOFAR-epoch of reionization
  experiment.
MNRAS, 389:1319–1335, September 2008.

1999A   A...345..380S
P. A. Shaver, R. A. Windhorst, P. Madau, and A. G. de Bruyn.
Can the reionization epoch be detected as a global signature in the
  cosmic background?
A&A, 345:380–390, May 1999.

2004MNRAS.355.1053D
T. Di Matteo, B. Ciardi, and F. Miniati.
The 21-cm emission from the reionization epoch: extended and point
  source foregrounds.
MNRAS, 355:1053–1065, December 2004.

2008MNRAS.388..247D
A. de Oliveira-Costa, M. Tegmark, B. M. Gaensler, J. Jonas, T. L.
  Landecker, and P. Reich.
A model of diffuse Galactic radio emission from 10 MHz to 100 GHz.
MNRAS, 388:247–260, July 2008.

2012MNRAS.419.3491L
A. Liu and M. Tegmark.
How well can we measure and understand foregrounds with 21-cm
  experiments?
MNRAS, 419:3491–3504, February 2012.

2017MNRAS.464.3486Z
H. Zheng, M. Tegmark, J. S. Dillon, D. A. Kim, A. Liu, A. R. Neben,
  J. Jonas, P. Reich, and W. Reich.
An improved model of diffuse galactic radio emission from 10 MHz to
  5 THz.
MNRAS, 464:3486–3497, January 2017.

2014ApJ...781...57S
J. R. Shaw, K. Sigurdson, U.-L. Pen, A. Stebbins, and M. Sitwell.
All-sky Interferometry with Spherical Harmonic Transit Telescopes.
ApJ, 781:57, February 2014.

2015PhRvD..91h3514S
J. R. Shaw, K. Sigurdson, M. Sitwell, A. Stebbins, and U.-L. Pen.
Coaxing cosmic 21 cm fluctuations from the polarized sky using m
  -mode analysis.
Phys. Rev. D, 91(8):083514, April 2015.

1982AandAS...47....1H
C. G. T. Haslam, C. J. Salter, H. Stoffel, and W. E. Wilson.
A 408 MHz all-sky continuum survey. II - The atlas of contour maps.
A&AS, 47:1, January 1982.

2007AJ....134.1245C
A. S. Cohen, W. M. Lane, W. D. Cotton, N. E. Kassim, T. J. W. Lazio,
  R. A. Perley, J. J. Condon, and W. C. Erickson.
The VLA Low-Frequency Sky Survey.
AJ, 134:1245–1262, September 2007.

1998AJ....115.1693C
J. J. Condon, W. D. Cotton, E. W. Greisen, Q. F. Yin, R. A. Perley,
  G. B. Taylor, and J. J. Broderick.
The NRAO VLA Sky Survey.
AJ, 115:1693–1716, May 1998.

2002ApJ...564..576D
T. Di Matteo, R. Perna, T. Abel, and M. J. Rees.
Radio Foregrounds for the 21 Centimeter Tomography of the Neutral
  Intergalactic Medium at High Redshifts.
ApJ, 564:576–580, January 2002.

2005ApJ...625..575S
M. G. Santos, A. Cooray, and L. Knox.
Multifrequency Analysis of 21 Centimeter Fluctuations from the Era
  of Reionization.
ApJ, 625:575–587, June 2005.

2018arXiv180505222B
A. Bonaldi, M. Bonato, V. Galluzzi, I. Harrison, M. Massardi,
  S. Kay, G. De Zotti, and M. L. Brown.
The Tiered Radio Extragalactic Continuum Simulation (T-RECS).
ArXiv:1805.05222, May 2018.

2015MNRAS.451.4311R
M. Remazeilles, C. Dickinson, A. J. Banday, M.-A. Bigot-Sazy, and
  T. Ghosh.
An improved source-subtracted and destriped 408-MHz all-sky map.
MNRAS, 451:4311–4327, August 2015.

2001MNRAS.325.1241V
P. Väisänen, E. V. Tollestrup, and G. G. Fazio.
Confusion limit resulting from galaxies: using the Infrared Array
  Camera on board SIRTF.
MNRAS, 325:1241–1252, August 2001.

2006MNRAS.369..281J
W.-S. Jeong, C. P. Pearson, H. M. Lee, S. Pak, and T. Nakagawa.
Far-infrared detection limits - II. Probing confusion including
  source confusion.
MNRAS, 369:281–294, June 2006.

2012ApJ...758...23C
J. J. Condon, W. D. Cotton, E. B. Fomalont, K. I. Kellermann,
  N. Miller, R. A. Perley, D. Scott, T. Vernstrom, and J. V. Wall.
Resolving the Radio Source Background: Deeper Understanding through
  Confusion.
ApJ, 758:23, October 2012.

1973ApJ...183....1M
H. S. Murdoch, D. F. Crawford, and D. L. Jauncey.
Maximum-Likelihood Estimation of the Number-Flux Distribution of
  Radio Sources in the Presence of Noise and Confusion.
ApJ, 183:1–14, July 1973.

2005ApJ...622..759G
K. M. Górski, E. Hivon, A. J. Banday, B. D. Wandelt, F. K.
  Hansen, M. Reinecke, and M. Bartelmann.
HEALPix: A Framework for High-Resolution Discretization and Fast
  Analysis of Data Distributed on the Sphere.
ApJ, 622:759–771, April 2005.

1976MNRAS.177..601C
J. L. Caswell.
A map of the northern sky at 10 MHz.
MNRAS, 177:601–616, December 1976.

1999AandAS..137....7R
R. S. Roger, C. H. Costain, T. L. Landecker, and C. M. Swerdlyk.
The radio emission from the Galaxy at 22 MHz.
A&AS, 137:7–19, May 1999.

2011AandA...525A.138G
A. E. Guzmán, J. May, H. Alvarez, and K. Maeda.
All-sky Galactic radiation at 45 MHz and spectral index between 45
  and 408 MHz.
A&A, 525:A138, January 2011.

1970AuJPA..16....1L
T. L. Landecker and R. Wielebinski.
The Galactic Metre Wave Radiation: A two-frequency survey between
  declinations +25^o and -25^o and the preparation of a map of the
  whole sky.
Australian Journal of Physics Astrophysical Supplement, 16:1,
  1970.

2001AandA...376..861R
P. Reich, J. C. Testori, and W. Reich.
A radio continuum survey of the southern sky at 1420 MHz. The atlas
  of contour maps.
A&A, 376:861–877, September 2001.

2013yCat..35560001T
C. Tello, T. Villela, S. Torres, M. Bersanelli, G. F. Smoot, I. S.
  Ferreira, A. Cingoz, J. Lamb, D. Barbosa, D. Perez-Becker,
  S. Ricciardi, J. A. Currivan, P. Platania, and D. Maino.
VizieR Online Data Catalog: The 2.3GHz continuum survey of the GEM
  project (Tello+, 2013).
VizieR Online Data Catalog, 355, May 2013.

2013ApJS..208...20B
C. L. Bennett, D. Larson, J. L. Weiland, N. Jarosik, G. Hinshaw,
  N. Odegard, K. M. Smith, R. S. Hill, B. Gold, M. Halpern,
  E. Komatsu, M. R. Nolta, L. Page, D. N. Spergel, E. Wollack,
  J. Dunkley, A. Kogut, M. Limon, S. S. Meyer, G. S. Tucker, and
  E. L. Wright.
Nine-year Wilkinson Microwave Anisotropy Probe (WMAP) Observations:
  Final Maps and Results.
ApJS, 208:20, October 2013.

2007ARNPS..57..285S
A. W. Strong, I. V. Moskalenko, and V. S. Ptuskin.
Cosmic-Ray Propagation and Interactions in the Galaxy.
Annual Review of Nuclear and Particle Science, 57:285–327,
  November 2007.

PhysRevD.82.092004
M. Ackermann, M. Ajello, W. B. Atwood, L. Baldini, J. Ballet, G. Barbiellini,
  D. Bastieri, B. M. Baughman, K. Bechtol, and F. Bellardi.
Fermi lat observations of cosmic-ray electrons from 7 gev to 1 tev.
Phys. Rev. D, 82:092004, Nov 2010.

2011A   A...534A..54S
A. W. Strong, E. Orlando, and T. R. Jaffe.
The interstellar cosmic-ray electron spectrum from synchrotron
  radiation and direct measurements.
A&A, 534:A54, October 2011.

1998astro.ph..1121S
G. F. Smoot.
Galactic Free-free and H-alpha Emission.
ArXiv:astro-ph/9801121, January 1998.

1998ApJ...505..473P
P. Platania, M. Bensadoun, M. Bersanelli, G. De Amici, A. Kogut,
  S. Levin, D. Maino, and G. F. Smoot.
A Determination of the Spectral Index of Galactic Synchrotron
  Emission in the 1-10 GHz Range.
ApJ, 505:473–483, October 1998.

1988A   AS...74....7R
P. Reich and W. Reich.
A map of spectral indices of the Galactic radio continuum emission
  between 408 MHz and 1420 MHz for the entire northern sky.
A&AS, 74:7–20, July 1988.

2001A   A...368.1123T
J. C. Testori, P. Reich, J. A. Bava, F. R. Colomb, E. E. Hurrel,
  J. J. Larrarte, W. Reich, and A. J. Sanz.
A radio continuum survey of the southern sky at 1420 MHz.
  Observations and data reduction.
A&A, 368:1123–1132, March 2001.

1967MNRAS.136..219B
A. H. Bridle.
The spectrum of the radio background between 13 and 404 MHz.
MNRAS, 136:219, 1967.

1974MNRAS.166..345S
G. Sironi.
The spectrum of the galactic non-thermal background radiation-1.
  Observations at 151-5 and 408 MHz.
MNRAS, 166:345–354, February 1974.

1974MNRAS.166..355W
A. S. Webster.
The spectrum of the galactic non-thermal background radiational
  Observations at 408, 610 and 1407 MHz.
MNRAS, 166:355–372, February 1974.

1979MNRAS.189..465C
H. V. Cane.
Spectra of the non-thermal radio radiation from the galactic polar
  regions.
MNRAS, 189:465–478, November 1979.

1987MNRAS.225..307L
K. D. Lawson, C. J. Mayer, J. L. Osborne, and M. L. Parkinson.
Variations in the Spectral Index of the Galactic Radio Continuum
  Emission in the Northern Hemisphere.
MNRAS, 225:307, March 1987.

1991MNRAS.248..705B
A. J. Banday and A. W. Wolfendale.
Fluctuations in the galactic synchrotron radiation. I - Implications
  for searches for fluctuations of cosmological origin.
MNRAS, 248:705–714, February 1991.

2016A   A...594A...9P
Planck Collaboration, R. Adam, P. A. R. Ade, N. Aghanim, M. Arnaud,
  M. Ashdown, J. Aumont, C. Baccigalupi, A. J. Banday, R. B.
  Barreiro, and et al.
Planck 2015 results. IX. Diffuse component separation: CMB maps.
A&A, 594:A9, September 2016.

1984ApJ...284...44C
J. J. Condon.
Cosmological evolution of radio sources found at 1.4 GHz.
ApJ, 284:44–53, September 1984.

2005PASA...22...36J
C. Jackson.
The Extragalactic Radio Sky at Faint Flux Densities.
pasa, 22:36–48, 2005.

2016ApJ...831..168K
K. I. Kellermann, J. J. Condon, A. E. Kimball, R. A. Perley, and Ž. Ivezić.
Radio-loud and Radio-quiet QSOs.
ApJ, 831:168, November 2016.

2017ApJ...842...95M
C. Mancuso, A. Lapi, I. Prandoni, I. Obi, J. Gonzalez-Nuevo,
  F. Perrotta, A. Bressan, A. Celotti, and L. Danese.
Galaxy Evolution in the Radio Band: The Role of Star-forming
  Galaxies and Active Galactic Nuclei.
ApJ, 842:95, June 2017.

1999AJ....117.1578B
D. C.-J. Bock, M. I. Large, and E. M. Sadler.
SUMSS: A Wide-Field Radio Imaging Survey of the Southern Sky. I.
  Science Goals, Survey Design, and Instrumentation.
AJ, 117:1578–1593, March 1999.

2003MNRAS.342.1117M
T. Mauch, T. Murphy, H. J. Buttery, J. Curran, R. W. Hunstead,
  B. Piestrzynski, J. G. Robertson, and E. M. Sadler.
SUMSS: a wide-field radio imaging survey of the southern sky - II.
  The source catalogue.
MNRAS, 342:1117–1130, July 2003.

1959MmRAS..68...37E
D. O. Edge, J. R. Shakeshaft, W. B. McAdam, J. E. Baldwin, and
  S. Archer.
A survey of radio sources at a frequency of 159 Mc/s.
Mem. RAS, 68:37–60, 1959.

1975MNRAS.171..475P
T. J. Pearson.
The 5C 5 survey of radio sources.
MNRAS, 171:475–505, June 1975.

1995ApJ...450..559B
R. H. Becker, R. L. White, and D. J. Helfand.
The FIRST Survey: Faint Images of the Radio Sky at Twenty
  Centimeters.
ApJ, 450:559, September 1995.

2015ApJ...801...26H
D. J. Helfand, R. L. White, and R. H. Becker.
The Last of FIRST: The Final Catalog and Source Identifications.
ApJ, 801:26, March 2015.

2003A   A...405...53O
R. A. Overzier, H. J. A. Röttgering, R. B. Rengelink, and R. J.
  Wilman.
The spatial clustering of radio sources in NVSS and FIRST;
  implications for galaxy clustering evolution.
A&A, 405:53–72, July 2003.

1975CRASM.280.1551M
B. Mandelbrot.
On a decomposable model of a hierarchical universe - Derivation of
  galactic correlations on the celestial sphere.
Academie des Sciences Paris Comptes Rendus Serie Sciences
  Mathematiques, 280:1551–1554, June 1975.

1995ApJS...96..401C
S. Colombi, F. R. Bouchet, and R. Schaeffer.
A count probability cookbok: Spurious effects and the scaling
  model.
ApJS, 96:401–428, February 1995.

2003AJ....125.2064G
R. R. Gal, R. R. de Carvalho, P. A. A. Lopes, S. G. Djorgovski, R. J.
  Brunner, A. Mahabal, and S. C. Odewahn.
The Northern Sky Optical Cluster Survey. II. An Objective Cluster
  Catalog for 5800 Square Degrees.
AJ, 125:2064–2084, April 2003.

2014MNRAS.440...10A
D. Alonso, A. Bueno Belloso, F. J. Sánchez,
  J. García-Bellido, and E. Sánchez.
Measuring the transition to homogeneity with photometric redshift
  surveys.
MNRAS, 440:10–23, May 2014.

2018A   A...613A..58V
C. L. Van Eck, M. Haverkorn, M. I. R. Alves, R. Beck, P. Best,
  E. Carretti, K. T. Chyży, J. S. Farnes, K. Ferrière, M. J.
  Hardcastle, G. Heald, C. Horellou, M. Iacobelli, V. Jelić,
  D. D. Mulcahy, S. P. O'Sullivan, I. M. Polderman, W. Reich, C. J.
  Riseley, H. Röttgering, D. H. F. M. Schnitzeler, T. W. Shimwell,
  V. Vacca, J. Vink, and G. J. White.
Polarized point sources in the LOFAR Two-meter Sky Survey: A
  preliminary catalog.
A&A, 613:A58, June 2018.

2010MNRAS.409.1647J
V. Jelić, S. Zaroubi, P. Labropoulos, G. Bernardi, A. G. de
  Bruyn, and L. V. E. Koopmans.
Realistic simulations of the Galactic polarized foreground:
  consequences for 21-cm reionization detection experiments.
MNRAS, 409:1647–1659, December 2010.

Asad:2015sda
Khan M. B. Asad et al.
Polarization leakage in epoch of reionization windows ? I. Low
  Frequency Array observations of the 3C196 field.
Mon. Not. Roy. Astron. Soc., 451(4):3709–3727, 2015.

2016MNRAS.462.4482A
K. M. B. Asad, L. V. E. Koopmans, V. Jelić, A. Ghosh, F. B.
  Abdalla, M. A. Brentjens, A. G. de Bruyn, B. Ciardi, B. K. Gehlot,
  I. T. Iliev, M. Mevius, V. N. Pandey, S. Yatawatta, and S. Zaroubi.
Polarization leakage in epoch of reionization windows - II. Primary
  beam model and direction-dependent calibration.
MNRAS, 462:4482–4494, November 2016.

2017A   A...597A..98V
C. L. Van Eck, M. Haverkorn, M. I. R. Alves, R. Beck, A. G. de Bruyn,
  T. Enßlin, J. S. Farnes, K. Ferrière, G. Heald,
  C. Horellou, A. Horneffer, M. Iacobelli, V. Jelić,
  I. Martí-Vidal, D. D. Mulcahy, W. Reich, H. J. A.
  Röttgering, A. M. M. Scaife, D. H. F. M. Schnitzeler, C. Sobey,
  and S. S. Sridhar.
Faraday tomography of the local interstellar medium with LOFAR:
  Galactic foregrounds towards IC 342.
A&A, 597:A98, January 2017.

2018MNRAS.476.3051A
K. M. B. Asad, L. V. E. Koopmans, V. Jelić, A. G. de Bruyn, V. N.
  Pandey, and B. K. Gehlot.
Polarization leakage in epoch of reionization windows - III.
  Wide-field effects of narrow-field arrays.
MNRAS, 476:3051–3062, May 2018.

Spinelli:2018kpm
M. Spinelli, G. Bernardi, and M. G. Santos.
Simulations of Galactic polarized synchrotron emission for Epoch of
  Reionization observations.
2018.

1987A   A...184....7T
L. Toffolatti, A. Franceschini, L. Danese, and G. de Zotti.
The local radio luminosity function of galaxies.
A&A, 184:7–15, October 1987.

1979ApJ...232..352S
A. Sandage, G. A. Tammann, and A. Yahil.
The velocity field of bright nearby galaxies. I - The variation of
  mean absolute magnitude with redshift for galaxies in a magnitude-limited
  sample.
ApJ, 232:352–364, September 1979.

2007MNRAS.375..931M
T. Mauch and E. M. Sadler.
Radio sources in the 6dFGS: local luminosity functions at 1.4GHz for
  star-forming galaxies and radio-loud AGN.
MNRAS, 375:931–950, March 2007.

2016A   A...594A..13P
Planck Collaboration, P. A. R. Ade, N. Aghanim, M. Arnaud,
  M. Ashdown, J. Aumont, C. Baccigalupi, A. J. Banday, R. B.
  Barreiro, J. G. Bartlett, and et al.
Planck 2015 results. XIII. Cosmological parameters.
A&A, 594:A13, September 2016.

1990MNRAS.242..318S
W. Saunders, M. Rowan-Robinson, A. Lawrence, G. Efstathiou,
  N. Kaiser, R. S. Ellis, and C. S. Frenk.
The 60-micron and far-infrared luminosity functions of IRAS
  galaxies.
MNRAS, 242:318–337, January 1990.

1989ApJ...338...13C
J. J. Condon.
The 1.4 gigahertz luminosity function and its evolution.
ApJ, 338:13–23, March 1989.


















§ COMPUTATION OF THE CONFUSION LIMIT



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

    B_ circ(θ,φ) = exp{ -4 ln2 θ^2/θ_0^2}  ,   θ_0 = 1.029 λ/D

where θ_0 is the FWHM resolution, λ is the observation wavelength,
D is the aperture size of the telescope, 
or an elliptic Gaussian function 

    B_ ell(θ, φ) = exp{ -4ln2 θ^2 [ (cosφ/θ_1)^2 + 
    			(sinφ/θ_2)^2 ] }


For a telescope with a beam response B(θ,φ), a source with flux S will induce an output  
x=S · B(θ,φ), the average number of response R(x) with an amplitude
between x and x+ dx in a solid angle Ω_b is then
<cit.>

    R(x) = ∫_Ω_bn( x/B(θ,φ))/B(θ,φ) dΩ

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

    R(x) = k  x^-γ π θ_0^2/(γ-1) 4 ln2

Let S_c=Qσ_c where Q is a factor of 3–5 <cit.>, the rms confusion σ_c is 
<cit.>

    σ_c^2 = ∫_0^S_c x^2 R(x)  dx

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

    S_c = [ Q^2  k π θ_0^2/(3-γ)(γ-1)  4ln2]^1/γ-1
	
For broader range of frequencies, we may model  n(S) using a piece-wise power law <cit.>, 

    S^5/2 n(S) = k_n  S^5/2 - γ__n,    n=1,2,3,4,5,⋯

The best-fit parameters are 


    [ γ_1=1.57,    k_1=29500                    for         S ≤ 10^-4   Jy;    γ_2=2.23,    k_2=60                    for  10^-4< S ≤ 10^-3   Jy;  γ_3=1.77,    k_3=1300                    for  10^-3< S ≤ 10^-1   Jy;   γ_4=2.20,    k_4=440                    for 10^-1< S ≤ 1        Jy;   γ_5=2.77,    k_5=390                    for     S    > 1        Jy ]


With these parameters we obtain the confusion limit, we find

    S_c = {[ 0.086   mJy· Q^3.51(θ_0/ arcmin)^3.51(ν/ GHz)^-2α ,                                         θ_0 ≤ 0.25';  0.12   mJy· Q^1.53(θ_0/ arcmin)^1.53(ν/ GHz)^-2α ,                                 0.25' < θ_0 ≤ 1.13'; 0.017   mJy· Q^2.56(θ_0/ arcmin)^2.56(ν/ GHz)^-2α ,                                 1.13' < θ_0 ≤ 7.18';  0.51   mJy· Q^1.63(θ_0/ arcmin)^1.63(ν/ GHz)^-2α ,                                   7.18' < θ_0 ≤ 30';     5   mJy· Q^1.18(θ_0/ arcmin)^1.18(ν/ GHz)^-2α ,                                           θ_0 > 30' ].

where α is the spectrum index of flux for the sources, typically α= 0.7–0.8 (Fig. <ref>), 
For example, for the NVSS+SUMSS catalogue, θ_0=45” =0.75', if choose Q=5, 
the confusion limit S_c=0.5 mJy and the rms σ_c = 0.1 mJy. So
point source whose flux density above 0.5 mJy at 1.4 GHz can be distinguished.













§ 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 
<cit.> and the STY maximum-likelihood method
<cit.>. 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: 

    ρ_m(L) = ∑_i=1^N ( 1/V_m)_i

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

    σ = √(∑_i=1^N ( 1/V_m)_i^2 )

We use the data from <cit.>, with the cosmological parameters from <cit.>. 

For start-forming galaxies, the luminosity function is of the form<cit.>

    ρ_m(L) = C [ L/L_*]^1-aexp{ -1/2[ log_10( 1+L/L_*)/σ]^2 }

The best-fit parameters of Eq. (<ref>) are

    C          =     10^-2.88 ± 0.74  mag^-1 Mpc^-3
    
    	L_*    =     10^21.22 ± 5.67   W  Hz^-1
    
    	a          =     1.04 ± 0.21                                                  
    σ    =     0.61 ± 0.14


For radio-loud AGNs, the luminosity function is of the form<cit.>

    log_10(ρ_m) = Y - 3/2log_10(L) - √( B^2 + [ log_10(L) - X/W]^2 )

The best-fit parameters of Eq. (<ref>) are

    Y     =     33.89 ± 7.34  
    
    	B     =      2.27 ± 0.83  
    
    	X     =     25.96 ± 2.27            
    
    	W     =      0.85 ± 0.12




<ref> shows the fitted local radio luminosity function.


For an isotropic source with spectrum index α≡- d[ln(S)]/ d[ln(ν)] at redshift z, its flux is given by 

    S = L/4π D_L^2 (1+z)^1+α = L/A (1+z)^1+α

where L is the luminosity, and the luminosity distance of this source is
D_L = c(1+z)/H_0∫_0^z  dz^'/√(Ω_Λ + Ω_m(1+z^')^3),
and we adopt the ΛCDM cosmological model with the parameters from Plank 2015 <cit.>.
Then we have

    dL/ dS = A(1+z)^1+α


In a comoving distance r to r+ dr, the number of sources with luminosities L to L+ dL is 
ρ_m(L,z)A dL dr. 
This number can also be expressed by the sources with flux densities S to S+ dS in redshift range z to z+ dz,  
that is η(S,z) dS dz. The comoving distance is given by d r = c dz/H(z).
Let η(S,z) dS dz = ρ_m(L,z)A dL dr, we obtain

    η(S,z) = ρ(L,z) A  dL/ dS dr/ dz = cA^2 (1+z)^1+αρ(L,z)/H_0 √(Ω_Λ + Ω_m(1+z)^3)

	
The weighted differential count is 

    S^5/2 n(S) = 1/4π∫_0^∞ S^5/2η(S,z)  dz

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” <cit.>,

    ρ_m(L,z) = g(z) ρ_m ( L/f(z),0 )

where

    f(z) = (1+z)^x  ,    g(z) = (1+z)^y   e^-(z/z_c)^q
	
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. <ref>) is shown in Fig. <ref>, and
the best-fit parameters are

    x     =      6.68 ± 1.33  
    
    	  y     =     -8.34 ± 1.86  
    
    	z_c     =      2.39 ± 0.17            
    
    	  q     =      2.53 ± 0.24

If approximate as a single power law, we dind n(S) ≈ 1300  S^-1.77. 














