Detecting Cosmic 21 cm Global Signal Using an Improved Polynomial Fitting Algorithm

Detecting the cosmic 21 cm signal from Epoch of Reionization (EoR) has always been a difficult task. Although the Galactic foreground can be regarded as a smooth power-law spectrum, due to the chromaticity of the antenna, additional structure will be introduced into the global spectrum, making the polynomial fitting algorithm perform poorly. In this paper, we introduce an improved polynomial fitting algorithm - the Vari-Zeroth-Order Polynomial (VZOP) fitting and use it to fit the simulation data. This algorithm is developed for the upcoming Low-frequency Anechoic Chamber Experiment (LACE), yet it is a general method suitable for application in any single antenna-based global 21 cm signal experiment. VZOP defines a 24-hour averaged beam model that brings information about the antenna beam into the polynomial model. Assuming that the beam can be measured, VZOP can successfully recover the 21 cm absorption feature, even if the beam is extremely frequency-dependent. In real observations, due to various systematics, the corrected measured beam contains residual errors that are not completely random. Assuming the errors are frequency-dependent, VZOP is capable of recovering the 21 cm absorption feature even when the error reaches 10%. Even in the most extreme scenario where the errors are completely random, VZOP can at least give a fitting result that is not worse than the common polynomial fitting. In conclusion, the fitting effect of VZOP depends on the structure of the error and the accuracy of the beam measurement.


INTRODUCTION
During the final phase of the Dark Age, the first generation of stars emerged under the gravitation effect (Castelvecchi 2019).Afterwards, the first stars emit UV and X-ray photons to heat and ionize the neutral hydrogen in the universe (Dayal & Ferrara 2018).Therefore, this period is called the Cosmic Dawn (CD) and Epoch of Reionization (EoR).The cosmic 21 cm signal which comes from the hyperfine splitting of the 1S ground state of hydrogen is an important probe to study the Dark Age, CD and EoR.Because the baryonic matter in the universe is mainly hydrogen, it is possible to detect 21 cm signal -a forbidden line in the laboratory -by differential brightness temperature between spin temperature of hydrogen and background brightness temperature (Morales & Wyithe 2010;Pritchard & Loeb 2012;Zaroubi 2012;Mesinger 2019).We always take Cosmic Microwave Background (CMB) as background to detect cosmic 21 cm signal.When  ≲ 200, the difference between spin temperature and CMB temperature kept changing due to several factors such as the cosmic expansion, different ionization and recombination mechanisms, heating and cooling processes, and the evolving of the number density of different kinds of ions, showing different absorption or emission features in different periods.Therefore, the details of the pictures of ★ E-mail: jhgu@nao.cas.cn(JG) † E-mail: guoquan@shao.ac.cn (QG) the universe can be inspected through the detection of the cosmic 21 cm signal.
However, there are several difficulties in detecting 21 cm signal.A widely accepted prediction about the biggest differential from CMB is only about 200 mK (Cohen et al. 2018;Reis et al. 2021), while Galactic foreground and extragalactic radio sources are stronger by 4-5 orders of magnitude.In addition to the foreground, there are many other contaminants such as radio frequency interference (RFI), ground radiation temperature, and thermal noise of the receiver.All of these contaminants are always much greater than the signal.
There are mainly three strategies to detect cosmic 21 cm signal: 1) global 21 cm signal detection, which is used by experiments including BIGHORNS (Broadband Instrument for Global HydrOgen ReioNisation Signal, Sokolowski et al. 2015a,b), EDGES (Experiment to Detect the Global EOR Signature, Bowman et al. 2008;Rogers & Bowman 2012;Mozdzen et al. 2016b;Bowman et al. 2018), LEDA (Large-Aperture Experiment to detect the Dark Ages, Greenhili et al. 2012), SCI-HI (Sonda Cosmológica de las Islas para la Detección de Hidrógeno Neutro, Voytek et al. 2014), PRI  M (Probing Radio Intensity at high-Z from Marion, Philip et al. 2019), SARAS (Shaped Antenna measurement of the background RAdio Spectrum, Singh et al. 2017;Subrahmanyan et al. 2021), CTP (Cosmic Twilight Polarimeter, Nhan et al. 2019), REACH (Radio Experiment for the Analysis of Cosmic Hydrogen de Lera Acedo 2019), 2) 21 cm signal power spectrum detection, which is performed by some SKA (Square Kilometer Array) pathfinders such as 21CMA (21 Centimeter Array, Zheng et al. 2016;Huang et al. 2016), GMRT (Giant Metrewave Radio Telescope, Paciga et al. 2013), MWA (Murchison Widefield Array, Beardsley et al. 2016), LOFAR (Low-Frequency Array, Patil et al. 2017), PAPER (Precision Array for Probe the Epoch of Reionization, Ali et al. 2015) and HERA (Hydrogen Epoch of Reionization Array, DeBoer et al. 2017), 3) 21 cm signal tomography, which is one of the main tasks of SKA in the future.However, no valuable results have been obtained by now except EDGES, which claimed an absorption profile centered at 78 MHz with half-maximum width 19 MHz and 0.5 K depth (Bowman et al. 2018).If this profile is indeed the cosmic 21 cm absorption feature, it will unveil new physics or astrophysics (Chen & Miralda-Escudé 2004;Nelson et al. 2013;Barkana 2018;Li et al. 2021;Yang 2021).Nevertheless, this result is not yet widely accepted because of serious divergence from theory (Cohen et al. 2018).SARAS studied this profile and found that it may be spectral distortion caused by EDGES low-band instrument (Singh et al. 2022).
Global 21 cm signal detection usually uses just a single antenna to obtain the global averaged spectrum of the sky brightness temperature but ignores its spatial variation.Benefiting from the portability of single antennas, it could be relatively easy to select the experiment site.Furthermore, only a few data acquisition channels are required, and no complex data exchange network is needed.Among various Global 21 cm signal detection experiments, only EDGES claimed that they successfully detected the signal.EDGES fitted a fifth-order polynomial to data in log space and sampled in the parameters space to detect the absorption feature (Bowman et al. 2018).Nth-order polynomial fitting is a common algorithm in many global cosmic 21 cm signal detection simulations (e.g.Pritchard & Loeb 2010;Voytek et al. 2014;Bernardi et al. 2016;Mozdzen et al. 2016b;Bowman et al. 2018;Gu & Wang 2020;Singh et al. 2022;Shi et al. 2022).This method is based on the assumption that the Galactic foreground spectrum can be well described by a power-law model as it is dominated by the synchrotron radiation because of free electrons spiraling around the Galactic magnetic field lines.
However, the frequency-dependent (chromatic) beam pattern introduces an artificial spectral structure, which significantly deviates from a simple power law model and hinders the reconstruction of the 21 cm signal.Researchers typically divide the global spectrum by a beam correction factor to account for the chromaticity of the antenna.The beam correction factor is defined by where   is the real sky temperature at frequency  and direction n,   is the antenna beam,   is the reference frequency and Ω is the solid angle (e.g.Mozdzen et al. 2016a;Monsalve et al. 2017;Mozdzen et al. 2019;Sims & Pober 2020;Shen et al. 2021;Anstey et al. 2022;Spinelli et al. 2022;de Lera Acedo et al. 2022).However, due to the varying spectral index of foreground temperatures in different directions, the effectiveness of the beam correction factor is often limited.Therefore, some researchers have parametrized the foregrounds (Anstey et al. 2021;Pagano et al. 2022).They divided the sky maps into N regions with similar spectral indices and assumed that the spectral index within each region is constant.Assuming the antenna beam can be precisely measured, they employed a Bayesian approach to estimate the spectral indices of each region, allowing for a more accurate foreground subtraction.This result can be anticipated since it incorporates the information of the antenna beam into the signal reconstruction model.Similarly, the algorithm described in this paper leverages the same principle and further takes into account the scenario of imprecise beam measurements.We develop the Nth-order polynomial fitting method above, i.e.Vari-Zeroth-Order Polynomial (VZOP hereafter) fitting, which will be used in a ground-based global 21 cm signal experiment named "Low-frequency Anechoic Chamber Experiment" (LACE Huang et al. 2021) in the future.LACE plans to build a lowfrequency underground anechoic chamber to provide a radio-quiet environment, the schematic of the chamber is shown in the Fig. 1 of Huang et al. (2021).We intend to apply the VZOP in LACE in the future.Consequently, some of the parameters chosen for our simulation are aligned with the context of LACE.Nevertheless, it is crucial to note that VZOP is not exclusively tailored for LACE; rather, it is a general algorithm and can be employed effectively in any single antenna-based global 21 cm signal experiments.The rest of this paper is organized as follows.Section 2 presents the foreground, cosmic 21 cm absorption feature and noise model used in the simulation, and details the VZOP algorithm.In Section 3 we describe the simulation process.We provide the fitting results assuming that the beam is accurately known in Section 4 and that the beam has errors in Section 5. We discuss and conclude in section 6.

Foreground model
In this work, we make use of the GSM (Galactic Global Sky Model, de Oliveira-Costa et al. 2008), which provides a sky maps generator that can generate maps at any frequency between 10 MHz and 10 GHz.GSM used the first 3 principal components out of 11 components obtained using PCA in the sky region with observed data at all 11 frequencies, then fitted 3 principal components with a cubic spline at any frequency to get all-sky maps.It should be noted that the foreground emission is slightly polarized, which is direction-and frequency-dependent (Spinelli et al. 2019).This could also create an artificial spectral structure in the measured global spectra.In this work, we do not consider this factor to keep the conciseness.Further investigation into characterising this factor will be the aim of future work.

Cosmic 21 cm signal models
This work focuses on the absorption profile in 50 − 120 MHz, which has the largest amplitude and is the most significant signature of the global EoR 21 cm signal.In this paper, we choose two signal models as input to simulate the expected 21 cm absorption feature: the Gaussian model refers to a simple Gaussian profile (Huang et al. 2021), and the EDGES model refers to a flattened Gaussian profile found by EDGES (Bowman et al. 2018).The 21 cm signal brightness temperature  eor () refers to the absorption or emission feature on CMB at frequency .
Gaussian model: This model is an approximation of the theoretical model, which can be written as where amplitude  = −0.150K, center frequency   = 78.3MHz, width  = 5 MHz.EDGES model: EDGES have found a flattened Gaussian profile centered at 78.3 MHz, which can be expressed as Teor.

Thermal noise
We adopt a simple thermal noise model (Shi et al. 2022) with the following standard deviation where   is the antenna temperature,  m = 1 is the number of independent measurements, and we have chosen channel bandwidth Δ = 1 MHz and integration time  int = 10 d, as such configurations ensure that the RMS remains below 10 mK, thereby exerting no substantial influence on the signal extraction.This noise was first obtained by Dicke (1946).For any specific system, the enhancement in the RMS noise cannot be better than as given in equation ( 4) (Wilson et al. 2009).

Algorithm
The first challenge to detecting the cosmic 21 cm signal is the interference from the Galactic foreground radiation.Within 50-120 MHz, the Galactic foreground temperature is as high as ∼ 10 3 K, and the maximum amplitude of the absorption dip of the 21 cm signal is just ∼ 10 −1 K. Subtracting contaminants that are 4 to 5 orders of magnitude higher than the signal is extremely difficult.The antenna temperature at a certain frequency is the weighted average of the summary of the foreground and 21 cm signals at different positions in the sky.Therefore, the antenna temperature   () at  contains a frequency-dependent antenna beam, which can be expressed as ( where T () is the arithmetic mean of the observed sky temperature (i.e.assuming a uniform isotropic antenna).
The common algorithm is to do a polynomial fit in logarithmic space.
where the hat indicates a modelled quantity, zeroth-order term T (  ) is the modelled value of the antenna temperature at reference frequency   , and   is the coefficient of the th-order term.If the antenna temperature spectrum does contain the EoR component, then the residuals from the antenna temperatures minus the fitted values will show some structures.Next, assuming a 21 cm absorption feature model, adding a polynomial model and refitting, we have If the residuals obtained by subtracting the fitted values from the antenna temperatures become small and uniform, the fitted signal parameters are considered to be true parameters.However, this algorithm usually does not perform well.The chromaticity of the antenna introduces additional structure in the temperature spectrum, which makes the reconstructed signal deviate from the expected signal (e.g.Wilensky et al. 2022;Mozdzen et al. 2016b).This requires researchers to carefully design the antenna and handle the effects of chromaticity (e.g.Shi et al. 2022;de Oliveira-Costa et al. 2008;Subrahmanyan et al. 2021;Philip et al. 2019;Scheutwinkel et al. 2022;Shen et al. 2022).The problem introduced by the chromatic antenna beam could be partly solved by containing the frequency dependence of the beam pattern in the solution process.
In this paper, we use VZOP to improve signal reconstruction by fitting chromaticity-generated structures to an improved polynomial model.Averaging the beam and sky temperature over 24 hours to define a 24-hour averaged beam model B(, ) and a 24-hour averaged temperature model T (, ), where  is the declination.Then, we can obtain the global averaged temperature as follows: where  is the right ascension, the directional antenna gains below the horizon are set to be 0. T (, ) and B(, ) are just dependent on frequency  and  because the integral relating to  is eliminated after 24 hours averaging.The Galactic foreground is regarded to be roughly a power law spectrum, and the antenna temperature without the signal can be approximated as a separable variable, i.e.  (, n) =  (/  ) ×  (n,   ).Therefore, the modeled value of T (, ) can be written as a polynomial in log space where Tgal (, ) is the modelled value of mean Galactic brightness temperature averaging over 24 hours, Teor () is the modelled value of 21 cm signal mean brightness temperature averaging over 24 hours (21 cm signal is independent of ), and we set   = 85 MHz.We define ) . (10) According to equation ( 8), ( 9) and ( 10), discretizing angle and frequency Writing equation ( 11) as matrix form where The  gal is a vector containing real Galactic temperatures at different declinations at frequency   .
In this paper, we take  bin bins at equal intervals between the declination interval [ ℎ ,   ] (If the observation site is in the northern hemisphere,  ℎ = 90 • ,   is the southernmost declination that can be seen; if it is in the southern hemisphere,  ℎ is the northernmost declination that can be seen,   = −90 • .).In this paper, we use a 5th-order VZOP to fit our mock data, and unless otherwise specified,  bin = 10.
We use a truncated likelihood as our posterior distribution (The upper and lower bounds for each parameter are shown in Table 1) where   = 71 is the number of frequency channels, and  is a diagonal matrix because we assume that the thermal noise of different frequency channels is independent, and the -th element on the diagonal is Even if the assumed frequency independence does not hold, VZOP can still be used because it does not rely on this assumption.In addition, the first term on the right side of equation ( 13) is the normalization coefficient.Although this term can be ignored during Bayesian parameter estimation, it is necessary to include it when computing the Bayesian evidence afterwards.Theoretically, VZOP would outperform polynomial fitting algorithms.The principle of VZOP is that since the chromaticity of the antenna cannot be avoided, the distribution of the foreground temperature is simulated by sampling 24-hour averaged temperatures in several declination bins.From a mathematical point of view, VZOP is equivalent to an improved polynomial fitting model with the zerothorder term that varies with frequency, and the simplified expression is (according to equation ( 9)) The degree of freedom of the zeroth-order terms is equal to the rank For the convenience of discussion, all matrices B constructed in this paper are full-rank matrices, and the number of columns is not greater than the number of rows (i.e. the number of declination bins is not greater than frequency channels).Formally, VZOP can also obtain the 24-hour averaged temperature within each declination bin, but we are not interested in this nuisance parameter.We use a python package emcee (Foreman-Mackey et al. 2013) and find the posterior is a multi-modal distribution.To find the global maximum quickly and efficiently, we use the python package ptemcee, which is an improved version of emcee using parallel tempering (Vousden et al. 2016).The package ptemcee simultaneously samples from tempered versions of the posterior distribution (), where  () and () are respectively likelihood and prior distributions.The probability distribution at the lowest temperature  = 1 corresponds to the target posterior distribution, and as the temperature increases, the distribution gradually approaches the prior distribution.At the highest temperature, ptemcee approximates sampling from the prior distribution.In this way, the Markov chain at low temperatures is responsible for exploring a specific mode, while the Markov chain at high temperatures can explore the entire posterior space.By exchanging information between different temperatures, ptemcee can effectively sample the high-dimensional and multimodal distribution.

SIMULATION
In our simulation, we positioned the antenna at the site of the 21CMA (42.93 • N, 86.68 • E), which is far away from artificial radio sources.
The antenna to be used in LACE has not been finalized at this time.An important principle in antenna design is to minimize the antenna chromaticity and to avoid the coupling of spatial structures into the frequency domain as much as possible.Three candidate antenna designs are to be tested: a nominal blade antenna, an optimized blade antenna and a bow-tie blade antenna.Each of them consists of two identical Perfect Electric Conductor (PEC) panels with a thickness to be 2 mm.They work in frequency bands 50 − 120 MHz corresponding to redshifts at  = 10 − 27 where the strong absorption feature is expected to exist.The ground plane is assumed to be an infinite uniform PEC.The structures of the three simulated antennas are shown in Fig. 2. The cross-section of antenna beam profiles, cut at zenith angle Θ = 30 • with different azimuth angle Φ, are shown in Fig. 3.
Under the assumption that the antenna beams do not change with time, we use the following method to generate mock observation data i.e. the antenna temperature.
1) Given an arbitrary frequency, the GSM model can generate a global map using healpix RING scheme with an angular resolution of approximately 6.87 ′ ( side = 512) (Gorski et al. 2005).2) We divide the sky map into 1024 bins along the declination of the celestial coordinate system and take the averaged temperature in each bin as the 24-hour averaged brightness temperature at the declination of the centre of the bin.3) We could perform the same for the antenna beam, but a preprocessing is required because the simulated beam is not equal in area but equal in angle.The zenith and azimuth angles of our generated antenna gains are both 2 • apart (the gains below the horizon are set to 0), so the resolution is much lower than 6.87 ′ .Firstly, therefore, we calculate the smallest  side = 16 at which there is at least one gain within each pixel.We obtain the arithmetic mean if there is more than one gain within each pixel.Secondly, we perform bilinear interpolation based on the nearest four points and expand it to a map with  side = 512.4) The antenna is placed at the observation site at a certain moment (the time can be chosen arbitrarily because it does not affect the final result).Next, similar to step 2), the gains within each bin are summed and averaged to obtain the 24-hour averaged beam.5) Different from the common polynomial fitting, if we want to generate the antenna temperature by multiplying the 24-hour averaged temperatures and the 24-hour averaged beam, we need to multiply the beam by the cosine of the corresponding declination and then normalize, and then the result of the multiplication is the antenna temperature as following, In Fig. 4, we show the sky brightness temperature weighted by the 24-hour averaged antenna beam of the optimized blade antenna at 85 MHz.The temperature distribution of the figure is independent of the longitude of the observation site.The figure shows the map using our defined 24-hour averaged antenna beam after taking into account the effects of Earth's shading.It tells us which part of the sky has made the main contribution to observational data.The darker area represents the area of the sky near the southern horizon that is only visible for part of the day.The grey area at the bottom represents the area below the horizon throughout the day.

FITTING RESULTS
In this section, we performed signal reconstruction at the 21CMA site assuming that the antenna beam can be precisely measured.

Common Polynomial Fitting
Before using the VZOP fitting, we first present the result of common polynomial fitting.Fig 5 shows the results obtained by fitting the antenna temperatures generated by different antennas using the polynomial fitting (Gaussian model), where a uniform isotropic antenna is added for comparison.When using the uniform isotropic antenna, it is not surprising that the reconstructed absorption feature agrees well with the input feature.When using other antennas, the results of signal reconstruction are relatively poor due to the existence of chromaticity.
Since the uniform isotropic antenna is achromatic, it is predictable that the reconstructed absorption feature and the input feature are perfectly coincident In Fig 5 .However, the fitting results based on the mock data obtained from the other three antennas have different degrees of deviation.The optimized blade antenna has the least deviation, while the nominal blade antenna has the strongest deviation from the input absorption feature.This is consistent with Fig 3, where the frequency dependence of the optimized blade antenna beam is the weakest, and that of the nominal blade antenna is the strongest.
It can be seen from the above results that the use of non-ideal antennas will lead to poorer performance of polynomial fitting to extract the 21 cm absorption feature.The stronger the frequency dependence of the antenna, the worse the performance of the polynomial fitting.Because of the different foreground temperatures at different points in the sky, the spatial variation of the beam with frequency will couple into the spectrum, which will eventually lead to additional structure in the spectrum and affect the final fitting result.

Mitigating the Chromatic Effect with VZOP
The chromaticity will introduce additional structure on the 24-hour averaged antenna temperature spectrum and have an impact on signal reconstruction.Therefore, we first fit the foreground model with only noise (i.e.no 21 cm signal) using VZOP.The fitting results are shown in Fig 6 .For comparison, we also introduce a uniform isotropic antenna.The lines are shifted vertically for clarity.Since we generate the noise using a random number generator with the same seed for all antennas, the shapes of these lines are very similar.It can be seen that the residuals after fitting the foreground are similar to that of input noise.The Root Mean Squares (RMSs) of the residuals of these lines are less than 2 mK, far less than ∼ 10 2 mK of the 21 cm signal.
The RMS of a uniform antenna is small and its residuals are uniform, which is reasonable because the uniform antenna has no chromaticity.The RMSs of the other three antennas are similar to that of the uniform antenna, indicating that VZOP effectively fits the structure introduced by chromaticity into the improved polynomial model.

VZOP Fitting
As described in Section 4.2, VZOP effectively mitigates the effects of chromaticity.Therefore, in the reconstruction of the cosmic 21 cm absorption feature, the performance of VZOP will be better than the common polynomial fit.We plot in Fig. 7 the fitting results for Gaussian model.We can see in these three panels that if the EoR signal is included, the RMS will decrease from over 20 mK to only about 2 mK.The reconstructed EoR absorption features almost completely coincide with the input features for both antennas, indicating that VZOP fitting can effectively reduce the impact of the antenna chromaticity for Gaussian model.
Next, we plot the fitting results for EDGES model in Fig. 8.We obtained similar results to Fig. 7, with the reconstructed absorption feature being highly coincident with the input feature.Compared with Gaussian model, the fitting result of EDGES model is slightly

= 500+]
= 850+] = 1200+] take the set of parameters with the largest posterior probability as the best fitting value, whether it is the Gaussian model or the EDGES model, the results obtained by each run of MCMC have some fluctuations, and the fluctuation of the EDGES model is greater.This may be because the EDGES model is deeper, wider, and has more parameters so that is more difficult to constrain.But all of their best-fit values are within the 68% confidence interval, so such fluctuations have little effect on the final results.The parameters distribution of the middle panel in Fig. 7 (i.e.Gaussian model based on the optimized blade antenna) is shown in Fig. 9.As can be seen from this figure, both the polynomial coefficients and the 21 cm signal parameters follow normal distributions approximately.But our model has just weak constraints on the 10 24-hour averaged temperature parameters.We found that the true 24hour averaged temperatures lie within 1 of the best-fit value, which means we obtain the true 24-hour averaged temperatures within large error ranges.But it does not make much sense because of too large errors, we treat them as nuisance parameters and ignore them.
We list the fitted 21 cm signal parameters for Gaussian model and EDGES model in Table 2 and Table 3, respectively.Except for  in EDGES model, all input parameters are within the 1 error interval of the recovered parameters.To verify the plausibility of the model themselves, we also calculated the Bayesian evidence using polychord (Handley et al. 2015a,b), which are also placed in two tables.It can be seen from these two tables that if the 21 cm absorption feature is considered during sampling, the evidence significantly increases while the RMSs significantly decrease.This indicates that the signal reconstruction results of VZOP are highly reliable.

Chromatic Beam Effect
It can be seen from the above results that the VZOP fitting algorithm has an exact signal reconstruction for all three antennas and can largely reduce the impacts of chromaticity.However, the three antennas used in this paper are all close to achromatic antennas.To better understand the effect of antenna frequency dependence on VZOP, we introduce four other antennas with higher frequency-dependency, the beam pattern of which are shown in

Number of Declination Bins
It can be seen from Fig. 7 and Fig. 8 that our reconstructed EoR absorption feature and input feature are highly coincident for all antennas when we sample 10 declination bins, whether it is Gaussian model or EDGES model.These results show that VZOP fitting reduces the effects of antenna chromaticity very well.
To determine at least how many bins obtain the best results, we refit with different  bin for two models.We use an optimized blade antenna as the fiducial antenna because of the weakest frequency correlation.Fig 12 shows the fitting results, where the common polynomial algorithm is added for comparison.
When we use only 1 bin, VZOP fitting will degrade into the polynomial fitting in theory.As can be seen from the two panels in Fig 12, as expected, the result obtained using 1 bin is highly coincident with that obtained by the common polynomial fit algorithm.However, with only 2 bins, the reconstructed absorption feature is almost identical to the input feature for both cases.There is no significant improvement in fitting even if continually increasing  bin .The results in Fig. 12 show that for the optimized blade antenna, only 2 bins are sufficient to recover the 21 cm absorption feature.For a well-designed antenna, VZOP will carry enough information of the beam with only a few declination bins.

INACCURATE ANTENNA BEAM PATTERN
We have seen that the VZOP fits very well for a wide range of commonly used antennas.The above results confirm our assumption, that as long as the information of the antenna beam is brought into the fitting by some method, the foreground can be deducted efficiently.Therefore, VZOP fits very well for a wide range of commonly used antennas when the beam can be accurately known, and usually can extract 21 cm absorption feature by using just 2 declination bins.Unfortunately, in actual observations, the antenna beam pattern can hardly be precisely measured due to various systematics such as reflections from surrounding objects, ionospheric effects, noise inside the instrument, etc.Even performing a systematic correction, it is expected to still have 1% -10% residual systematics.Therefore, in this section, we will study the performance of VZOP when there are errors in the beam.Specifically, we will use an accurate beam in simulation and add a Gaussian error to each pixel of the beam model used in VZOP.Note that the random error will be smoothed out in the process of arithmetic averaging, the error of the 24-hour averaged beam will be much less than that of the origin antenna beam.
The errors at different pixels and frequencies are not completely random, as they arise from systematics.In other words, the errors at different pixels and frequencies exhibit correlation.Therefore, we assume that the errors at different pixels are random, but the errors at different frequencies within the same pixel are correlated.We studied four models: constant, linearity, quadratic and cosine.The details of them are shown in Table 4.We consider these four error models and simulate observations using three fiducial antennas respectively (using the Gaussian model).For all the situations, using the VZOP considering 10 declination bins can successfully extract the 21 cm absorption feature, even if the error reaches 10% which is a value the current correction technology can completely make the remaining error of the beam less than it.
Compared with polynomial fitting, VZOP has a higher degree of freedom, so VZOP performs well even with a relative error of up to 10%.Only considering the cosine error model with 10% error, we show the effect of the number of declination bins on the fitting results under the three fiducial antennas in Fig. 13 (Gaussian model) and Fig. 14 (EDGES model).From the six panels in the two figures, on the whole, when only two declination bins are used, the fitting results are close to that of the polynomial fitting.As the number of bins increases, the fitting results gradually become better, and the output absorption feature mostly coincides with the input feature when 5 bins are used.But there are some exceptions.Such as the middle panel of the Fig. 13 shows an abnormal situation that the output feature of  bin = 3 is more accurate than  bin = 4.This may just be due to the error in our model.When  bin = 3, the error is offset by a large part.When the number of bins continues to increase, the error increases again.On the other hand, comparing different columns of the two figures (i.e.different antennas), using 10 bins gives an exact fit for any antenna.When the number of bins is too small to extract accurately the input absorption feature, the fitting results of the optimized blade antenna that is weaker frequency dependence are better.
Based on the aforementioned results, VZOP works well even in the presence of 10% errors.Furthermore, we evaluate the performance of VZOP under the worst-case scenario, where the errors at different pixels and different frequencies are completely random.We added different Gaussian errors to the optimized blade antenna and reextracted the 21 cm absorption feature, and the signal reconstruction results are shown in Fig. 15.The output absorption features are highly consistent with the input features for both two signal models when there are no errors.As the error increases, the reconstruction results gradually become worse and closely approximate to that of polynomial fitting when the error magnitude reaches 0.6%.In other words, 0.6% error is enough to cause VZOP to lose its advantage.However, it is worth noting that as the error continues to increase, the fitting performance of VZOP does not deteriorate further.After testing, no matter how large the error is, VZOP will not be worse than the polynomial fitting, because VZOP has more degrees of freedom.We have even set the gain on each pixel to an arbitrary value, and the results did not become worse than that of polynomial fitting.Hence, on the one hand, the actual errors of the antenna beam are certainly not entirely random.On the other hand, even in the worst-case scenario, the fitting results of VZOP do not deteriorate compared to that of common polynomial fitting.

DISCUSSIONS AND CONCLUSIONS
In this work, we reconstruct the cosmic 21 cm absorption feature from foreground contamination using VZOP -an improved polynomial fitting algorithm -by defining a 24-hour averaged beam model.We plan to apply VZOP in the forthcoming LACE.Nevertheless, in its capacity as a general algorithm, VZOP can be employed in any single antenna-based global 21 cm signal experiment.VZOP is introduced to address a specific problem that the chromaticity of the antenna couples the spatial structure of the Galactic bright temperature map into the frequency structure, making the common polynomial fitting work poorly and hindering the reconstruction of the cosmic 21 cm signal.If we know the antenna beam, we can invert the true global spectrum, which is the idea utilized by VZOP that brings the information of the beam into the reconstruction model such as polynomial fitting.From a mathematical viewpoint, VZOP is equivalent to an improved "polynomial fitting" model, in which the zeroth-order item is a quantity that varies with frequency and beam-related information is included in the zeroth-order item.The spectral structure introduced by the antenna chromaticity can be deducted by using the beam information to fine-tune the polynomial "constant term".
Assuming the antenna beam can be accurately measured, VZOP will accurately deduct the extra structure introduced by chromaticity (even with only 2 declination bins), and the reconstructed 21 cm absorption feature is highly coincident with the input feature, which is consistent with expectations.The fact that VZOP works with only 2 bins suggests that additional structures may be described with only one parameter.In addition, VZOP performs well even for extremely frequency-dependent antennas.
In actual observation, we cannot accurately measure the antenna beam due to various system errors.It is not realistic to assume an accurate antenna beam pattern in the solving process.Therefore, taking into account that the errors are from systematics, such as the antenna dimensions errors, we assume that they are random at different pixels and follow several models outlined in  presence of errors, more declination bins may be required to ensure that VZOP performs well.For example, under the cosine error model, at least 5 bins are required to recover the Gaussian model and the EDGES model.
Finally, we investigate the fitting results of VZOP under the worstcase scenario, where the errors are completely random across different pixels and frequencies.As long as the error exceeds 0.6%, it is enough to interfere with VZOP's identification of additional structures and make VZOP lose its advantage.However, an interesting phenomenon is that no matter how large the error is, the fitting results of the VZOP are not worse than that of the polynomial fit.The reason is that the essence of VZOP is to fine-tune the polynomial "constant term" to identify the extra structure introduced by chromaticity, increasing the error only reduces the ability to identify the extra structure.Based on empirical observations, it can be seen that as the error increases, the variation of the zeroth-order term to frequency becomes progressively smaller, eventually converging to a "constant".This feature ensures that even in the worst-case scenario, the performance of VZOP will not be worse than the common polynomial fitting, let alone the fact that the errors are certainly not completely random.
In summary, the fitting performance of VZOP depends on the accuracy of beam measurements and the error structure.More bins are required to recover the 21 cm absorption feature when errors exist due to the error reducing the amount of information carried by each bin, which is why VZOP's performance deteriorates in the case of completely random errors.When the error is completely random, the degree of freedom of the error becomes too high, resulting in little information carried by each bin.It is worth noting that using too many bins can lead to overfitting.To address this issue, using a narrower frequency bandwidth could be considered.Because then there will be more frequency channels, which will enable the use of more declination bins and improve the performance of VZOP in the presence of beam errors.This could be an area of future research.
In the simulation, VZOP can successfully extract the cosmic 21 cm absorption feature, but the real observation may be more complicated.Whether it is the Gaussian model or the EDGES model, it is only the approximate description of the real cosmic 21 cm signal.On the other hand, the true residual error may not be described by a few simple models.In the future, we will continue to study the performance of VZOP in more realistic observations.

Figure 1 .,
Figure 1.Brightness temperature profiles of two cosmic 21 cm absorption feature models.The two lines represent the Gaussian model (blue solid) and the EDGES model (orange dashed), respectively.

Figure 2 .
Figure 2. The simulated images of nominal blade antenna (left), optimized blade antenna (middle) and bow-tie blade antenna (right).

Figure 6 .
Figure6.The residuals after fitting foreground with only noise for different antennas.The dashed lines represent the input noise and the solid lines are residuals.From the bottom up, the lines with different colours represent residuals for uniform isotropy (orange), nominal blade (vermilion), optimized blade (blue) and bow-tie blade antenna (reddish purple), respectively.The lines are shifted vertically for clarity.
Fig 10.For the convenience of comparison, panels (a), (b) and (c) are all optimized blade antennas with different  ( = 1.25, 1.5, 2 m, respectively), and panel (d) is the bow-tie blade antenna with  = 3 m.Compared with Fig 3, the four antenna beams in Fig 10 vary significantly with frequency.When the  = 120 MHz, the beams in panels (b), (c), and (d) have obvious changes.The antenna in panel (b) has very obvious side lobes; the main lobe of the antenna in panel (c) is rotated 90 • ; the beam in panel (d) also has a great change.Fig 11 shows the fitting results for Gaussian model based on the four extreme frequency-dependent antennas.panels (a) and (b) still show near-perfect signal reconstruction results.The reconstructed EoR absorption feature in panel (c) deviates obviously from the input feature and the ranges of the parameter distribution of both panel (c) and (d) are significantly wider.Fig 11 shows that, assuming the beam is accurately known, VZOP still performs well even with extreme frequency dependence.

Figure 9 .
Figure 9. Distributions of all parameters for Gaussian model based on the optimized blade antenna.The parameters from left to right (from top to bottom) are 5 polynomial coefficients ( 1 −  5 ), 10 24-hour averaged temperature parameters at each bin and signal parameters (,   , ).The true values of cosmic 21 cm signal parameters are also shown in a black solid line.

Figure 11 .
Figure 11.Fitting results for Gaussian model with different antennas (optimized blade antennas with  = 1.25, 1.5, 2 m and bow-tie blade antenna with  = 3 m).The reddish-purple dash-dotted line represents the residuals after fitting and removing both foreground and absorption features.The yellow solid line denotes the input absorption feature, while the black solid is the recovered best-fit line.Also, we plot dozens of lines to represent different error levels with different shades of colour (see the colour bar to the right of the figure).For clarity, we intentionally plot two 1 lines in (dashed black lines).

Figure 12 .
Figure 12.Fitting results using different numbers of declination bins for Gaussian model (top panel) and EDGES model (bottom panel) with optimized blade antenna.

Figure 13 .
Figure 13.Fitting results for Gaussian model using different number of declination bins with nominal blade (left), optimized blade (middle) and bow-tie blade antenna (right).For comparison, we also plot the results of the polynomial fit.

Figure 14 .
Figure 14.Fitting results for EDGES model using different number of declination bins with nominal blade (left), optimized blade (middle) and bow-tie blade antenna (right).For comparison, we also plot the results of the polynomial fit.

Figure 15 .
Figure 15.Fitting results when adding different Gaussian errors to the optimized blade antenna for Gaussian model (top panel) and EDGES model (bottom panel) assuming the errors are completely random.We also plot the results from polynomial fitting for comparison.

Table 1 .
The upper and lower bounds for each parameter in the uniform prior distribution.When the number of frequency channels is fixed (i.e. the number of rows of B is fixed), if the rank of B increases with declination bins (i.e. the number of columns of B increases), the fitting result will be better theoretically.The maximum number of declination bins should not exceed the number of frequency channels.
Table. 4 at different frequencies.VZOP can still successfully extract the 21 cm absorption feature even if the error reaches 10%.However, in the