Calibration Error in 21-cm Experiments Global Spectrum Experiment

1. Introduction

The redshifted 21 cm line of neutral hydrogen, potentially observable at low radio frequencies (50-200 MHz), should be a powerful probe of the physical conditions of the inter-galactic medium during Cosmic Dawn (CD) and the Epoch of Reionization (EoR) [1,2,3,4,5].

The global (all-sky averaged) spectrum of the redshifted 21 cm brightness temperature can be measured with a single antenna and wide band spectrometer. It is a direct approach with high precision instrument. Experimental efforts include the Experiment to Detect the Global EoR Signature (EDGES, [6,7]), the Broadband Instrument for Global HydrOgen ReioNisation Signal (BIGHORNS, [8]), the Shaped Antenna measurement of the background RAdio Spectrum (SARAS, [9,10,11]), the Probing Radio Intensity at high-Z from Marion (${\mathrm{PRI}}^{\mathrm{Z}}\mathrm{M}$,12]), Radio Experiment for the Analysis of Cosmic Hydrogen (REACH, [13,14]), the Large-aperture Experiment to Detect the Dark Age (LEDA, [15]), Sonda Cosmologica de las Islas para la Deteccion de Hidrogeno Neutro (SCI-HI, [16]), Cosmic Twilight Polarimeter (CTP, [17]) and Mapper of the IGM Spin Temperature (MIST, [18]), which use interferometric measurements to help estimate instrumental and foreground parameters. The EDGES experiment has reported the detection of a 500 mK deep absorption feature centered at 78 MHz [19], which may be associated with the signature of the cosmic dawn. However, the measurement taken by the SARAS-3 experiment has not detect such an absorption feature [20]. Considering the measurement difficulty from most ground sites at low frequencies due to ionosphere refraction and reflection of broadband radio frequency interference (RFI), there are some proposals to go to the far side of the moon, for example the Dark Ages Polarimetry PathfindER in low lunar orbit (DAPPER, [21]), Dark Ages Radio Explorer project (DARE, [22]), and the Discovering the Sky at Longest wavelength project (DSL, [23]).

Observing the global 21 cm signal is very challenging, as the low-frequency radio sky is dominated by intense synchrotron emission from our own Galaxy which is more than four orders of magnitude brighter than the signal. There are also Galactic free-free emission, numerous compact radio sources such as supernova remnants, quasars and radio galaxies which contribute to the total received power. Extremely high sensitivity and large dynamic ranges are required to discern the small 21 cm signature in the spectrum. Moreover, the cosmic dawn 21 cm signal has an unknown but generally broad shape, its detection requires a good understanding of the instrument response, which must be determined by a calibration procedure. Any frequency-dependence in instrumental gain, noise spectrum, antenna beam shape, and ionospheric effect, if not properly accounted for, can affect the measurement result. In particular, the reflections of electric signal at the interfaces of the various components of the measurement system are frequency-dependent and could not be neglected. Simple estimates of system temperature and noise factor would be inadequate for the proper characterization of such a precision experiment.

The noise wave concept [24] provides a framework for describing such circuit system at high precision [25]. In this formulation, the linear radio frequency (RF) devices are characterized by a few fundamental wave parameters, which can be determined experimentally by measuring the reflection coefficients. It is then possible to relate the noise power measured by the instrument to the sky temperature, and give an estimate of the error in measured results.

To assess the uncertainties of the measurement, besides the thermal noise, one also needs to consider possible systematic errors which may bias the result. This is the main aim of the present work. Previously, Monsalve et al. [26] and Monsalve et al. [27] presented the calibration strategy of the EDGES experiment, and analyzed the propagation of uncertainties in the receiver calibration. In a recent work [28], a Bayesian approach to the noise wave formulation problem (the one used by the REACH experiment [13]) is presented. In this work, we investigate this problem with realistic models, and a detailed characterization of the systematic errors. In particular, we consider the Vector Network Analyzer (VNA) measurement error which depends on frequency, and adopt a comprehensive model for VNA measurement uncertainties, include all measurement errors in VNA measurement procedure and in calibration standards, to estimate theirs impacts on the spectrum measurement system. Furthermore, we also study how the error affect the extracted 21cm signal, i.e. the non-smooth component of the sky spectrum.

The rest of this paper is organized as follows: Section 2 describes the foreground model and the 21 cm signal model, introduces the noise wave formulation (some derivation given in Appendix A), and the calibration scheme based on it, and establishes a framework for cosmological spectrum measurement system simulation, with an emphasis on measurement uncertainties from the Vector Network Analyzer (VNA) (details given in Appendix B) and its impact on cosmological spectrum measurement. Section 3 presents results of a simulation of the global spectrum measurement system, adding VNA measurement uncertainties and other generic systematic errors. How different errors impact the signal recovery are then discussed. Section 4 summarizes our conclusions and discusses future work.

2. Models and Simulation Methods

2.1. Foreground and 21 cm Spectrum Model

The Galactic foreground radiation comes from a variety of complex physical processes [29]. At low frequency, the synchrotron radiation from cosmic ray electrons spiraling in the Galactic magnetic field dominates. In this paper, following Bowman et al. [19], the foreground received by ground based experiment is modeled using polynomials of frequency as:

$$\begin{array}{cc}\hfill T_F\left(\nu \right)\approx & a_0{\left(\frac{\nu }{{\nu }_c}\right)}^{-2.5}+a_1{\left(\frac{\nu }{{\nu }_c}\right)}^{-2.5}log\left(\frac{\nu }{{\nu }_c}\right)+\hfill \\ \hfill & a_2{\left(\frac{\nu }{{\nu }_c}\right)}^{-2.5}{\left[log\left(\frac{\nu }{{\nu }_c}\right)\right]}^2+a_3{\left(\frac{\nu }{{\nu }_c}\right)}^{-4.5}+a_4{\left(\frac{\nu }{{\nu }_c}\right)}^{-2}\hfill \end{array}$$

(1)

where $T_F\left(\nu \right)$ is the brightness temperature of the foreground emission, $\nu$ is the frequency, ${\nu }_c$ is the centre frequency of the observed band and the coefficients $a_n$ are fitted to the data.

The brightness temperature of the redshifted 21 cm signal from the early Universe is frequency dependent [1]. The predicted spectral signature is a relatively broadband absorption (during dark ages and cosmic dawn) and emission (during EoR) signal from 50-200 MHz ($30\le z\le 6$), with a peak absolute absorption amplitude between 10 and 250 mK, which is dependent on the particular model of the first stars/galaxies [30]. In this paper, the 21 cm absorption profile during the cosmic dawn is modeled as a negative flattened Gaussian format [19]:

$$T_{21}\left(\upsilon \right)=-A\left(\frac{1-e^{-{\tau }^B}}{1-e^{-\tau }}\right)$$

(2)

where A is the absorption amplitude, and

$$B=\frac{4{(\upsilon -{\upsilon }_0)}^2}{w^2}log\left[-\frac{1}{\tau }log\left(\frac{1+e^{-\tau }}{2}\right)\right]$$

(3)

${\upsilon }_0$ is the centre frequency, w is the full-width at half-maximum (FWHM) and $\tau$ is a flattening factor. This is a mathematical description of the line profile, not a physical model. We adopt the EDGES measured 21 cm absorption spectrum as our 21 cm spectrum fiducial model. However, we argue that the results of this paper, i.e. the impact of systematic error, is not sensitive to a particular 21 cm signal model or foreground model. The input value of the parameters of the foreground and 21 cm absorption signal model used in this work are listed in Table 1. The foreground, 21 cm absorption, and the combined foreground and 21 cm absorption signal are shown as Figure 1.

2.2. The Noise Wave Formulation

Unlike the artificially generated signals which often have regular wave forms, the naturally produced astronomical signals often appear to be varying randomly and can be well characterized by its power or temperature. In global spectrum measurement instruments, where the radio waves are converted to electrical signals, the astronomical signals are joined by the noise generated within the instruments. At the interfaces of different components of the electronic system, part of the signal propagates down to the next component, and part of the signal is reflected back if the impedance of the two parts are not perfectly matched. Moreover, the receiver is an active device, it does not only receive the signal passively at its front end, but also emits noises.

The noise wave concept [24] provides a description of the propagation of such noise-like signal in the electronic system. This formulation describes the noise behavior of linear RF devices with a few measurable parameters of the system, such as the reflection coefficients and physical temperatures of the calibrators. In this formulation, we quantify the power of the electrical signal by a noise temperature. As the signal propagates through the system and gets reflected back and forth at the junctions of the various components of the system due to impedance mismatches, the noise power and temperature are accordingly transferred. The signal power at the receiver as induced by the signal from the antenna can be written as $P=k{T}_{\mathrm{ant}}\Delta \nu$, where k is the Boltzmann constant, $\Delta \nu$ is the bandwidth, and $T_{\mathrm{ant}}$ is the so called antenna temperature, which is related to the averaged sky temperature $T_{\mathrm{sky}}$ by:

$$\begin{array}{cc}\hfill T_{\mathrm{ant}}& =T_{\mathrm{sky}}(1-|{\mathsf{\Gamma }}_{\mathrm{ant}}{|}^2\left)\right|F_{\mathrm{ant}}{|}^2+T_u|{\mathsf{\Gamma }}_{\mathrm{ant}}{|}^2{\left|F_{\mathrm{ant}}\right|}^2+\hfill \\ \hfill & (T_{c}cos\left(\varphi \right)+T_{s}sin\left(\varphi \right))|{\mathsf{\Gamma }}_{\mathrm{ant}}\left|\right|F_{\mathrm{ant}}|+T_0,\hfill \end{array}$$

(4)

where

$$F_{\mathrm{ant}}=\frac{(1-|{\mathsf{\Gamma }}_{\mathrm{rec}}{{|}^2)}^{\frac{1}{2}}}{1-{\mathsf{\Gamma }}_{\mathrm{ant}}{\mathsf{\Gamma }}_{\mathrm{rec}}},$$

(5)

$${\mathsf{\Gamma }}_{\mathrm{ant}}=\frac{Z_{\mathrm{ant}}-Z_0}{Z_{\mathrm{ant}}+Z_0},$$

(6)

$${\mathsf{\Gamma }}_{\mathrm{rec}}=\frac{Z_{\mathrm{rec}}-Z_0}{Z_{\mathrm{rec}}+Z_0},$$

(7)

Here $Z_0$ denotes the characteristic impedance of transmission line and can be viewed as a reference impedance, and $T_0$ represents the receiver noise offset. ${\mathsf{\Gamma }}_{\mathrm{ant}}$ and ${\mathsf{\Gamma }}_{\mathrm{rec}}$ are the reflection coefficients of the antenna and receiver respectively, and $F_{\mathrm{ant}}$ is the transfer function of the noise wave reflected between the antenna and the receiver. The noise reflected from the antenna reenters the receiver with phase $\varphi$, which is equal to the phase of ${\mathsf{\Gamma }}_{\mathrm{ant}}{F}_{\mathrm{ant}}$. $T_u$ is the uncorrelated portion of the receiver noise reflected back from the antenna, and $T_c$ and $T_s$ are the cosine and sine components of the correlated receiver noise reflected back by the antenna. The derivation of this formula is given in Appendix A.

The reflection coefficient in this formula depends on the impedance of the source and the network connected to the source. For a typical two port network with finite impedance, $|{\mathsf{\Gamma }}_s|=1$ for an open circuit source ($Z=\infty$) or a shorted source ($Z=0\mathsf{\Omega }$), while $|{\mathsf{\Gamma }}_s|=0$ for a matched load with the same impedance as the network. $T_u$, $T_c$, $T_s$ and $\varphi$ can be determined from measurements.

2.3. The Global Spectrum Measurement System

The block diagram of the global spectrum measurement instrument analyzed in this paper is shown in Figure 2. This system can be switched to one of three inputs: an antenna for observation of the sky signal, an open cable, and a 50 $\mathsf{\Omega }$ resistance load, the latter two inputs are used for system calibration. The receiver is modeled as a two-port network system. We use a method similar to the one used by the EDGES experiment to do calibration and signal reconstruction. Here we use a 50 $\mathsf{\Omega }$ ambient resistor load to calibrate the system response, and assume a stable and precise 300 K ambient temperature load for absolute temperature calibration. We note that both an ambient and a hot calibrator are used in the EDGES experiment, which can provide calibration for two parameters–the frequency-dependent temperature scale factor and the offset. In the present work we focus on the impact of the VNA measurement error, so here we only use an ambient load for temperature reference, which simplifies the description, and does not significantly affect our estimates. In actual experiments, however, it may be necessary to introduce both the ambient and hot calibrators.

An open cable is used to solve the noise wave parameters in this simulation. Here we adopted a simple model of the system, with $T_u,T_c,T_s$ as linear functions of frequency. Results show that the fitting residue is less than 3 mK when using only an open cable to solve the noise wave parameters. In a real experiment, the system could be more complicated, e.g. the system temperature could have a more complicated dependence on the frequency. An additional shorted cable can provide an independent calibration measurement, and thus help improving the accuracy of the noise wave parameter determination.

We focus on how the systematic errors propagate in the global spectrum measurement system, and their impact on the recovered signal. We first consider the time-independent case, the time variation of the system parameters such as the gain, the noise excess temperature and the ambient temperature are discussed in latter sections.

When an antenna is connected to the receiver, the sky signal is fed into the system, but if the antenna output impedance and the receiver input impedance are not in perfect match, a fraction of the signal power would be reflected back to the antenna. Moreover, the receiver also produces noise which propagates back to the antenna. The receiver noise can be separated into a component which is correlated with the receiver output, and a component which is uncorrelated [25]. The uncorrelated noise power depends only on the magnitude of the antenna reflection, while the correlated noise power depends also on the phase. The antenna temperature is related to the sky noise power as given by Eq. (4).

For receiver calibration, the receiver input can also be switched to a 50 $\mathsf{\Omega }$ ambient resistor load or an open cable. With the 50 $\mathsf{\Omega }$ resistor load as input, the receiver noise temperature can be written as:

$$\begin{array}{cc}\hfill T_{\mathrm{load}}& =T_{\mathrm{amb}}(1-|{\mathsf{\Gamma }}_{\mathrm{load}}{|}^2\left)\right|F_{\mathrm{load}}{|}^2+T_u|{\mathsf{\Gamma }}_{\mathrm{load}}{|}^2{\left|F_{\mathrm{load}}\right|}^2+\hfill \\ \hfill & (T_{c}cos\left(\varphi \right)+T_{s}sin\left(\varphi \right))|{\mathsf{\Gamma }}_{\mathrm{load}}\left|\right|F_{\mathrm{load}}|+T_0,\hfill \end{array}$$

(8)

where

$$F_{\mathrm{load}}=\frac{(1-|{\mathsf{\Gamma }}_{\mathrm{rec}}{{|}^2)}^{\frac{1}{2}}}{1-{\mathsf{\Gamma }}_{\mathrm{load}}{\mathsf{\Gamma }}_{\mathrm{rec}}},$$

(9)

$${\mathsf{\Gamma }}_{\mathrm{load}}=\frac{Z_{\mathrm{load}}-Z_0}{Z_{\mathrm{load}}+Z_0},$$

(10)

${\mathsf{\Gamma }}_{\mathrm{rec}}$ is given by Eq (7) since the same receiver is used, $\varphi$ is the phase of ${\mathsf{\Gamma }}_{\mathrm{load}}{F}_{\mathrm{load}}$, and $T_{\mathrm{amb}}$ is the ambient temperature of the resistance load. For a realistic experiment system, $Z_{\mathrm{load}}\approx 50\mathsf{\Omega }$, then ${\mathsf{\Gamma }}_{\mathrm{load}}\approx 0$. Equation (8) is reduced to

$$\begin{array}{c}\hfill T_{\mathrm{load}}=T_{\mathrm{amb}}(1-|{\mathsf{\Gamma }}_{\mathrm{rec}}{|}^2)+T_0\end{array}$$

(11)

Similarly, when connected to the open cable, the receiver noise temperature is given by:

$$\begin{array}{cc}\hfill T_{\mathrm{open}}& =T_{\mathrm{amb}}(1-|{\mathsf{\Gamma }}_{\mathrm{open}}{|}^2\left)\right|F_{\mathrm{open}}{|}^2+T_u|{\mathsf{\Gamma }}_{\mathrm{open}}{|}^2{\left|F_{\mathrm{open}}\right|}^2+\hfill \\ \hfill & (T_{c}cos\left(\varphi \right)+T_{s}sin\left(\varphi \right))|{\mathsf{\Gamma }}_{\mathrm{open}}\left|\right|F_{\mathrm{open}}|+T_0,\hfill \end{array}$$

(12)

where

$$F_{\mathrm{open}}=\frac{(1-|{\mathsf{\Gamma }}_{\mathrm{rec}}{{|}^2)}^{\frac{1}{2}}}{1-{\mathsf{\Gamma }}_{\mathrm{open}}{\mathsf{\Gamma }}_{\mathrm{rec}}},$$

(13)

$${\mathsf{\Gamma }}_{\mathrm{open}}=L_c\frac{Z_{\mathrm{in}}-1}{Z_{\mathrm{in}}+1},$$

(14)

where $L_c$ is the cable loss factor $L_c={10}^{-\frac{l_c}{10}}$ , $l_c$ is the one-way cable loss, it is a function of frequency $\nu$. For open cable, $Z_{\mathrm{in}}=-j{Z}_{0}cot\beta l$. The cable acts just like an antenna looking at an isotropic sky with temperature equal to the physical temperature of the cable. The noise wave parameters $T_u$, $T_c$ and $T_s$ can be solved from these calibration measurements.

In Figure 3 we plot the block diagram of the sky signal reconstruction procedure for this global spectrum measurement system. The 21 cm spectrum signal is recovered from the observation and calibration mock data in 4 steps: (1) Calibration for receiver response using the load calibration (mock) data, the receiver reflection coefficient measurement data, with assumed ambient temperature for the load, and receiver noise temperature. (2) Recover the open cable temperature $T_{\mathrm{open}}$, and fitting it to the noise wave formulation with the reflection characteristics of the open cable and receiver, then solve for noise wave parameters $T_u,T_c,T_s$. (3) Using the system parameters obtained in the first two steps to recover antenna temperature from the raw data, then reconstruct the sky temperature. (4) Fitting the sky temperature with foreground model only, and with both foreground and 21 cm model simultaneously, then derive the foreground and 21 cm model parameters.

2.4. Hardware Models

We assume the antenna to be a wideband blade dipole similar to the one used by the EDGES experiment, though the exact parameters may differ. The antenna response is obtained by an electromagnetic field simulation using the CST STUDIO SUITE ([31]) software.

To obtain concrete estimates of the noise, we assume the following design and components. The custom-made Low Noise Amplifier (LNA) receiver covers the 30-250 MHz band, and is optimized for low noise figure and flatness of $S_{11}$ and $S_{22}$. The Wantcominc WHM0003AE [32] transistor is used for the input stage as well as the second amplification stage. This is a wideband, high linearity SMT packaged amplifier with exceptional gain flatness design. The gain and noise figure of the LNA are plotted in Figure 4 as a function of frequency. The noise figure is typically 0.70 dB. Measurements of the LNA show an average $S_{11}$ of -21.6 dB, with a flatness better than 0.4 dB within the band, and the average $S_{22}$ is -20.0 dB, with a flatness better than 0.1 dB across 50-100 MHz. The LNA offers typically 40 dB gain with a flatness of 0.05 dB, and the noise figure is less than 0.95 dB across the band 50-100 MHz.

An open-terminated cable of 5 meters length is used for internal calibration of the system. The open cable model is based on a precision test cable assembly, its impedance is 49.6 $\mathsf{\Omega }$, and the velocity of propagation in the cable is 83% of speed of light, a first-order polynomial of frequency is used as the attenuation coefficient, i.e. 0.24 dB/m @ 50 MHz, 0.30 dB/m @ 100 MHz. The load used for receiver calibration is modeled as a pure impedance terminator, its impedance is set as 50 $\mathsf{\Omega }$.

The reflection coefficients ($S_{11}$) of the antenna, the LNA receiver, and the open cable are shown in Figure  as a function of frequency. These are complex quantities so both the magnitudes and phases are shown, and the expected measurement errors (to be discussed in the next section) are also plotted.

2.5. Vna Measurement Error Model

A Vector Network Analyzer (VNA) can be used to measure the S parameters of a one-port or two-port system. In a 21 cm global spectrum measurement, a VNA is employed to measure the reflection coefficients [25]. The measurement accuracy of the reflection coefficients ($S_{11}$) of antenna, cable, receiver are very important for our purpose. As a measurement instrument, VNA has errors that are caused by imperfections within the network analyzer, these errors are referred to as measurement errors. Measurement errors can be classified into three groups: Systematic, Random, Drift & Stability, more detailed description of the source of VNA error is given in Appendix B.

For a given set of VNA, standard calibrator & cable, the measurement error or uncertainty in the reflection coefficient ($S_{11}$) can be estimated using the VNA error model as shown in Figure A2, with the magnitude/phase given in Eq. (A18) and Eq. (A25). Below, we take the Keysight^TM N5247B PNA-X series Vector Network Analyzer and the corresponding 85050C precision mechanical calibration kit as the components used in our model. Table 2 list typical specifications for directivity, source match, load match, reflection tracking, and dynamic accuracy of this type of VNA. These specifications are for typical measurement setups: 10 Hz IF bandwidth, isolation calibration with an averaging factor of 8, and the specifications of the systematic effects are guaranteed values and they should be the maximum values in the frequency domain. These are estimated with the Vector Network Analyzer Uncertainty Calculator script supplied by its maker [33].

Figure 6 shows the calculated magnitude and phase uncertainties of the reflection coefficient measurement for the specific VNA and the calibrator standard set. Below, we use these results as the VNA measurement error in our simulation. Distribution of the systematic error is usually unknown, and the uniform distribution is assumed. The standard deviation is calculated as $\left(\mathrm{maximum}\mathrm{value}\right)/\sqrt{3}$, and the total uncertainty is estimated by taking the root sum square (RSS) of the systematic and random errors, though the errors introduced through cable movement or temperature drift are not included. As shown in the figure, the uncertainty in the phase is large when the magnitude of the reflection coefficient is small, because in that case the signal is too weak to be measured precisely.

3. Error Propagation

We now consider how the errors affect the measured signal. We first consider the errors arising from the receiving system, which we shall call receiver errors, and then those induced by the VNA measurement errors, which we shall call VNA errors. We use the reconstruction simulation pipeline to simulate and analyze the impact of each error. We also consider the application of attenuator in the system, which in some cases could reduce the system error.

3.1. Receiver Error Simulation

The noise figure of the receiver in our example system is about 1 dB in the observation band (c.f. Figure 4), which corresponds to a 75 K receiver noise temperature for a load physical temperature of 290 K, while the transmission line loss between the antenna and the receiver is negligible. The input signal is amplified by the receiver, and then fed into a digital data acquisition system, here we neglect the errors in the digital system. Inspired by actual experimental results, we model the noise wave parameters $T_u,T_c,T_s$ as linear functions of frequency, $T_u=0.04{\nu }_{\mathrm{MHz}}+31\mathrm{K}$, $T_c=0.04{\nu }_{\mathrm{MHz}}+6\mathrm{K}$, $T_s=0.06{\nu }_{\mathrm{MHz}}+6\mathrm{K}$, where ${\nu }_{\mathrm{MHz}}$ is the frequency in units of MHz. We assume the sampling rate is 100 MSPS for 50 MHz bandwidth, the FFT length is 16384, to produce 8192 output spectral bins. A total of 500,000 seconds (138.89 hours) of mock data for antenna input and the same amount of mock data for the load and open cable inputs are generated. Thermal noise is also added in each terminator, with $\delta T=(T+T_{\mathrm{rec}})/\sqrt{\Delta \nu \tau }$, where $\tau$ is integration time (set to 10 seconds), $\Delta \nu$ is bandwidth, $T_{\mathrm{rec}}$ is receiver noise temperature, and T is the temperature of the terminator, i.e. when switched to antenna $T=T_{\mathrm{sky}}$, when switched to open cable or load, T equals to the ambient temperature. The raw spectrum accumulated from each input is averaged, which reduces random noise. Figure 7 shows the simulated receiver output for the three input sources with an integration period of 10 seconds. The data for the open cable shows an oscillation with respect to the frequency, as the reflection at the end of the cable set up a standing wave in the system. If the reflection coefficients are known, we can then derive the noise wave parameters $T_u,T_c,T_s$ as a function of frequency by least square fitting of the calibration data.

We plot the calibrated sky spectrum in the top panel of Figure 8. To reduce scatter, the data of every 64 frequency channels are rebinned to one, the resulting spectrum has 128 frequency bins with a width of 393.7 kHz. This is fitted with the five term foreground model given by Eq. (1) in the middle panel, and the foreground plus 21 cm signal (Eq. 2) model in the bottom panel. As expected, in the foreground-only model the residue appears to have spectral features, while the foreground+21 cm model yields flat residues. Figure 9 shows the distributions of the foreground model parameters ($a_0$, $a_1$, $a_2$, $a_3$, $a_4$) and the 21 cm model parameters (A, ${\upsilon }_0$, w, $\tau$) derived from the mock observation. Although there is some degeneracy in the foreground polynomial coefficients, the 21 cm model parameters are basically uncorrelated amongst themselves.

In the above we have assumed that the system is time-independent, ignoring time variation of the system characteristics. However, the gain of the system will drift over time, especially if there is temperature variations in the circuit. If the overall gain drifts ±0.5 dB, while keeping other system parameters stable and bandpass shape unchanged, simulations shows that it cause about 20 mK deviation in the result for our system. The gain could be relatively calibrated with a noise source (usually a diode) relatively, and absolutely with a thermal load. The calibration error is given by $\delta g=\delta T_{\mathrm{cal}}/T_{\mathrm{cal}}$. If we assume the noise source temperature fluctuation is given by $\delta T_{\mathrm{cal}}=T_{\mathrm{cal}}/\sqrt{\Delta \nu \tau }$, then for an integration time of 10 seconds and bandwidth of 393.7 kHz, we have $\delta g=1/\sqrt{\Delta \nu \tau }\sim 5\times {10}^{-3}$ or 0.002 dB.

The output of the load in the circuit also has thermal fluctuations. This is used to calibrate the $T_0$ by applying Eq. 11, so the fluctuation would induce an error of the size $\delta T_0\sim T_{\mathrm{amb}}/\sqrt{\Delta \nu \tau }$. We set the integration time to be 10 seconds as before, while we use three cycle’s thermal load data (30 seconds) for receiver calibration. This will induce a variation in the calibrated gain of the system. Figure 10 shows the probability distribution of the deviation from the true signal in our simulation, indicating an error with a mean value of 45 mK. This error could be further reduced if longer integration time is used. An optimal time scale could be achieved by equalizing the thermal noise and the error induced by the drift, the latter depends on the stability of the system.

3.2. Vna Error Simulation

The system calibration depends on the VNA measurement results. However, even with the finest test equipment, there are still residue errors, as shown in Figure 6, which could bias the result. We now simulate the impact of the VNA measurement. The model is the same as above but without thermal noise.

For the reflection coefficients of the antenna, the load and the open cable, we assume the VNA measurements have errors as given by the VNA specification and shown above in Figure 6. However, such specifications are often not given in terms of standard deviations, to be on the safe side, the uniform distribution is assumed. The standard deviation is calculated as $\left(\mathrm{maximum}\mathrm{value}\right)/\sqrt{3}$. Such error is also likely correlated over the frequencies, though the exact variation is unknown. We have plotted the magnitude and phase of the antenna, receiver, and the open cable reflection coefficients together with their error bars in Figure , computed as prescribed in Section. 2. In each panel, the solid and dashed curves show the magnitude (read number from vertical axis on the left) and phase (read number from vertical axis on the right) respectively, and the uncertainty is shown in error bars. Note here the scales are different for each panel, as in each case the variation with frequency is different.

The sky temperature can be derived from the observed quantities as:

$$\begin{array}{cc}\hfill T_{\mathrm{sky}}=& \frac{T_{\mathrm{ant}}-T_0}{(1-{\left|{\mathsf{\Gamma }}_{\mathrm{ant}}\right|}^2){\left|F_{\mathrm{ant}}\right|}^2}-\frac{T_u{\left|{\mathsf{\Gamma }}_{\mathrm{ant}}\right|}^2}{1-{\left|{\mathsf{\Gamma }}_{\mathrm{ant}}\right|}^2}-\frac{(T_{c}cos\left(\varphi \right)+T_{s}sin\left(\varphi \right))\left|{\mathsf{\Gamma }}_{\mathrm{ant}}\right|}{(1-{\left|{\mathsf{\Gamma }}_{\mathrm{ant}}\right|}^2)\left|F_{\mathrm{ant}}\right|}\hfill \end{array}$$

(15)

The errors in each quantity will propagate into the sky temperature, and affect the reconstruction of the 21 cm temperature. We can use simulations to illustrate the magnitude and pattern of each source of error.

The error induced in the sky temperature by measurement errors can be estimated by the error propagation formula:

$$\begin{array}{c}\hfill \delta T_{\mathrm{sky}}=\frac{\partial T_{\mathrm{sky}}}{\partial {\mathsf{\Gamma }}_{\mathrm{ant}}}\delta {\mathsf{\Gamma }}_{\mathrm{ant}}+\frac{\partial T_{\mathrm{sky}}}{\partial {\mathsf{\Gamma }}_{\mathrm{rec}}}\delta {\mathsf{\Gamma }}_{\mathrm{rec}}+\frac{\partial T_{\mathrm{sky}}}{\partial {\mathsf{\Gamma }}_{\mathrm{open}}}\delta {\mathsf{\Gamma }}_{\mathrm{open}}\end{array}$$

(16)

where $\delta {\mathsf{\Gamma }}_{\mathrm{ant}}$, $\delta {\mathsf{\Gamma }}_{\mathrm{rec}}$ and $\delta {\mathsf{\Gamma }}_{\mathrm{open}}$ denotes the VNA measurement error on the antenna, receiver and open cable respectively. ${\mathsf{\Gamma }}_{\mathrm{ant}}$ is present in the expression of sky temperature Eq. (15) (also through the F term, c.f. Eq. (5)) so the derivative can be computed directly. The sky temperature is not directly dependent on ${\mathsf{\Gamma }}_{\mathrm{open}}$, but $T_u,T_c,T_s$ are derived from the calibration measurement given in Eq.(12), so

$$\begin{array}{c}\hfill \frac{\partial T_{\mathrm{sky}}}{\partial {\mathsf{\Gamma }}_{\mathrm{open}}}=\frac{\partial T_{\mathrm{sky}}}{\partial T_u}\frac{\partial T_u}{\partial {\mathsf{\Gamma }}_{\mathrm{open}}}+\frac{\partial T_{\mathrm{sky}}}{\partial T_c}\frac{\partial T_c}{\partial {\mathsf{\Gamma }}_{\mathrm{open}}}+\frac{\partial T_{\mathrm{sky}}}{\partial T_s}\frac{\partial T_s}{\partial {\mathsf{\Gamma }}_{\mathrm{open}}}\end{array}$$

(17)

As $T_u,T_c,T_s$ are derived from a numerical fitting procedure of the measured $T_{\mathrm{open}}\left(\nu \right)$ over a frequency range, we obtain the partial derivatives numerically. Finally, ${\mathsf{\Gamma }}_{\mathrm{rec}}$ is present in the sky expression directly through $F_{\mathrm{ant}}$ (Eq. (5)), and indirectly by affecting $T_u,T_c,T_s$ (Eq.(12)), and $T_0$ (through Eq.11). So we have

$$\begin{array}{ccc}\hfill \frac{\partial T_{\mathrm{sky}}}{\partial {\mathsf{\Gamma }}_{\mathrm{rec}}}& =& {\left.\frac{\partial T_{\mathrm{sky}}}{\partial {\mathsf{\Gamma }}_{\mathrm{rec}}}\right|}_{\mathrm{direct}}+\frac{\partial T_{\mathrm{sky}}}{\partial T_u}\frac{\partial T_u}{\partial {\mathsf{\Gamma }}_{\mathrm{rec}}}+\frac{\partial T_{\mathrm{sky}}}{\partial T_c}\frac{\partial T_c}{\partial {\mathsf{\Gamma }}_{\mathrm{rec}}}+\frac{\partial T_{\mathrm{sky}}}{\partial T_s}\frac{\partial T_s}{\partial {\mathsf{\Gamma }}_{\mathrm{rec}}}+\frac{\partial T_{\mathrm{sky}}}{\partial T_0}\frac{\partial T_0}{\partial {\mathsf{\Gamma }}_{\mathrm{rec}}}\hfill \end{array}$$

(18)

where the first term on the R.H.S. of the equation is obtained by calculating the partial derivative using the explicit function form directly, while the remaining terms are calculated in a way similar to the ${\mathsf{\Gamma }}_{\mathrm{open}}$ case.

We estimate the VNA measurement errors according to the value of reflection coefficients for the antenna, receiver and open cable as shown. The error in magnitude is at the level of a few times ${10}^{-3}$, and the error in the phase is at the level of $1.5^{\circ }$. In the frequency range of 50 MHz-100 MHz, these errors only have slight variations.

In Figure 11a–c we plot separately the deviation in reconstructed sky temperature caused by the VNA error when it is used to measure the reflection coefficients of the antenna, receiver and open cable respectively. As the errors given in the manual of VNA are $\left(\mathrm{maximum}\mathrm{values}\right)/\sqrt{3}$ and are perhaps correlated over the frequencies, we plot the induced errors with the error set in the range of 0.1 to 1.0 times the nominal uncertainty, which vary smoothly over the frequencies, and the multiple curves remind us the possible size of the error. From this figure, we see the induced errors are large for the antenna and receiver measurements, while the open cable measurement error has less impact on results.

In Figure 12a we plot the recovered sky temperature deviation with the VNA measurement error on antenna, receiver and open cable simultaneously–this is closer to the real case as the same VNA are switched to make these measurements. Results show that for an error equal to 10% (0.1) of the nominal error, there is a deviation of about 0.5 K in reconstructed sky temperature, while for the nominal value the deviation is about 6 K in sky temperature at 63 MHz.

The sky temperature spectrum is fitted with the foreground and 21 cm model spectra, so the error in the reconstructed 21 cm signal is a nonlinear function of the VNA measurement error which does not have an analytical expression. However, one can apply the fitting algorithm to the sky temperature spectrum with the error, and then derive the residue in 21 cm signal reconstruction. We plot in Figure 12b the 21 cm reconstruction results. The error in the reconstructed 21 cm signal is not directly proportional to the whole deviation in sky temperature spectrum, as the smooth components are removed as foregrounds. The residue deviations in the 21 cm spectrum have much smaller magnitude than the whole sky temperature spectrum shown in Figure 12a, but are less smooth, showing more oscillatory features, and the deviation magnitude also increases as the VNA error grows. The deviation level can reach $\pm \sim 750$ mK maximum in this model. We also calculated the RMS deviation value of reconstructed 21 cm signal for each magnitude uncertainty level, e.g. for half (0.5) the maximum uncertainty in magnitude, the mean uncertainties of the reflection coefficient of the antenna, the receiver, and the open cable over the 50 MHz-100 MHz band are about 0.001, 0.001, 0.002 (linear) respectively, and the RMS deviation value of reconstructed 21 cm signal is 200 mK.

The global spectrum experiment considered in the present paper measures the total received power from one antenna, the phase information of the wave is not recorded. Despite this, the phase error in the reflection coefficient will also affect the total power spectrum reconstruction. In the above we have considered the error induced by the magnitude of the VNA measurement. Next we consider the effect of the error of VNA phase measurement. The phase measurement error also depends on the reflection coefficient. As shown in Figure 6, for a reflection coefficient of larger than 0.4 ($S_{11}=-7.95\mathrm{dB}$), the phase uncertainty curve is quite flat. In the simulation, we set a series of 10 relative error values, for example, the uncertainty for open cable is ranging from the relatively small $0.{03}^{\circ }$, to the maximum uncertainty of $0.3^{\circ }$, while for receiver, the uncertainty is ranging from the relatively small $0.{15}^{\circ }$, to the maximum uncertainty of $1.5^{\circ }$. Figure 13 shows the reconstruction results. As the error grows the deviation become greater, while the absorption shape of the reconstruction result is similar to the residual model. We also calculated the RMS deviation value of the reconstructed 21 cm signal for each phase uncertainty level, e.g. for a half (0.5) of the maximum phase uncertainty, the mean uncertainties of the antenna, the receiver, and the open cable over the 50 MHz-100 MHz band are about $0.{48}^{\circ }$, $0.{78}^{\circ }$, $0.{15}^{\circ }$ respectively, and the RMS deviation value of reconstructed 21 cm signal is 40 mK.

To evaluate whether the length of the calibrator open cable has an impact on the calibration accuracy and 21 cm signal reconstruction when VNA measurement errors are added in, we also simulated the cases with open cables of length from 3 meters to 30 meters. The result for the case of a 15 meters open cable is plotted in Figure 14. The errors on the magnitude and phase are added to the terminators’ reflection characteristics, then the same simulation and reconstruction procedures are applied as before. It seems that there is very little difference in the result, showing that the cable length does not have a significant impact on the precision of the experiment under the premise that the noise wave parameters were solved with high accuracy.

3.3. Effect of an Attenuator

In some global spectrum setups, an attenuator is inserted between the antenna and receiver, to improve the effective matching between the antenna and the receiver, reducing the standing wave amplitude, and alleviating the impact of VNA measurement uncertainty. However, adding a passive device before the LNA will also increase the thermal noise of the system.

Below we treat the antenna terminated with the attenuator as an equivalent antenna, as depicted in Figure 15. In this approach, the parameters associated with the receiver, including $T_u$, $T_c$, $T_s$ are unchanged, just as it is replaced by the open cable during calibration. The reflection coefficient of the effective antenna, i.e. the combination of antenna and attenuator is

$${\mathsf{\Gamma }}_{\mathrm{a}}^{\prime }=S_{22}+\frac{S_{12}{S}_{21}{\mathsf{\Gamma }}_{\mathrm{a}}}{1-S_{11}{\mathsf{\Gamma }}_{\mathrm{a}}}$$

(19)

where the S-parameters are those of the attenuator. For an ideal attenuator, $S_{11}=S_{22}=0,S_{21}=S_{12}$, $G=S_{21}^2$, then

$$\begin{array}{c}\hfill {\mathsf{\Gamma }}_{\mathrm{a}}^{\prime }=G{\mathsf{\Gamma }}_{\mathrm{a}}.\end{array}$$

(20)

Compared with the original reflection coefficients of the antenna, the magnitude of the reflection coefficient is reduced, as shown in Figure 16, with attenuation of $1\sim 10$ dB.

Figure 17 shows the sky temperature deviation and 21 cm signal deviation for the equivalent antenna and the equivalent receiver model respectively. Here we add both magnitude and phase measurement error(1.0× uncertainty) in the antenna, the receiver and the open cable. Results show that with the attenuator, the residue error is indeed reduced. For example, the maximum residue is reduced from 0.63 K to 0.44 K for 21 cm signal for attenuation of 9 dB.

The original $T_0$ corresponds to the output noise temperature generated by the LNA, which is 75 K for a load with physical temperature of 290 K. In a global spectrum experiment, the thermal noise in each terminator can be expressed as $\delta T=(T+T_{\mathrm{rec}})/\sqrt{\Delta \nu \tau }$, so an increased of receiver noise temperature will result in the increase of thermal noise in the system, which will need a longer integration time to achieve the desired signal to noise ratio.

Figure 18 shows the recovered residue after foreground subtraction for without attenuator, and cases of 3 dB, 6 dB, 9 dB attenuation. To keep the RMS of the residue at 0.05 K level, different integration times are required. Results show that for the 0 dB, 3 dB, 6 dB and 9 dB attenuator cases, the integration time required are 138.89, 277.78, 694.44,1111.11 hours respectively.

From the above analysis, we see the addition of an attenuator between the antenna and LNA could reduce the impact of VNA measurement uncertainty, but increase the thermal noise. To get the best performance, a compromise should be made by choosing a moderate value of attenuation. For the experimental setup considered here, it seems that an attenuation of ∼ 3 dB is near the optimal, while adopting a greater attenuation would only reduced the deviation slightly, but would require considerably longer integration time.

3.4. Discussions

We estimate the induced errors on the final 21cm signal. The central frequency of 21cm signal ${\nu }_0$ is not influenced much by the pure VNA magnitude error, while the width of the absorption profile W is not influenced much by the pure phase error. However, the typical VNA magnitude measurement errors of 0.001, 0.001 and 0.002 (linear) on antenna, receiver and open cable measurements respectively could result in a 200 mK deviation in the 21cm absorption depth A, and a 2 MHz deviation on the absorption profile width W. The typical phase measurement errors of $0.{48}^{\circ }$, $0.{78}^{\circ }$, $0.{15}^{\circ }$ on the antenna, receiver, and open cable measurements respectively will result in a 40 mK deviation on absorption depth and a 2 MHz deviation on the central frequency.

Since $A\sim 30(1-\frac{T_R}{T_S})$ mK, where $T_R$ and $T_S$ are the radiation field and spin temperature, suppose $T_S=7$ K which is the adiabatic-cooling temperature at $z\sim 17$ for the gas [34,35,36,37], then $\delta A=200$ mK (40 mK) is equivalent to $\delta T_R\sim 50$ K (10 K), or $\delta T_S\sim 7$ K (1.5 K). Therefore, the deviation on the 21cm absorption depth caused by the VNA error, if not properly corrected in data reduction, might confuse the astrophysical interpretation of the observed 21cm signal, e.g., whether extra radio background or exotic cooling mechanism is essential, and/or whether weak X-ray heating and strong Ly$\alpha$ coupling is inevitable.

We foresee that techniques could be developed to improve the VNA measurement accuracy, so that it would not adversely affect the 21 cm global spectrum experiment. For example, at present partly the error given above is for an instrument calibrated with a standard calibration kit, and the parameters for the kit are provided by the manufacturer. More precise measurement of the individual calibrator kit can be made in the laboratory, especially by characterizing and fitting the impedance and the voltage reflection coefficient of as functions of frequency. The error in the phase of the VNA measurement could be reduced by laboratory calibration with an air-dielectric line which has more precisely known delay characteristics [38]. The VNA error model can also be refined, so that a number of parameters could be individually determined by additional calibration measurements [39]. The temperature variation of the instrument also could be compensated. In each calibration measurement, beyond what is well understood by the physical model, the residue could also be checked and reduced by empirical methods, e.g. the ripple method [40].

4. Conclusions

Detecting the 21 cm global spectrum signal from the cosmic dawn and the epoch of reionisation is a very challenging experiment, as it requires extremely high measurement precision. In this paper we present a simulation of the observation and reconstruction of a 21 cm global spectrum experiment, mainly focus on systematic error propagation and its impact on recovering the 21 cm signal. We model the foreground with a smooth power-law spectrum, and the 21 cm absorption as a negative flattened Gaussian function centered at 78 MHz, as suggested by the EDGES observation result. Our simulation confirms that thermal noise can be reduced to a low level if the system is stable and allows reasonably long integration time. Furthermore, a stable and constant deviation of the system gain or receiver noise offset temperature only have a slight impact on the recovered signal.

Over the broad frequency band needed for the cosmic dawn 21 cm signal measurement, mismatch in the impedance of the various RF components inevitably generate reflections and standing waves inside the instrument. If not properly accounted for, such features can be confused with the 21 cm signal. Based on the noise wave formulation, this effect can be corrected, but it is necessary to measure the system parameters by performing a high precision calibration. Since the system calibration accuracy is determined by the VNA measurement results of reflection coefficients of the antenna, the receiver and the calibrator terminators, any measurement error due to the imperfection of the instrument may affect the calibration and the 21 cm reconstruction result. In order to evaluate this possible source of error, we modeled the uncertainty of the reflection coefficient measurement according to the VNA signal flow graph, for a typical VNA model and calibrator type with one of the highest accuracy specifications presently available from industry, the model includes all known possible errors in VNA measurement procedure and in calibration standards, and realistic VNA uncertainty values are used. We propagate the error to the deviation in sky temperature, and the extracted 21 cm signal from the foreground removal procedure. The results show that even a relatively small error may cause distortion in both the recovered total sky signal and the 21 cm signal, with amplitude comparable with or even larger than the true 21 cm signal. e.g. a typical VNA magnitude measurement error of 0.001, 0.001, 0.002 (linear) exist in the antenna, the receiver, the open cable respectively would result in a 200 mK deviation on the expected signal. A typical phase measurement error of $0.{48}^{\circ }$, $0.{78}^{\circ }$, $0.{15}^{\circ }$ in the antenna, the receiver, and the open cable respectively would cause a 40 mK deviation. Results also show that an attenuator will benefit to the effect that the VNA measurement uncertainties caused. If this kind of systematic errors can not be reduced in calibration procedures or corrected in data analysis, it may confuse the astrophysical interpretation of observed 21 cm signal. It is of vital importance to develop new methods or techniques to further improve the accuracy of calibration.