Calibration Error in 21-Centimeter Global Spectrum Experiment

Abstract

The redshifted 21 cm line signal is a powerful probe of the cosmic dawn and the epoch of reionization. The global spectrum can potentially be detected with a single antenna and spectrometer. However, this measurement requires an extremely accurate calibration of the instrument to facilitate the separation of the 21 cm signal from the much brighter foregrounds and possible variations in the instrument response. Understanding how the measurement errors propagate in a realistic instrument system and affect system calibration is the focus of this work. We simulate a 21 cm global spectrum observation based on the noise wave calibration scheme. We focus on how measurement errors in reflection coefficients affect the noise temperature and how typical errors impact the recovery of the 21 cm signal, especially in the frequency domain. Results show that for our example set up, a typical vector network analyzer (VNA) measurement error in the magnitude of the reflection coefficients of the antenna, receiver, and open cable, which are 0.001, 0.001, and 0.002 (linear), respectively, would result in a 200 mK deviation on the detected signal, and a typical measurement error of $0.48°$, $0.78°$, or $0.15°$ in the respective phases would cause a 40 mK deviation. The VNA measurement error can greatly affect the result of a 21 cm global spectrum experiment using this calibration technique, and such a feature could be mistaken for or be combined with the 21 cm signal.

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 the 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 a high precision instrument. The experimental efforts include the Experiment to Detect the Global EoR Signature (EDGES, [6,7]), the Broadband Instrument for Global HydrOgen ReioNization 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. There are also experiments using short spacing interferometers, such as the Short spacing Interferometer Telescope probing cosmic dAwn and epoch of ReionizAtion (SITARA) [19]. 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, [20]), Dark Ages Radio Explorer project (DARE, [21]), and the Discovering the Sky at Longest wavelength project (DSL, [22]).

The EDGES experiment has reported the detection of a 500 mK deep absorption feature centered at 78 MHz [23], which may be associated with the signature of cosmic dawn. This absorption strength is much stronger than the prediction of the standard model, so its cosmological interpretation may require new physics or astrophysics mechanisms, e.g., cooling mechanism with exotic dark matter particles [24], or extra radio background [25,26]. However, it has also been questioned whether this feature is true or arises from some systematical errors [27,28], and the measurement taken by the SARAS-3 experiment has not detected such an absorption feature [29].

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 emissions, 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.

To derive an accurate global spectrum, a calibration with very high precision is required. The calibration of an instrument can be based on either the measurement of an external standard source, or by the internal calibration mechanisms designed into the system. For the global spectrum experiment, calibration using an external standard source is difficult, because its antenna has a very wide beam, covering essentially the whole sky above the horizon. Although some crude calibration could be achieved with a model of sky radiation and helped by the variation in sky due to Earth rotation, as, for example, used by the SCI-HI experiment [16], the precision of such calibration is limited by the large uncertainties in the sky model and the beam model. It is also conceivable that one could try to use an artificial external source, but again, this is complicated by the antenna response, which is not isotropic. To mimic the response to the whole sky, the whole apparatus can, in principle, be enclosed inside an anechoic chamber kept at a constant temperature. Ideally, the instrument is thus well-isolated and immersed in the black body radiation of the chamber, which can serve as the standard source for calibration 1. However, it is very challenging to build such an anechoic chamber with sufficient shielding and absorbing power at low frequency and keeping all its parts at the same temperature.

Presently, a more practical approach to high-precision calibration is based on internal calibrators and the noise wave formulation of the electric circuit system [30]. This approach was pioneered by the EDGES group [31] and subsequently adopted by most experiments. 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 using the vector network analyzer (VNA). 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 the measured results.

To assess the uncertainties of the measurement, in addition to the thermal noise, one also needs to consider the possible systematic errors that may bias the result. This is the main aim of the present work. Previously, Monsalve et al. [32] and Monsalve et al. [33] presented the calibration strategy of the EDGES experiment and analyzed the propagation of uncertainties in the receiver calibration. In a recent work [34], 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 VNA measurement error, which depends on frequency, and adopt a comprehensive model for VNA measurement uncertainties, including all measurement errors in the VNA measurement procedure and in calibration standards, to estimate their impacts on the spectrum measurement system. Furthermore, we also study how the error affects the extracted 21 cm signal, i.e., the non-smooth component of the sky spectrum. In the simulation, we have adopted the parameter values we found in our actual experiment work. While there are many possible instrument designs, and many hardware components with different characteristics, the differences are not that large. We have used typical values, although there is still room for refinement.

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 the 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 the 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 [35]. At low frequency, the synchrotron radiation from cosmic ray electrons spiraling in the galactic magnetic field dominates. In this paper, following Bowman et al. [23], the foreground received by the ground-based experiment is modeled using polynomials of frequency, with five terms based on the known spectral properties of the galactic synchrotron spectrum and Earth’s ionosphere [36], 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 center 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 [37]. In this paper, the 21 cm absorption profile during the cosmic dawn is modeled as a negative flattened Gaussian format [23]:

$$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 center 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 description of the physics that create the 21 cm absorption profile. 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, are 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 input value of the foreground and 21 cm signal model parameters.

$a_0$	$a_1$	$a_2$	$a_3$	$a_4$
1284	570	−1240	753	98
A	${\upsilon }_0$	w	$\tau$
0.52	78.3	20.7	6.5

. 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 their 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 is 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 [30] provides a description of the propagation of such noise-like signals 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, according to the noise wave model, is related to the averaged sky temperature $T_{\mathrm{sky}}$ by

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

(4)

where the transfer function of the noise wave reflected between the antenna and receiver is

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

(5)

and the reflection coefficients of the antenna and receiver are, respectively,

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

(6)

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

(7)

here, $Z_0$ denotes the characteristic impedance of the transmission line and can be viewed as a reference impedance, and $T_0$ represents the receiver noise offset. The noise reflected from the antenna reenters the receiver with phase $\varphi$, which is equal to the phase of ${\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, $|{\Gamma }_s|=1$ for an open circuit source ($Z=\infty$) or a shorted source ($Z=0 \Omega$), while $|{\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 $\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 perform calibration and signal reconstruction. Here, we use a 50 $\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; therefore, 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,$ and $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 improve the accuracy of the noise wave parameter determination.

We focus on how the systematic errors propagate in the global spectrum measurement system, as well as their impact on the recovered signal. We first consider the time-independent case; the time variation in the system parameters such as the gain, the noise excess temperature, and the ambient temperature are discussed in the 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 a perfect match, a fraction of the signal power would be reflected back to the antenna. Moreover, the receiver also produces noise that propagates back to the antenna. The receiver noise can be separated into a component that is correlated with the receiver output and a component that is uncorrelated [31]. 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 Equation (4).

For receiver calibration, the receiver input can also be switched to a 50 $\Omega$ ambient resistor load or an open cable. With the 50 $\Omega$ resistor load as input, the sky temperature $T_{\mathrm{sky}}$ is replaced by the ambient temperature $T_{\mathrm{amb}}$, and the receiver noise temperature can be written as

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

(8)

where

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

(9)

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

(10)

${\Gamma }_{\mathrm{rec}}$ is given by Equation (7) since the same receiver is used, $\varphi$ is the phase of ${\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\Omega$, then ${\Gamma }_{\mathrm{load}}\approx 0$. Equation (8) is reduced to

$$\begin{array}{c}\hfill T_{\mathrm{load}}=T_{\mathrm{amb}}(1-|{\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-|{\Gamma }_{\mathrm{open}}{|}^2\left)\right|F_{\mathrm{open}}{|}^2+T_u|{\Gamma }_{\mathrm{open}}{|}^2{\left|F_{\mathrm{open}}\right|}^2+\hfill \\ \hfill & (T_{c}cos\left(\varphi \right)+T_{s}sin\left(\varphi \right))|{\Gamma }_{\mathrm{open}}\left|\right|F_{\mathrm{open}}|+T_0,\hfill \end{array}$$

(12)

where

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

(13)

$${\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, which is a function of frequency $\nu$. For open cable, $Z_{\mathrm{in}}=-j{Z}_{0}cot\beta l$. The cable acts 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) Calibrate for the 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}}$, 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,$ and $T_s$. (3) Using the system parameters obtained in the first two steps to recover antenna temperature from the raw data, reconstruct the sky temperature. (4) Fit 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 ([38]) software.

We use some typical and representative models and parameters for receiver and system settings to be as close as possible with the actual observation systems, though it is understood that for different brands and types, the parameter value could vary. 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 the low noise figure and flatness of $S_{11}$ and $S_{22}$. The Wantcominc WHM0003AE [39] 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 typically offers a 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 m 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 $\Omega$, and the velocity of propagation in the cable is 83% of the 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 $\Omega$.

The reflection coefficients ($S_{11}$) of the antenna, the LNA receiver, and the open cable are shown in Figure 6 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 [31]. The measurement accuracy of the reflection coefficients ($S_{11}$) of the antenna, cable, and 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, and drift and stability. A more detailed description of the source of a VNA error is given in Appendix B.

For a given set of VNA, standard calibrator, and 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 Equations (A17) and  (A24). Below, we take the ${\mathrm{Keysight}}^{\mathrm{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. Corrected system performance of the ${\mathrm{Keysight}}^{\mathrm{TM}}$ N5247B PNA-X series vector network analyzer with the 85050C standard calibration kit (50–500 MHz). Dynamic accuracy is specified by the condition of −10 dBm input power at 50 MHz.

Directivity	Source Match	Load Match
Mag	Mag	Mag
48 dB	40 dB	47 dB
Reflection Tracking		Dynamic Accuracy
Mag	Phase	Mag	Phase
$\pm 0.0030$ dB	$\pm 0.020°$	0.009 dB	$0.08°$

lists the 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 [40].

Figure 5 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. The 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, although 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. Figure 6 shows the reflection coefficients of the antenna, the LNA receiver, and the open cable with measurement errors.

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 call receiver errors, and then those induced by the VNA measurement errors, which we 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 approximately 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,$ and $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}$, and $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, and the FFT length is 16,384, to produce 8192 output spectral bins. A total of 500,000 s (138.89 h) 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 s), $\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 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 s. The data for the open cable show an oscillation with respect to the frequency, as the reflection at the end of the cable sets 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,$ and $T_s$ as a function of frequency by the 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 Equation (1) in the middle panel, and the foreground plus 21 cm signal (Equation (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 among themselves.

In the above, we have assumed that the system is time-independent, ignoring the time variation in the system characteristics. However, the gain of the system drifts over time, especially if there are temperature variations in the circuit. If the overall gain drifts ±0.5 dB, while keeping other system parameters stable and bandpass shape unchanged, the simulations show that it causes approximately 20 mK deviation in the results 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 s 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 Equation (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 10 s, as before, while we use three cycle’s thermal load data (30 s) for receiver calibration. This induces 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 a 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

3.2.1. Error Propagation in Sky Temperature

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 5, 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 5. 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, although 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 6, 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 the 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|{\Gamma }_{\mathrm{ant}}\right|}^2){\left|F_{\mathrm{ant}}\right|}^2}-\frac{T_u{\left|{\Gamma }_{\mathrm{ant}}\right|}^2}{1-{\left|{\Gamma }_{\mathrm{ant}}\right|}^2}-\frac{(T_{c}cos\left(\varphi \right)+T_{s}sin\left(\varphi \right))\left|{\Gamma }_{\mathrm{ant}}\right|}{(1-{\left|{\Gamma }_{\mathrm{ant}}\right|}^2)\left|F_{\mathrm{ant}}\right|}\hfill \end{array}$$

(15)

The errors in each quantity 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 {\Gamma }_{\mathrm{ant}}}\delta {\Gamma }_{\mathrm{ant}}+\frac{\partial T_{\mathrm{sky}}}{\partial {\Gamma }_{\mathrm{rec}}}\delta {\Gamma }_{\mathrm{rec}}+\frac{\partial T_{\mathrm{sky}}}{\partial {\Gamma }_{\mathrm{open}}}\delta {\Gamma }_{\mathrm{open}}\end{array}$$

(16)

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

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

(17)

As $T_u,T_c,$ and $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, ${\Gamma }_{\mathrm{rec}}$ is present in the sky expression directly through $F_{\mathrm{ant}}$ (Equation (5)) and indirectly by affecting $T_u,T_c,$ and $T_s$ (Equation (12)) and $T_0$ (through Equation (11)). Therefore, we have

$$\begin{array}{ccc}\hfill \frac{\partial T_{\mathrm{sky}}}{\partial {\Gamma }_{\mathrm{rec}}}& =& {\left.\frac{\partial T_{\mathrm{sky}}}{\partial {\Gamma }_{\mathrm{rec}}}\right|}_{\mathrm{direct}}+\frac{\partial T_{\mathrm{sky}}}{\partial T_u}\frac{\partial T_u}{\partial {\Gamma }_{\mathrm{rec}}}+\frac{\partial T_{\mathrm{sky}}}{\partial T_c}\frac{\partial T_c}{\partial {\Gamma }_{\mathrm{rec}}}+\frac{\partial T_{\mathrm{sky}}}{\partial T_s}\frac{\partial T_s}{\partial {\Gamma }_{\mathrm{rec}}}+\frac{\partial T_{\mathrm{sky}}}{\partial T_0}\frac{\partial T_0}{\partial {\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 ${\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°$. In the frequency range of 50–100 MHz, these errors only have slight variations.

In Figure 11, 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 varies 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 VNAs are switched to make these measurements. The results show that for an error equal to 10% (0.1) of the nominal error, there is a deviation of approximately 0.5 K in reconstructed sky temperature, while for the nominal value, the deviation is approximately 6 K in sky temperature at 63 MHz.

3.2.2. Error Propagation in 21 cm Signal

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 with 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 they 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 the 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–100 MHz band are approximately 0.001, 0.001, and 0.002 (linear), respectively, and the RMS deviation value of the 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 also affects 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 in the VNA phase measurement. The phase measurement error also depends on the reflection coefficient. As shown in Figure 5, for a reflection coefficient 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 ranges from the relatively small $0.03°$ to the maximum uncertainty of $0.3°$, while for the receiver, the uncertainty ranges from the relatively small $0.15°$ to the maximum uncertainty of $1.5°$. 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–100 MHz band are approximately $0.48°$, $0.78°$, and $0.15°$, respectively, and the RMS deviation value of the 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 the VNA measurement errors are added in, we also simulated the cases with open cables of length from 3 m to 30 m. The result for the case of a 15 m 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 also increases 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$, and $T_s$ are unchanged, just as when 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

$${\Gamma }_{\mathrm{a}}^{\prime }=S_{22}+\frac{S_{12}{S}_{21}{\Gamma }_{\mathrm{a}}}{1-S_{11}{\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}$, and $G=S_{21}^2$, then

$$\begin{array}{c}\hfill {\Gamma }_{\mathrm{a}}^{\prime }=G{\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 an 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. The 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 the 21 cm signal for an 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 a 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 increase in the receiver noise temperature results in an increase in thermal noise in the system, which needs a longer integration time to achieve the desired signal-to-noise ratio.

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

From the above analysis, we see the addition of an attenuator between the antenna, and the LNA could reduce the impact of VNA measurement uncertainty but increase the thermal noise. To obtain 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, whereas adopting a greater attenuation would only reduce the deviation slightly but would require a considerably longer integration time.

3.4. Discussion

We estimate the induced errors on the final 21 cm signal. The central frequency of the 21 cm 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 the antenna, receiver, and open cable measurements, respectively, could result in a 200 mK deviation in the 21 cm absorption depth A and a 2 MHz deviation on the absorption profile width W. The typical phase measurement errors of $0.48°$, $0.78°$, and $0.15°$ on the antenna, receiver, and open cable measurements, respectively, 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 [41,42,43,44], 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 in the 21 cm absorption depth caused by the VNA error, if not properly corrected in the data reduction, might confuse the astrophysical interpretation of the observed 21 cm 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 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 [45,46]. The VNA error model can also be refined so that a number of parameters could be individually determined by additional calibration measurements [47]. 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 [48]. Moreover, some novel VNA calibration methods have been developed, with potential possibilities that can be applied to the global spectrum measurement experiments [49].

4. Conclusions

Detecting the 21 cm global spectrum signal from the cosmic dawn and the epoch of reionization 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 focusing 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 for a reasonably long integration time. Furthermore, a stable and constant deviation in the system gain or receiver noise offset temperature only has a slight impact on the recovered signal.

In addition to the broad frequency band needed for the cosmic dawn 21 cm signal measurement, a mismatch in the impedance of the various RF components inevitably generates 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 the 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, or 0.002 (linear) existing in the antenna, the receiver, or the open cable, respectively, would result in a 200 mK deviation on the expected signal. A typical phase measurement error of $0.48°$, $0.78°$, and $0.15°$ in the antenna, the receiver, and the open cable, respectively, would cause a 40 mK deviation. Results also show that an attenuator benefits to the effect that the VNA measurement uncertainties caused. If these kind of systematic errors cannot be reduced in calibration procedures or corrected in the data analysis, they may confuse the astrophysical interpretation of the observed 21 cm signal. It is of vital importance to develop new methods or techniques to further improve the accuracy of the calibration.