Impact of extragalactic point sources on the low-frequency sky spectrum and cosmic dawn global 21-cm measurements

Contribution of resolved and unresolved extragalactic point sources to the low-frequency sky spectrum is a potentially non-negligible part of the astrophysical foregrounds for cosmic dawn 21-cm experiments. The clustering of such point sources on the sky, combined with the frequency-dependence of the antenna beam, can also make this contribution chromatic. By combining low-frequency measurements of the luminosity function and the angular correlation function of extragalactic point sources, we develop a model for the contribution of these sources to the low-frequency sky spectrum. Using this model, we find that the contribution of sources with flux density $>10^{-6}\,$Jy to the sky-averaged spectrum is smooth and of the order of a few kelvins at 50--$200\,$MHz. We combine this model with measurements of the galactic foreground spectrum and weigh the resultant sky by the beam directivity of the conical log-spiral antenna planned as part of the Radio Experiment for the Analysis of Cosmic Hydrogen (REACH) project. We find that the contribution of point sources to the resultant spectrum is $\sim0.4\%$ of the total foregrounds, but still larger by at least an order of magnitude than the standard predictions for the cosmological 21-cm signal. As a result, not accounting for the point-source contribution leads to a systematic bias in 21-cm signal recovery. We show, however, that in the REACH case, this reconstruction bias can be removed by modelling the point-source contribution as a power law with a running spectral index. We make our code publicly available as a Python package labelled epspy.


INTRODUCTION
The cosmological 21-cm signal holds a wealth of information about the dark ages, cosmic dawn and epoch of reionisation (Furlanetto et al. 2006;Pritchard & Loeb 2012;Liu & Shaw 2020;Shimabukuro et al. 2022;Bera et al. 2023).It is wellknown, however, that the detection of this signal is extremely challenging due to a number of reasons, such as the rapidly changing ionosphere (Datta et al. 2016;Jordan et al. 2017;Shen et al. 2021), radio frequency interference, noise, and systematic errors (Offringa et al. 2015;Scheutwinkel et al. 2022;Leeney et al. 2023;Anstey & Leeney 2024), chromatic response of the radio antenna (Vedantham et al. 2014;Mozdzen et al. 2016;Wang et al. 2024), and most importantly the astrophysical foregrounds (Oh & Mack 2003;Gnedin & Shaver 2004;Santos et al. 2005;Bowman et al. 2006;Jelić et al. 2008).While the strength of the cosmological signal of interest is expected to be of the order of ∼ 100 mK, the astrophysical foregrounds can easily be 10s or 100s of kelvins.These foregrounds are composed of galactic and extragalactic contributions.Primarily, synchrotron radiation due to accel-⋆ E-mail: shikhar.mittal@tifr.res.in,sm2941@cam.ac.uk erating relativistic electrons in the galactic magnetic field and free-free emission from the ionized hydrogen gives rise to diffuse galactic emissions (Condon 1992;Shaver et al. 1999), which can be ∼ 75% of the total foregrounds (Bernardi et al. 2009) at 150 MHz.The remaining part is mostly due to discrete extragalactic radio sources, only some of which are resolved by radio surveys.
In this work we model the contribution of the extragalactic radio point sources to the foregrounds for cosmic dawn 21-cm experiments.Such contamination by extragalactic point sources has previously been studied in the context of CMB experiments (González-Nuevo et al. 2005, and references therein).We also investigate the effect that point sources have on the detectability of cosmic dawn global 21-cm signal.In this work, we focus on the case of the Radio Experiment for the Analysis of Cosmic Hydrogen (REACH ) 1 (de Lera Acedo 2019; de Lera Acedo et al. 2022) experiment.Nonetheless, our methods are also applicable to other global 21-cm experiments such as Experiment to Detect the Global EoR Signal (EDGES, Bowman & Rogers 2010), Shaped Antenna measurement of the background RAdio Spectrum (SARAS, Patra et al. 2013;Singh et al. 2022), Large Aperture Experiment to Detect the Dark Ages (LEDA, Price et al. 2018), Probing Radio Intensity at high-Z from Marion (PRIzM, Philip et al. 2019), and Mapper of the IGM Spin Temperature (MIST, Monsalve et al. 2024), as well as to experiments targeting the dark ages 21-cm signal by means of lunar-based telescopes such as Dark Ages Radio Explorer (DARE, Burns et al. 2012Burns et al. , 2017)).
This paper is organized as follows.In section 2, we describe our empirical model of extragalactic point sources.In section 3 we investigate the contribution of these sources to the low-frequency sky spectrum as seen by REACH.In section 4 we discuss the implications of point sources contribution on the 21-cm signal extraction.We summarise our conclusions in section 5.

A MODEL FOR POINT SOURCES
Our aim is to model the contribution Tps to the sky spectrum by extragalactic point radio sources as a function of position on the sky, n, and frequency, ν.Thus, we seek Tps = Tps (n, ν).We construct a model for this by incorporating three observational inputs, namely, i) measurements of the flux density distribution, dn/dS, ii) measurements of the clustering of point sources on the sky, C(χ), and iii) the spectral energy distribution, S = S(ν), of the sources.

Flux density distribution
Several measurements of the flux density distribution of point sources exist in literature at different survey frequencies and sky coverage.Baldwin et al. (1985), Hales et al. (1988), andMcGilchrist et al. (1990) obtained some of the earliest differential source counts based on radio catalogues by 6C and 7C surveys at 151 MHz with the faintest source resolved at 0.1 Jy.While some modern surveys have a similar sensitivity to the older survey, a few of the new generation lowfrequency telescopes are uncovering radio sky with unprecedented depth and sensitivity.Intema et al. (2017)  In this work we follow the fitting function for dn/dS reported by Gervasi et al. (2008) based on 150-MHz 6C and 7C surveys  The blue curve is a sum of two inverse double-power-laws (see equation 1).The green and red points are obtained from Mandal et al. (2021) and Franzen et al. (2019), respectively.The uncertainty bars are too small to be visible in this figure .where the 8 numbers, A1, B1, . .., are empirical fitting parameters summarised in Table 1 and S and dn/dS are in units of Jy and Jy −1 sr −1 , respectively.See blue curve in Fig. 1.As evident the fitting function (1) fits reasonably to the newer, TGSS-ADR1 and GLEAM data points, which are shown in light green squares and light red circles, respectively.The uncertainty bars are too small to be visible on the figure.
The fitting function (1) implies that the total flux density converges to a finite value as S approaches 0 but the number of sources increase indefinitely.Since we need to scale the flux density for individual sources to a range of frequencies we avoid modelling of infinite number of sources -without sacrificing the accuracy -by setting the lower limit of flux density of the point sources to Smin = 10 −6 Jy (Di Matteo et al. 2002;Wang et al. 2006;Liu et al. 2009).
Previous works investigated impact of only the unresolved point sources.Accordingly, they chose cut-off value above which the point sources can be resolved.For example, Wang et al. (2006) and Liu & Tegmark (2011) adopt Smax to be 10 −4 and 0.1 Jy, respectively.In this work we do not make any distinction between unresolved and resolved point sources, and account for their contribution alike.For our choice of dn/dS it is extremely rare to find sources above 10 2 Jy.Stated differently, extending our S range beyond 10 2 Jy does not impact the collective foregrounds by more than a few percent (Gervasi et al. 2008).Hence, we chose Smax = 10 2 Jy.We do not consider the effect of brighter sources here.

Clustering
At zeroth level one expects the brightness temperature contributed by the point sources to be isotropic with the sources having a random distribution and fluctuation in the number count to be Poissonian.However, a closer look suggests that the anisotropies in the large-scale structure of the Universe will be imprinted on the distribution of sources (Peebles 1993;Blake et al. 2004;Wake et al. 2008).Consequently, we account for the clustering of point sources.
The anisotropy in the point source positioning, or simply the clustering, can be quantified by the 2-point angular correlation function (2PACF) of the fluctuation of the overdensity δps = δps(n), i.e., where nps is the mean number of point sources per pixel and nps = nps(n) is the number of point sources per pixel on the clustered sky.Here n represents the coordinates of a pixel on the sky.Under the assumption of mean isotropy of the Universe, the 2PACF depends only the relative positions represented by n′ and n.Thus, where χ is the angle between n′ and n (given by cos χ = n′ •n) (Peebles 1993;Peacock 2010).
A number of 2PACF have been derived at different frequencies and for different flux density range.For a summary of existing and ongoing surveys reporting the clustering see Rana & Bagla (2019).See also the latest work by Hale et al. (2024) who derive 2PACF based on LoTSS data at 144 MHz.It is evident from the inspection of results that the clustering depends on frequency of observation.Moreover, clustering law is expected to change in different flux density ranges as the brighter sources -possibly residing in heavier dark matter haloes -will cluster more strongly than the fainter sources (Cress et al. 1996;Overzier et al. 2003).As we employ the dn/dS function at 150 MHz we adopt the 2PACF at the same frequency for consistency.Rana & Bagla (2019) derive the 2PACF based on the TGSS-ADR1 survey at ν0 = 150 MHz.For the threshold flux density of 50 mJy the 2PACF is where γ = 0.821, A = 7.8 × 10 −3 and χ is in degrees.In this work we assume the same 2PACF is applicable for our choice of flux density range, which goes from Smin = 1 µJy to Smax = 10 2 Jy.As we show later, the choice of 2PACF has a negligible impact on our sky-averaged temperature.Thus, the assumption of a uniform 2PACF across our choice of flux density range is an excellent approximation.

Spectral energy distribution
Radio experiments will observe the sky for the global signal at a range of frequencies.We thus need to know how the flux density of a source evolves with frequency.The extragalactic radio sources are believed to be active galactic nuclei (AGNs), radio galaxies and relics, star-forming galaxies and haloes that give free-free emission and intergalactic medium (Gleser et al. 2008;Niţu et al. 2021).Among these, AGNs dominate the total extragalactic foregrounds.Unlike in the case of galactic foregrounds, range of emission mechanisms for different types of extragalactic sources makes their collective flux density a complicated function of frequency.As we will see later, a power law function with a running index explains the collective brightness temperature (even the contribution to antenna temperature for a conical log-spiral and a hexagonal dipole antenna) spectrum sufficiently well.See Fig. 7 and Table 3.
At an individual source level, we assume that the spectral energy distribution (SED) for all types of extragalactic sources is very-well fit by a power law for frequencies of our interest, so that S(ν) ∝ ν −α .However, note that besides these standard non-thermal power-law models, there are other SED models in radio band proposed in literature such as the curved power-law model (Callingham et al. 2017).We investigate such models in a future work.
In terms of brightness temperature, S(ν) ∝ ν −α translates to T (ν) ∝ ν −β , where β = α + 2. Often the radio sources are categorised into flat-spectrum (β < 2.5) and steep-spectrum (β > 2.5).Steep-spectrum are more common at low frequencies -such as the frequencies of interest in this work -while flat-spectrum sources are more common at high frequencies, typically ν ≳ 2 GHz (Peacock & Gull 1981).In order to capture the wide spread in β values we follow the strategy by Liu et al. (2009); we assume the indices to be normally distributed around β0, i.e., where we take β0 = 2.681, inferred by Gervasi et al. (2008) and a spread of σ β = 0.5 (Tegmark et al. 2000).A Gaussian distribution of indices is motivated by TGSS-ADR1 and NVSS surveys data (Intema et al. 2017;Tiwari 2019).Our chosen β0 value is consistent with steep-spectrum which dominate the point source population at low frequencies.
We summarise our choice of free parameters in Table 2.

Modelling point sources
In this section we put together the 3 observational inputs to simulate the foregrounds contributed by the point sources.
There are three main steps to simulate point sources.In brief they are as follows: 1) Find the total number of point sources given an S distribution 2) Distribute these sources on the sky based on the clustering law 3) Pixel-wise compute the flux density at any frequency and convert it to brightness temperature We now go over each step in detail.For our first step we find the total number of point sources on the whole sky as which comes out to be ≈ 3.26 × 10 5 Jy.As mentioned previously, a lower Smin and a higher Smax than our current choice changes this number by ∼ 1%.
Next we distribute Nps sources on the sky given the 2PACF, equation ( 4).We begin by converting 2PACF to angular power spectrum, C ℓ .The following is the conversion relation where P ℓ are the Legendre polynomials.We use the python package transformcl2 (Tessore et al. 2023) for this transformation.Note that the angular power spectrum (APS) obtained in the above equation is dimensionless.
With the APS so obtained we find the corresponding map δps using the synfast algorithm from HEALPix (Górski et al. 2005).We work with nside = 2 9 so that the number of pixels is Npix = 3145728.Finally, the point source distribution in terms of number of sources per pixel is where nps = Nps/Npix, which for our chosen pixelisation and flux density range comes out to be 1402.
Figure 2 shows the number density (number per pixel) of the extragalactic point sources.The top panel shows an isotropic sky (Poissonian) and the bottom panel shows a clustered sky following the law (4)3 .The contrast between the two skies is immediately apparent by a visual comparison.The sky in the bottom panel shows patchiness which is the result of the clustering of the sources.For a quantitative conclusion one might compare the APS of the overdensity for the two skies.These are shown in Fig. 3.As evident, for an isotropic sky APS is a horizontal line and 0 (albeit with some white noise) but not for the clustered sky.
Having populated our sky, we assign flux density and spectral index to each of the Nps sources.Let nps,i be the number of sources on the i th pixel or a patch on the sky.We draw nps,i flux densities (they will be at the reference frequency ν0) from the S distribution (equation 1).We convert the flux density to brightness temperature via the standard Rayleigh-Jeans limit of Blackbody function.Putting together everything, the brightness temperature at the i th pixel due to j th point source is   (n).Dotted red curve is for the isotropic sky and blue dashed is for a clustered sky (corresponding to the 2-point angular correlation function in equation 4), which also serves as the input for synfast.The green solid curve shows the power spectrum obtained for the simulated map of overdensity, in agreement with the input power spectrum.The corresponding maps of number density are shown in Fig. 2.
where k B is the Boltzmann constant, Sr is a random flux density drawn from the distribution and Ωpix = 4π/Npix ≈ 4 × 10 −6 sr (this corresponds to a resolution of ≈ 0.11 • ) is the solid angle subtended by each pixel.
Next we draw nps,i spectral indices from the β distribution (equation 5).Finally, translate the individual brightness temperatures to a frequency ν (according to the chosen β) and sum them to get the brightness temperature on the i th pixel, i.e., Tps,i(ν), We do the above exercise for all pixels i = 1, 2, . . .Npix, and a range of frequencies ν to obtain the brightness temperature map, Tps(n, ν). Figure 4 shows an example of point sources brightness temperature at a frequency of 150 MHz which corresponds to a redshift of z = 8.4.For our realisation, the global average is 1.28 K.In Fig. 5 the green-dotted curve shows the global average of the point sources brightness temperature in the frequency range 50 to 200 MHz.We represent the global average of the point sources as ⟨Tps⟩ which can be calculated as The curve seen in Fig. 5 can be described by a power law of index −2.68(additionally with a small running index).
It would be interesting to consider extragalactic point sources as a potential source of the excess radio background (ERB) reported by ARCADE2 /LWA1 experiments (Fixsen et al. 2011;Dowell & Taylor 2018) partly or wholly.(For the latest report on ERB see Singal et al. (2023)).The spectrum of ERB has a spectral index of −2.58 ± 0.05, which is close to what point sources imply.This is not a surprising result though, since AGNs dominate our extragalactic source population and it is well-known that the accretion-induced emission from AGNs or AGN-like objects have a strong radio emission which have power law spectrum with index ≈ −2.6 (Ewall- Wice et al. 2018Wice et al. , 2019;;Mittal & Kulkarni 2022a).Todarello et al. (2024) investigate extragalactic point sources and ERB and concluded that about 20% of ERB must be of extragalactic origin that traces the large-scale structure.
To complete the model for a simulated sky data we account for the galactic emissions.Additionally, for testing our inference pipeline we also inject to this a Gaussian 21-cm signal.Section 4 gives these details.
Note that for results shown in figures 4 and 5 the effect of chromaticity has not yet been taken into account.This is discussed in the next section.

CONTRIBUTION TO THE REACH ANTENNA TEMPERATURE
Antenna chromaticity arises because of the dependence of sensitivity of the antenna to the signal on the frequency and direction of observation.As a result of antenna beam chromaticity, the 'effective' brightness temperature, or commonly known as the antenna temperature, registered on the antenna in the direction n and frequency ν is given by weighting the sky map with directivity pattern of antenna, D(n, ν), so that T A (n, ν) ∼ D(n, ν)Ttot(n, ν).Thus, the spectrum of sky as seen by the antenna or simply the beam-weighted foregrounds is given by (Anstey et al. 2021, hereafter A21) where we also add an uncorrelated Gaussian noise with a standard deviation of 0.025 K representing the antenna noise, σ A .Section 4 gives the details of galactic foregrounds and the total sky temperature, Ttot = Ttot(n, ν).
For our beam directivity pattern we assume a conical logspiral antenna as appropriate for REACH (Cumner et al. 2022).We do not work with time varying data because of which we do not have a time integral.Stated differently, we work for a fixed snapshot of the Global Sky Model which we choose to be at UTC 0 h: 0 m: 0 s 1 st January 2019 at which the Galaxy is above the horizon.The antenna location, height above sea level and orientation (angle between antenna's xaxis and North) are (30.71• S, 21.45 • E), 1151 m and 0 • , respectively (de Lera Acedo et al. 2022).
Our focus of investigation in this work will be on the point sources contribution to T A , which may be evaluated as Figure 6 shows the total antenna temperature for a conical log-spiral antenna given the total sky brightness temperature map, Ttot(n, ν).Here we have used the fiducial model parameters for point sources given in Table 2.The red crosses in the top panel show the antenna temperature which does not include the point sources and the solid black curve shows the antenna temperature that includes the point sources emission.The difference between the two, represented as T A,ps , goes from 0.6 to 27.3 K.The total antenna temperature (with point sources contribution) goes from 151 to 6851 K. Thus, the point source contribution is ∼ 0.4 percent of the total antenna temperature throughout the frequency range 50 to 200 MHz.
In Fig. 7 we show the difference, T A,ps , for various properties of point sources.The thick solid blue curve is repeated from the bottom panel of Fig. 6.Models which differ only in β0 and σ β with respect to the fiducial model will all have the same temperature at the reference frequency (150 MHz in our case), which is why the blue, light blue, cyan, green and orange curves pass through the same point.This temperature is ∼ 1.3 K.However, this temperature will scale differently with frequency and hence the shape is different.Allowing sources to have a higher S hardly makes difference.This is because such sources are very rare to find given the S distribution and thus the collective flux remains close to the fiducial value.Increasing Smin to a value of, say, 10 −4 Jy though results in a smaller number of sources (∼ 1.6 × 10 8 ), sources with higher S are relatively more likely to be found.This results in stronger extragalactic foregrounds, and hence the antenna temperature, which is why the red curve is above the fiducial model curve.Finally, the variation in clustering law has the least impact on the antenna temperature; pink and magenta curves nearly overlap the fiducial model curve.This is because the different number density distributions we consider for different 2PACFs are only mildly different from each other.As a result we do not see much difference in Tps  maps.Even in the worst case scenario if the beam directivity was a delta function peaked on the pixel of least number of sources (for fiducial model we have ∼ 500; see Fig. 2), number of sources would still be large enough so as to have the same functional form for the temperature as that for the whole sky.Thus, we get nearly the same trend for T A,ps when we change 2PACF.
We find that T A,ps = T A,ps (ν) can be described by a powerlaw function with a running spectral index as follows where we fix ν0 to 150 MHz.We fit this form to the different spectra seen in Fig. 7; we find an excellent fit to the data.We have not shown the best-fitting curves as they perfectly overlap, however, Table 3 gives the best-fitting parameter values T f , β f and ∆β f .The 1σ uncertainty on all the reported numbers is of the order of 10 −5 .From Fig. 7 and Table 3 it is evident that the flux density range controls the overall amplitude while the spectral index distribution of the radio sources controls the shape of the antenna temperature spectrum.The clustering of sources has a negligible impact.Another interesting result is to compare the point sources spectrum and its contribution to antenna temperature thus explicitly bringing out the chromatic distortions.This is shown in Fig. 8 for the fiducial model.The blue solid curve is T A,ps (which is repeated from the bottom panel of Fig. 6) and the green dotted curve (repeated from Fig. 5) shows sky average of the point sources, which can also be thought of as T A,ps for a perfectly achromatic antenna for the full sky, i.e., with D(n, ν) = 1 for all n and ν.The red curve in the bottom panel shows the difference T A,ps −⟨Tps⟩.Maximum chromatic distortion is of the order of 100 mK seen at 50 MHz.

BIAS IN SIGNAL RECONSTRUCTION DUE TO POINT SOURCES
So far in this paper we have developed a model for the distribution of extragalactic point sources on the sky.We have also studied the contribution of these sources to the REACH beam.We now investigate its effect on the 21-cm signal reconstruction.To do this, we consider a simulated data set which includes the point sources, galactic emissions and a Gaussian 21-cm signal.From this mock dataset, we attempt to extract the 21-cm signal using a Bayesian inference pipeline.We will show that the presence of point sources can bias the 21-cm signal reconstruction.The strongest component in the low-frequency spectrum measured by a global 21-cm experiment is due to the galactic emissions.These can be modelled using the Global Sky Model (GSM) maps (de Oliveira-Costa et al. 2008;Zheng et al. 2017).These maps have been de-sourced for bright extragalactic sources such as giant elliptical galaxies, radio galaxies and quasars.Using GSM, we construct the pure galactic emission map as (A21), where β gal is the spectral index for the galactic emissions.The value of β gal on each pixel of the sky can be computed given Tgsm(n, ν) at two reference frequencies 230 MHz and  ( We emphasize that galactic emission map T gal does not include the CMB. In addition to the foregrounds we have a constant backdrop of CMB, T cmb = 2.73 K.And finally, we add the 21-cm signal.Literature on theoretical modelling of 21-cm signal suggests that at cosmic dawn the signal shape is Gaussian-like (Mirocha et al. 2013;Fialkov et al. 2014;Mittal & Kulkarni 2021;Mittal et al. 2022).While there exist theoretical models of the global 21-cm signal in which the absorption feature has a non-Gaussian shape (Mittal & Kulkarni 2022b), here we adopt a Gaussian form given by where we set the amplitude to T 0 21 = 0.155 K, the central frequency to νc = 85 MHz and the Gaussian width to σ21 = 15 MHz.
Weighing the total sky temperature, Ttot, by the telescope beam directivity will result in our mock observation of the sky spectrum.(Note that only the galactic and point sources component are spatially varying, while CMB and the 21-cm signal are functions of ν only.)The goal of global 21-cm signal reconstruction is to extract the 21-cm signal from the observed sky spectrum.
Given the antenna temperature data we attempt to extract the 21-cm signal.For our inference we follow a fully Bayesian analysis pipeline.In Bayesian statistics, a model M characterised by parameter set θ is fit to a data set D for a given choice of priors.Using Bayes' theorem the posterior distribution, or the probability of getting a certain parameter set can be computed as P(θ|D, M) = P(D|θ, M)P(θ|M) where P(D|θ, M) = L is the likelihood, P(θ|M) = π is the prior distribution and P(D|M) = Z is the Bayesian evidence.
We work with a Gaussian likelihood and write where σ ′ is a uniform noise value across the entire frequency band (Scheutwinkel et al. 2023).The sum runs over 151 values corresponding to frequencies ranging from 50 to 200 MHz on interval of 1 MHz.Note that σ ′ is a free parameter and is therefore part of θ.We implement our Bayesian pipeline using the python package PolyChord4 (Handley et al. 2015a,b), using the code's default settings throughout.

Reconstruction in the absence of point sources
We first consider the scenario when there is no point-source contribution to the sky.In this case we have Weighing Ttot by the beam directivity using equation ( 13) gives us the mock data, TD.
For our Bayesian inference pipeline we require an antenna temperature model, TM.A21 showed that non-smooth distortions introduced by chromatic distortions cannot be captured by smooth polynomial foreground fits.Thus, in this work we work with the parametrized sky model, in the same spirit of simulated antenna temperature generation, for the inference procedure.We construct a fitting model for the galactic foregrounds map similar to the construction used for the simulated data; we scale the brightness temperature at a known reference or a base frequency to the required frequency using the spectral index map.These spectral indices will be the free parameters which need to be inferred.However, for a faster likelihood inference while maintaining accuracy, we work with a course-grained version of spectral index map with effectively only Nreg regions or pixels on the sky with different β's.Throughout we work with Nreg = 9.Thus, due to galactic component we introduce 9 free parameters.We hereafter refer to the course-grained galactic foregrounds model by A21 as the 'N -region' model.
Just as in the case of simulated data construction, we next add the constant CMB piece to the model.This does not have any associated free parameter.
Finally, we assume a Gaussian 21-cm signal characterised by 3 parameters, namely, amplitude, width and the central frequency.Thus, we have our first model to fit to the data TM = (beam directivity) × (N -region galactic + CMB+ Gaussian 21-cm signal) , where, as in the case of simulated data generation, we have weighed the net temperature by the perfectly-known beam directivity.(It might be possible to consider the case of a parametrized beam directivity but we do not explore it in this work.) Thus, for the inference we have 9 parameters for the foreground model, 3 parameters for the Gaussian 21-cm model and 1 additional parameter which is the uncorrelated antenna noise, making a total of 13 parameters.For all parameters we use uniform priors.
The scenario considered here is a reproduction of the results from A21 (which is the default REACH pipeline).In this case we only have galactic emission as part of the foregrounds.Consequently, for inference procedure a course-grained Nregion foregrounds model suffices as already established by A21. Figure 9 shows the recovered or the best-fitting 21-cm signal and its posterior in blue, the residuals after subtraction of N -region model of foregrounds in red and the true injected signal in green.As evident the recovered signal is in excellent agreement with true injected signal.

Mock data with point sources
We now consider the scenario when we have a finite point sources contribution to the total sky brightness.In this case the total sky data is given by As before, weighing Ttot by the beam directivity, D, gives us the mock antenna temperature data, TD.The point sources contribution that we consider above follows the fiducial set of parameters given in Table 2.Note that this contribution to antenna temperature data can be described by the power-law function with a running spectral index with parameter values T = 1.234K, β = 2.680 and ∆β = 0.126 (see 'Fiducial' row in Table 3).Just as in the case of 21-cm signal, these values serve as the true 'injected' parameter values allowing for a pipeline testing, as we demonstrate in subsection 4.2.2.The Ttot constructed above contains galactic as well as extragalactic emission, CMB and 21-cm signal.The brightness temperature data so obtained is a result of a variety of astrophysical and cosmological processes.It is a simple yet sufficiently realistic model representative of the real picture considered for the first time.
For fitting procedure when we have point sources contribution to the data, we consider two sub-cases as described below.

Signal recovery without correcting for point sources
For our first sub-case we use exactly the same fitting model as in the previous section so that our model is TM = (beam directivity) × (N -region galactic + CMB+ Gaussian 21-cm signal) .
Thus, we have 13 free parameters and we set the same uniform prior ranges as in the previous case.
As evident from Fig. 10, the signal recovery is quite poor in this case.While the data has four components, namely a galactic, an extragalactic, CMB and a 21-cm signal, fitting model used in this pipeline has all but the extragalactic part.Thus, there is no piece in the fitting model to account for the extragalactic foregrounds and as a result the residuals after removing best-fitting foregrounds from simulated data are quite large, as shown by red posteriors.Also, note that value of noise parameter inferred is σ ′ = (0.1 ± 4.7 × 10 −5 ) K which is of the order of the 21-cm signal strength.
The inferred values for the 21-cm signal are νc = 82.09± 0.81 MHz, σ21 = 10.50 ± 0.42 MHz and T 0 21 = 0.248 ± 0.002 K, which are more than 3-sigma away from the true value, especially the signal depth T 0 21 .

Improved reconstruction of signal with correction for point sources
For the final case, in our modelling we account for the extragalactic point sources contribution in the data.As we have demonstrated via Fig. 7 and Table 3, the point sources contribution to the antenna temperature is a smooth function of frequency -a power-law-with-a-running-index function.Motivated by this we propose equation ( 15) to account for the point sources contribution present in the antenna temperature data.Thus, TM = (beam directivity) × (N -region galactic+ power-law-with-running-index extragalactic + CMB+ Gaussian 21-cm signal) .
Correspondingly, we have three additional free parameters T f , β f and ∆β f .So now we have 9 + 3 = 12 parameters for the foreground model, 3 parameters for the Gaussian 21-cm model and 1 additional parameter which is the uncorrelated  9, the mock data contains the extragalactic point sources contribution.However, the fitting model remains the same having a galactic term, CMB and 21-cm signal.As evident we do not recover the 21-cm signal.
antenna noise, making a total of 16 parameters.For all parameters we use uniform priors, either in linear scale or log scale.We have a log uniform prior [10 −3 , 10 2 ] K for T f , uniform prior [0, 5] for β f and uniform prior [0, 5] for ∆β f .
Note that, as opposed to spectral inhomogeneity in the Nregion model to capture the galactic emissions in the data, there is no spectral inhomogeneity associated with the point sources model.Simply stated there is no n dependence.
We show our results for this setup in Fig. 11.As evident the result is quite promising showing that the net sky model, which includes point sources contribution, can be fit reliably using a course-grained N -region in conjunction with a 3-parameter power-law-with-a-running-index foregrounds model.The inferred values of point source model parameters are T f = 1.072 ± 0.041 K, β f = 2.333 ± 0.102 and ∆β f = 0.067 ± 0.045.Thus, the true 'injected' values of the point sources model are within 3-sigma of the inferred values.(Recall that the true 'injected' values are T = 1.234K, β = 2.680 and ∆β = 0.126).Also, the true 21-cm signal parameters are within 1-sigma of the inferred values, νc = 85.38 ± 1.61 MHz, σ21 = 14.47 ± 1.68 MHz and T 0 21 = 0.192 ± 0.023 K.The function power-law-with-a-running-index (equation 15) that we introduce in this work, when expressed in log-log units takes the familiar 'polynomial' (in log T -log ν space) form (precisely a second degree polynomial).Global 21-cm experiments such as EDGES (Bowman et al. 2018) and SARAS (Singh et al. 2022) have used a polynomial function to model the total foregrounds.However, Anstey et al. (2021) have shown that using a polynomial, in general, for the foregrounds (only galactic) can mask the 21-cm signal completely.Nevertheless, in this work it appears that extragalactic point sources contribution to foregrounds admits a polynomial description leading to a reliable inference of cosmological global 21-cm signal, as evidenced by our Fig.11.We emphasize that we use a polynomial only to model the extragalactic foreground component.10, we have an extragalactic contribution to the data.However, the fitting model has an extragalactic foregrounds model in addition to the galactic term, CMB and 21-cm signal.As evident we find that the signal recovery is excellent compared to the scenario shown in Fig. 10.
the galactic and extragalactic component and allows for an excellent signal recovery.
It remains to be investigated how the additional presence of ionosphere will affect our extraction of the signal (Shen et al. 2022).Other important effects to account for are the horizon (Pattison et al. 2024), time-varying systematics (Kirkham et al. 2024), sinusoidal-like instrumental systematics (Scheutwinkel et al. 2022) and amplitude errors in radio maps (Pagano et al. 2024).We leave the study of impact of these effects on the signal extraction in the presence of extragalactic point source emission for a future work.

CONCLUSIONS
In this work, we built a model for the contribution of extragalactic point sources to the low-frequency sky spectrum.Using measurements of the luminosity function and the angular correlation function of point sources, our model can compute a full sky with clustered point sources.For a flux limit of Smin = 10 −6 Jy, we find the all-sky spectrum has approximately a power-law distribution with typical values of a few kelvins at 50-200 MHz.
Further, we combined our point source model with a model of the galactic foregrounds, CMB and 21-cm signal to simulate the total sky temperature.We then weighed this by the beam pattern of a conical log-spiral antenna to simulate the antenna temperature as measured by the REACH experiment.We find that, at a representative frequency of 85 MHz which corresponds to a redshift of z = 15.7,contribution to the antenna temperature by the extragalactic point sources is 5.9 K while the total antenna temperature is 1584 K. Thus, extragalactic contribution is less than a percent compared to the total sky temperature data, as it is throughout the frequency range of interest, 50 to 200 MHz.While it seems insignificant compared to the total brightness, extragalactic contribution is still more than an order of magnitude stronger than the standard cosmological 21-cm signal.
We also find that the antenna temperature spectrum corresponding to a conical log-spiral antenna (REACH beam) is quite smooth for a wide range of properties of the point sources.Nevertheless, the chromatic distortions induced can be as high as 0.1 K, which is of the order of 21-cm signal strength.For any of these models we noted that a power law function with a running spectral index provides an excellent fit, with uncertainty levels of the order of 10 −5 , to the point sources spectrum.If unaccounted, in the presence of point sources in the data, the 21-cm signal recovery suffers with severe systematic bias.This bias can be successfully removed by incorporating our point-source model in the reconstruction procedure.This work paves the way forward for a more accurate inference for the current and upcoming global and interferometric 21-cm experiments.
The currently deployed antenna for the REACH experiment is the hexagonal dipole antenna (HDA).Anstey et al. (2022) showed that conical log-spiral antenna is far superior at recovering the 21-cm signal than HDA for various foreground and signal models because HDA is much more chromatic than a conical log-spiral.(However, conical log-spiral is costlier and harder to build than HDA).Phase I of REACH experiment operates an HDA.It will, thus, be interesting to consider the antenna temperature as observed through an HDA.
In this appendix we present our results shown in figures 6 to 8 for an HDA.For this antenna we have beam directivity data for a frequency range of 50 -150 MHz.
Figure A1 shows the mock antenna temperature with and without the extragalactic point sources contribution.The bottom panel shows the difference between the two data.The minimum, maximum and mean of the difference in T A,ps  for a hexagonal dipole and a conical log-spiral antenna are (approximately) 60, 3 and 22 mK, respectively.Thus, for the same point sources on the sky an HDA will report a higher antenna temperature than a conical log-spiral antenna throughout the frequency range.
Figure A2 shows antenna temperatures corresponding to a hexagonal dipole antenna for different point sources properties.As was the case with conical log-spiral antenna, all curves conform to the same functional form which is a powerlaw-with-a-running-index.Table A1 shows the result of a least-squares fitting to these curves.The uncertainty on all the numbers is of the order of 10 −4 .
Figure A3 shows an explicit comparison between averaged sky temperature due to point sources (⟨Tps⟩) and the corresponding antenna temperature (T A,ps ).The maximum chromaticity induced by hexagonal dipole is about 150 mK while that for conical log-spiral antenna is about 100 mK.
obtained the flux density distribution based on the first alternative data release of the TIFR GMRT Sky Survey (TGSS ADR1 ) by the Giant Metrewave Radio Telescope (GMRT ) covering 90% of the sky at 150 MHz with the faintest source resolved having flux density of 0.1 Jy.Mandal et al. (2021) did the same based on deep LOFAR Two Meter Sky Survey -LoTSS Deep Fields which cover the entire northern sky at 150 MHz with the faintest source resolved having flux density of 2.2 × 10 −4 Jy.As per our knowledge LoTSS Deep Fields resolves the faintest, ∼ 0.1 mJy, of the sources till date at 150 MHz.Franzen et al. (2019) did the same based on GaLactic and Extragalactic All-sky MWA (GLEAM ) at several frequencies between 72 and 231 MHz which covers entire sky south of declination 30 • .See also the latest work by Hale et al. (2024) and Tompkins et al. (2023) who have compiled data from a number of surveys at various frequencies.

Figure 1 .
Figure 1.Distribution of flux density of the extragalactic point sources at our chosen reference frequency of ν 0 = 150 MHz.The blue curve is a sum of two inverse double-power-laws (see equation 1).The green and red points are obtained from Mandal et al. (2021) and Franzen et al. (2019), respectively.The uncertainty bars are too small to be visible in this figure.
for dn/dS we get Nps ≈ 4.4 × 10 9 .Thus, we have 4.4 × 10 9 point sources on the whole sky in the flux density range 10 −6 and 10 2 Jy at a frequency of ν0 = 150 MHz.The total flux density is given by

Figure 2 .
Figure 2. Top panel shows the number of point sources per pixel for an isotropic sky.In this case the distribution is Poissonian.(Poisson fluctuations are small and hence not visible for the chosen colour range.)Bottom panel shows the number of sources for a clustered sky.The average number of point sources per pixel is 1402 in either case.The colour bar is in linear scale.

Figure 3 .
Figure3.The angular power spectrum of the number density contrast δps(n).Dotted red curve is for the isotropic sky and blue dashed is for a clustered sky (corresponding to the 2-point angular correlation function in equation 4), which also serves as the input for synfast.The green solid curve shows the power spectrum obtained for the simulated map of overdensity, in agreement with the input power spectrum.The corresponding maps of number density are shown in Fig.2.

Figure 4 .a
Figure 4.The brightness temperature due to extragalactic point sources at ν 0 = 150 MHz.The global or the sky average is 1.23 K.The colour bar is in logarithmic scale.The maximum value observed is ∼ 88 K but since the majority of pixels have very low brightness, we set the colour bar maximum to 5 K for better visualisation.

Figure 5 .
Figure 5. Sky average of the brightness temperature due to the extragalactic point sources in the frequency range 50 to 200 MHz.

Figure 6 .SFigure 7 .
Figure6.The top panel shows the antenna temperature for two different models of foregrounds.Solid black curve includes and the curve with red crosses does not include the contribution of point sources to the foregrounds.We follow the beam directivity of a conical log-spiral antenna; the REACH case.Bottom panel shows the difference between the data represented by solid black and crossed red curves.

Figure 8 .
Figure8.The thick blue solid line is the point sources contribution to antenna temperature (repeated from the bottom panel of Fig.6).Curve with green circles shows the sky average of point sources (repeated from Fig.5).The bottom panel shows the difference between the two curves.We see the chromatic distortions induced by the antenna beam chromaticity.

Figure 9 .
Figure9.The red region shows the posterior with different sigma levels of the residual after removing a foregrounds model (only galactic) from the mock antenna temperature data.The blue region shows the the posterior with different sigma levels of the bestfitting 21-cm signal.Finally, the green curve shows the true 21-cm signal that was injected to the mock data.The mock data does not contain the extragalactic point sources contribution.The fitting model has a galactic emissions term, CMB and 21-cm signal.The signal recovery is excellent in this case.

Figure 10 .
Figure10.Contrary to the scenario shown in Fig.9, the mock data contains the extragalactic point sources contribution.However, the fitting model remains the same having a galactic term, CMB and 21-cm signal.As evident we do not recover the 21-cm signal.

Figure 11 .
Figure11.Similar to the case shown in Fig.10, we have an extragalactic contribution to the data.However, the fitting model has an extragalactic foregrounds model in addition to the galactic term, CMB and 21-cm signal.As evident we find that the signal recovery is excellent compared to the scenario shown in Fig.10.

Figure A1 .
Figure A1.Same as Fig. 6 but for a hexagonal dipole antenna.

Figure A2 .Figure A3 .
Figure A2.Same as Fig. 7 but for a hexagonal dipole antenna.The solid blue curve is repeated from Fig. A1.

Table 1 .
Gervasi et al. (2008)meters for the flux density distribution function at reference frequency of ∼ 150 MHz.Taken fromGervasi et al. (2008).We have adopted the weighted average values for a 1 , b 1 , a 2 , b 2 and B 2 /B 1 .

Table 3 .
Best-fitting parameters T f , β f and ∆β f for power-lawwith-a-running-index (equation 15) against the extragalactic component of the antenna temperature data using a least-squares fitting.In this table we refer to the various antenna temperature data models by the legend labels shown in Fig.7.The 1σ uncertainty on all the reported numbers is of the order of 10 −5 .

Table A1 .
Same as Table3but for a hexagonal dipole antenna.The uncertainty on all the reported numbers is of the order of 10 −4 .