Quantifying Ionospheric Effects on Global 21-cm Observations

We modelled the two major layer of Earth's ionosphere, the F-layer and the D-layer, by a simplified spatial model with temporal variance to study the chromatic ionospheric effects on global 21-cm observations. From the analyses, we found that the magnitude of the ionospheric disruptions due to ionospheric refraction and absorption can be greater than the expected global 21-cm signal, and the variation of its magnitude can differ, depending on the ionospheric conditions. Within the parameter space adopted in the model, the shape of the global 21-cm signal is distorted after propagating through the ionosphere, while its amplitude is weakened. It is observed that the ionospheric effects do not cancel out over time, and thus should be accounted for in the foreground calibration at each timestep to account for the chromaticity introduced by the ionosphere.


INTRODUCTION
A direct measurement of neutral hydrogen (HI) by measuring the 21-cm hyperfine transition line can provide a way to trace the intergalactic medium over time. However, the strong galactic foreground compared with the extremely weak signal is one of the main limitations of current 21-cm experiments, and in some cases their effects can be up to five orders of magnitude greater than the cosmological signal. In 2018, Bowman et al. (2018) reported to have found an absorption profile in the form of a flattened Gaussian centred at 78 MHz in the integrated spectrum by fitting it with the foreground model and a model for the 21-cm signal simultaneously. However, concerns have been raised for the unphysical parameters implied by their data, and the non-uniqueness of their solution (Hills et al. 2018;Singh & Subrahmanyan 2019;Bevins et al. 2020;Sims & Pober 2020).
The strong galactic foreground compared with the extremely weak signal is one of the main difficulties in achieving a detection of the 21-cm signal. Other astronomical radio sources emitting in the frequency band 30-240 MHz together are ∼ 4-5 orders of magnitude larger than the 21-cm signal in brightness. These foreground sources include galactic synchrotron radiation, thermal free-free emission from galactic dust, and emissions from extragalactic point sources. Thus, precise understanding of frequency-dependent or chromatic effects of instrumental gain, instrumental noise contribution, antenna beam shape and ionospheric effects, coupled with spatially and spectrally varying foregrounds is critical to measuring the global 21-cm signal.
These effects are particularly challenging for detecting the global 21-cm signal because its relevant band lies in the lower frequency range; galactic foreground brightness temperature, with galactic synchrotron radiation being the major constituent, increases with de-★ E-mail: yhs24@cam.ac.uk creasing frequency as a power law sync ∝ − with a spectral index of ∼ −2.54 (Bowman et al. 2018). These foregrounds are considered spectrally smooth. It seems plausible that the signal could be extracted simply by fitting a power law. However, the inevitable chromaticity introduced by the antenna negates this assumption. Moreover, the increased fractional bandwidth in the lower frequency band leads to an increased variation of antenna beams across the measurement bandwidth, resulting in larger chromatic effects. Neither chromatic antenna beam nor the sky temperature is smooth as a function of frequency; after integration, the convolution would break the assumed smoothness.
The ionosphere further complicates chromatic mixing. The ionosphere is a magnetised plasma situated at the upper atmosphere spanning from ∼ 50-600 km above the Earth's surface. This ionised layer can corrupt a propagating signal. The most significant propagation effects relevant to the sky-averaged spectra are absorption and refraction. These effects increase with decreasing frequency, which scale approximately as −2 , and are 2-3 orders of magnitude larger than the global 21-cm signal (Vedantham et al. 2014;Datta et al. 2016). This poses an especially large problem to experiments that adopt a wide field of view. In a global experiment, in order to achieve an achromatic antenna beam, it is necessary to adopt a wide field of view, as narrow beams often come with the penalty of chromatic effects. Also, a narrow-beam global experiment would be subject to cosmic variance (Muñoz & Cyr-Racine 2020). The global 21-cm signal could be measured directly by subtracting the modelled foreground contribution from the antenna temperature were our prior knowledge on the foregrounds and chromatic mixing accurate enough, without which we are left to make simple assumptions on those properties.
This work aims to address the chromatic effects due to the signal propagation through the temporally varying ionosphere, to which a satisfying solution has yet to be proposed. It analyses the chromatic ionospheric effects on global 21-cm observations by simulating the two major ionospheric layers, the F-layer and the D-layer, by a sim- plified spatial model with temporal variance. The analysis focuses on the chromatic corruptions introduced by the ionosphere after applying a chromaticity correction factor. The model incorporates actual ionospheric data during nighttime, the time of day favourable to observations for Radio Experiment for the Analysis of Cosmic Hydrogen (REACH) (de Lera Acedo 2019).

IONOSPHERIC MODEL
The model adopted in this work is based on Vedantham et al. (2014); in their paper, it is shown that chromatic refraction and absorption on incoming rays due to the ionosphere would introduce a significant distortion to the foregrounds. Likewise, the model adopted in this paper considers only ionospheric refraction and absorption. Dynamic effects such as scintillation induced by ionospheric turbulence (Crane 1977), refraction resulted from large scale traveling ion-disturbances (Bougeret 1981), etc. are assumed to cancel out over time. Effects due to magnetic fields are smaller than day to day changes in electron density and thus are also ignored. Furthermore, Faraday rotation is assumed to be unpolarised on scales comparable to the antenna beam (Vedantham et al. 2014). The ionospheric layers are assumed to be homogeneous, and have no spatial variation across the azimuth.
In this paper, we simulate the ionosphere with time-varying parameters using actual data for the F-layer (Fig. 1). The data is collected from Lowell GIRO Data Center at station Louisvale, South Africa (https://ulcar.uml.edu/DIDBase/), whose geological location is close to where REACH is situated (the Karoo, South Africa). There is currently no available data required to simulate the D-layer, so the parameters are altered based on the F-layer data. The details will be elaborated in section 5.3. Moreover, the model is simulated with an actual beam whose gain pattern is different at each frequency. We investigate the change in the detected global 21-cm signal and foregrounds when subjected to the ionosphere. Also, chromaticity correction factors are applied when analysing the data.

Refractive Index
The refractive index of the ionosphere is given by the generalised Appleton-Hartree equation (Shkarofsky 1961), in which is described by electron density, magnetic field and wave propagation path. Because the change in refraction index due to the Earth's magnetic field is insignificant (Vedantham et al. 2014), the magnetic field term is left out to simplify the complex calculation. The simplified Appleton-Hartree equation without the magnetic field term is: where p is the electron plasma frequency, and c is the electron collision frequency. The electron plasma frequency can be determined by electron density e by the equation where is the electron charge, e is the effective mass of the electron, and 0 is the permittivity of free space. The electric field of a plane wave traveling in a homogeneous ionospheric layer is given by where is the speed of light, Δ is the distance measured along the propagation path, and 0 is the initial electric field. The real part of the refractive index depends mostly on the electron density; it introduces a change in the phase velocity, resulting in refraction. On the other hand, the imaginary part depends mostly on the electron collision rate, which is determined by the electron density, gas density and temperature. Being imaginary, it serves to dampen the wave amplitude exponentially, resulting in absorption. The ionospheric model in this work is composed of two distinct parts (see Fig. 2): one is the F-layer, accounting for refraction, and the other is the D-layer, accounting for absorption. The F-layer is characterised by low atmospheric gas density and high electron density, while the D-layer by high atmospheric gas density and low electron density. Thus, the F-layer is modelled with the real and frequencydependent , and the D-layer is modelled with the imaginary and also frequencydependent , The parameters used to model the D-layer are drawn from a typical ionospheric condition (Mitra 1959;Evans & Hagfors 1968).

F-layer Refraction
Two different equations used to describe electron density distribution within the F-layer have been tested in the model: the first one assumes a homogeneous F-layer shell spanning between 200 km and 400 km from the Earth's surface. The other one assumes an inhomogeneous layer at the same height, in which the electron density e varies parabolically with height, its peak lying in the middle of the layer (Bailey 1948). The parabolic electron density as a function of height ℎ is then where 0 is the maximum electron density at mean height ℎ m , and is the semi-thickness of the layer. It, if modelled with the same maximum electron density, would, imaginably, result in a lesser degree of refraction than the homogeneous model. Due to the curvature of the Earth, the angle of the incoming ray is subject to deviation while passing through the boundaries of the F-layer, whose total deviation can be calculated by Snell's law. For the parabolic layer, the total deviation experienced by an incoming ray of frequency can be written as where is the Earth's radius. For the homogeneous model, this integral can be simplified as a discretised sum, considering only the contribution from the two surfaces on which e suffers a change. By Snell's law, any form of electron distribution can be easily calculated. Fig. 3 shows the deviation angle of incoming rays of different incidence angles as a function of frequency, calculated by equation (7), as is done in Vedantham et al. (2014).
To not further complicate the matter at hand, for the rest of this paper, the simpler homogeneous F-layer is adopted. The parameters used to model the F-layer are changed with time (every fifteen minutes) according to the data retrieved from Lowell GIRO Data Center.
The deviation angle will result in the change of direction for the incoming rays: it is shown in Vedantham et al. (2014) that a chromatic increase in sky area of a few per cent, a disruption ∼ 2-3 orders of magnitude higher than the global signal, is introduced by the F-layer refraction. To model this effect, a lensing effect is applied to the modified beam, and in so doing, the beam integrates to larger than unity over the entire sky. It also means different regions of the sky will be weighted differently according to their respective deviation angle. Some radiation close to the horizon coming from below may be able to acquire a new zenith angle smaller than 90 • after refraction, and hence could be captured by the antenna. The antenna temperature is expected to increase by ∼ 1-10 K due to refraction. The effective antenna beam due to the F-layer refraction can thus be written.

D-LAYER ABSORPTION
The model assumes a homogeneous D-layer shell spanning between 60 km and 90 km from the Earth's surface with a constant electron density = 2.5 × 10 8 m −3 . The high atmospheric gas densities at these heights yield large electron collision frequencies that can result in significant absorption. The electron collision frequency is assumed to be 10 MHz (Nicolet 1953). By equation (3), the amplitude of the electric field of a plane wave traveling in the D-layer is Since the intensity of an electromagnetic wave is proportional to the square of its amplitude, the loss factor that describes the loss of signal intensity can be written as The effective beam subject to absorption in the D-layer is then Fig. 4 shows the loss factor of incoming rays of different incidence angles as a function of frequency calculated by equation (10). While the total distance propagation path can be integrated numerically by Snell's law, it can also be approximated as where Δℎ is the width of the D-layer, and ℎ is the mean D-layer height. The absorption level in the relevant frequency range 50-200 MHz suggested by our model can be as large as 15%. Combining the effects of both the F and the D-layer, the effective beam can be written: The simulated antenna temperature A can then be calculated by the modified integral, where f is the sky temperature at a given time.

All-Sky Map
In this paper, an all-sky map of base frequency 0 = 230 MHz structured in HEAL is used. The sky map is an instance of the 2008 Global Sky Model (https://arxiv.org/abs/0802.1525). It is scaled to a different frequency by the equation where CMB ∼ 2.725 K is the average temperature of the cosmic microwave background radiation, and is the spectral index. The spectral index varies across a realistic sky. However, the chromaticity function (equation 16) can correct the beam's chromaticity perfectly only when a constant spectral index is assumed. Therefore, in order to completely remove the chromatic effects introduced by the beam in the case where there is no ionosphere, with which the ionospherically compromised cases are compared, most of the cases shown in this paper will assume a constant spectral index = 2.5, a value that lies within a theoretical range derived by Mozdzen et al. (2017). Nevertheless, it is also critical to examine the signal in a more realistic scenario, so cases with a spatially varied set of spectral indices are also tested.

Beam Chromaticity Correction
The beam adopted in this model is the chromatic log-spiral antenna beam. In order to remove the chromaticity distortion from sky observations introduced by the chromatic antenna beam ( , , ), the integrated antenna temperature (equation 14) should be divided by a correction factor. The following equation for such a factor is suggested in Mozdzen et al. (2017): where 0 denotes the reference frequency. This correction factor works under the assumption that the spectral index is constant across the entire sky. Moreover, the constant CMB in equation (15), if not removed from the sky temperature, would introduce a component that does not scale to the frequency dependent power law. Therefore, in order to remove the beam chromaticity perfectly in our model, the following equation is adopted: .
The difference between the integrated antenna temperature corrected by equation (16) and (17) exceeds 0.3 K at lower frequencies, which is significant enough to corrupt a signal as weak as the global 21-cm line. This difference is also apparently chromatic. A similar modification should apply also to the antenna temperature integration; it should be subtracted before it is corrected by the chromaticity correction, and then added back afterwards: When these assumptions hold, the integrated antenna temperature at each frequency will be readjusted towards the corresponding beam gain at the chosen reference frequency, and it will remove the chromatic effect introduced by the originally chromatic beam. In order to see how well this equation could apply to the chromatic ionospheric effects, the modified equation is also tested: where an ionospherically affected beam is adopted: In the cases where simulations span over a period of time, the chromaticity function is calculated and applied at each time step before time integration. It should be noted, however, that this simple correction is only used in this paper to demonstrate the degree of change the ionosphere can cause. In practice, the REACH experiment (de Lera Acedo 2019) adopts a model using full Bayesian pipeline technique that allows a more sophisticated parameter estimation (Anstey et al. 2020).

Data Structure and Gridding
The total antenna temperature is integrated in the ring pixel scheme adopted in the HEAL coordinate instead of the conventional spherical coordinates suggested in equation (14). In our model, there are 12 × 512 2 ( side = 512) pixels in total, corresponding to a resolution of ∼ 0.002 in radian, or ∼ 6.87 minutes of arc. High resolution is required for this model, as the F-layer deviation can be very small at small zenith angle and at higher frequencies, as shown in Fig. 3.
While the sky map is stored in the HEAL structure, the logspiral beam is designed in spherical polar coordinates, and thus needs to be transformed and remapped into the HEAL structure by interpolation before the pixel-wise multiplication of the sky map and the effective antenna beam (equation 20). Before being subjected to the ionospheric expansion imposed by the F-layer, the beam is renormalised to make sure the beam integrates to unity over the entire sky at every frequency.

Modelled Global 21-cm Signal
In the simulations presented in this work, a 21-cm signal, assumed to be a normal Gaussian: is inserted to the sky map to test the effect the ionospheric disturbances have on the existing method to extract the signal. The modelled 21-cm signal is assumed to be constant across the sky with no spatial variation. Gaussian signals motivated by Bowman et al. (2018) as well as those of various amplitudes 0 , central frequencies , and width are tested. The sky temperature with an inserted global 21-cm signal is given by where f ( , , ) is f, in equation (15). This is the sky temperature we use in the model when simulating with a global 21-cm signal.

Data Analysis
It is essential to know the change in the detected antenna temperature the ionospheric effects would introduce, for it is what affects the values we get from direct observations. As the simulation considers a ionosphere that varies in time, it should be noted that the change in antenna temperature is dependent on the sky region at which the beam is pointing, as well as the ionospheric condition at the time of observation. It is speculated that the chromatic ionospheric disturbances would cancel out over time, so, in such cases, in order to make a fair comparison between different integration periods, the integrated antenna temperature needs to be scaled to a common ground.
A similar problem arises when it comes to the chromaticity correction mentioned earlier. Equation (16) shows that when this chromaticity correction is applied, it is corrected with respect to the beam gain ( 0 , Ω) at the reference frequency 0 , where Ω is the angular position of the sky area at which the beam is pointing, which means the final product of this correction will differ according to the reference frequency; if the beam gain ( , Ω) at 1 is smaller than at 2 by a factor of at Ω, (Ω) ( 1 , Ω) = ( 2 , Ω), with ∫ ( , Ω) Ω = 1 for all frequencies: where is a constant. As shown in the derivation, the difference is not a simple constant but a factor convolved with CMB and it is also dependant on the angular location. This correction is adopted in some of the calculations just to show a fairer comparison when a different reference frequency is applied, but a fair comparison might not be needed in practice since the raw, unscaled, value is perhaps what really matters in the experiment.
Besides the difference in magnitude, it is perhaps more critical to see how the chromatic profile is altered; it is the intermixed chromatic effect that confuses the signal measurement. As the foreground can be said to obey a power law sync ∝ − , the th order log-polynomial fitting of the following form is adopted to fit the foregrounds: In this paper, we adopt a 5 th log-polynomial to fit the data, as including higher order terms would risk modelling out the cosmological signal. This technique has been proposed in the literature (Pritchard & Loeb 2008;Harker et al. 2012) and is widely used in analysing global signal. However, the merit of this kind of fitting works well under the assumption of a smooth foreground; it falls short when a chromatic beam is used, where the spatial variation in the sky temperature is coupled into the spectral structure, especially for a wide frequency band. necessarily imply a Since the simulations include a signal in the form of a Gaussian, the complete form of the model that is used to fit the simulated data is: with constraints on each parameter, the values of which depend on the signal added in the simulation. If our model for the foregrounds is accurate enough, the latter form (equation 25) should fit the simulated data better than the previous form (equation 24) that only takes the foregrounds into account. Both forms of fitting will be shown in section 5. Also, one must stress that the yielded residuals after fitting do not necessarily imply an ionospheric feature, as any distortion to a originally smooth foreground would introduce an undulation to the residual of a log-polynomial fit.

SIMULATION RESULTS
In this section, the plots and analyses of the simulated data generated by the aforementioned model will be shown. The effects coming from the two different layers will be looked into separately before moving on to the simulations where both layers are in effect.

Separate Effects
Fig . 5 shows that the F-layer alone increases the foreground temperature by ∼ 10 K at 50 MHz, resulting in a ∼ 0.2% increment compared to the original antenna temperature, while the D-layer reduces the signal by ∼ 800 K at 50 MHz, resulting in a ∼ 16% deficit compared to the original antenna temperature. The D-layer dominates the ionospheric effect in magnitude, about 15 times that of the F-layer. However, we are more interested in how it affects the shape, which, when unaffected by the ionosphere or other effects, should resemble the power law that describes the foregrounds dominated by synchrotron radiation. Fig. 6 shows the residual after being fitted using equation (24) and (25). The fitted Gaussian signal is closer to the simulated signal in terms of amplitude and the central frequency for the case where only the F-layer is simulated. Moreover, the residual of the F-layer is weaker than that of the D-layer. It suggests that the D-layer has a greater effect on the overall chromatic variance of the integrated antenna temperature. In the presence of only the F-layer, a global  (16) with reference frequency 0 = 150 MHz. It is shown in the absolute difference (with respect to the case without an ionosphere) on a log scale. The F-layer increases while the D-layer decreases the overall antenna temperature. The D-layer (absorption) dominates the ionospheric effect in magnitude, about 15 times that of the F-layer (refraction). Table 1. The table shows the root mean square (RMS) of the residuals after fitting and the amplitude | 0 |, central frequency , and width (standard deviation) of the fitted Gaussian signal. The reference (in red) is the original input signal used in the model. The larger residuals in the D-layer cases suggest that it contributes more to the ionospheric distortion on the detected signal.

Residual
Fitted Gaussian Signal signal of amplitude ∼ 0.5 K is still detectable by the simple logpolynomial fitting (equation 25). However, in other cases simulated with a different set of parameters, where the residual yielded by the fitting exceeds 0.2 K, the F-layer alone can introduce significant distortions to the foregrounds. Although the signal in this simulation has bigger amplitude, ∼ 0.5 K, in reality, it is possible that the signal is as small as 0.1 K (Cohen et al. 2017). In such cases, a 0.2 K disturbance can still overwhelm the signal. Fig. 6 also shows the residual of the data integrated over a longer time interval (five nights), and a longer integration time appears to reduce the residual in the D-layer case. However, it is most likely due to the lesser effect at other time-steps that the overall distortion averages out when integrated, and this result is therefore not sufficient enough to be taken as an evidence that the ionospheric effect cancels out over time. The discussion will be elaborated in section 5.3.

Combined Effect
To get a grasp of the range of the total ionospheric effect can have on the integrated antenna temperature within the ionospheric parameter space, the parameters including electron density e , electron collision frequency c , height, and thickness, twelve sets of randomised parameters within respective realistic ranges in a normal ionospheric weather were simulated. This result suggests that the integrated antenna temperature subjected to ionospheric distortions modelled within a given scope can have a 2500 K difference at 50 MHz and the range narrows as the frequency increases. The percentage of difference also decreases with growing frequency.
It is also critical to know how the inserted global 21-cm signal could be affected due to the ionosphere. Fig. 7 shows the change due to the ionosphere simulated using randomised parameters. The plots show that, for a signal centring at 78.1 MHz, the amplitude can be weakened by ∼ 15%, which is similar to the D-layer absorption rate shown in Fig. 5. The change in the Gaussian shape shows the chromatic distortions introduced by the ionosphere, whose effects increase with lowering frequency.

Different Integration Time Length
A case with a constant and another with a time-varying D-layer, both with a time-varying F-layer, are shown in this section. There is currently no available data required to simulate the D-layer. Nevertheless, we are still interested in investigating the behaviour of a time-varying D-layer. To achieve that, we instead change the D-layer parameters according to the actual time-varying F-layer data by scaling, so that each scale factor of the same parameter of the D-layer and the F-layer corresponds to the respective scale factor suggested in the theoretical values adopted in Vedantham et al. (2014) (the same values are also adopted in Fig. 3 and 4). We have also made sure that those values lie within the theoretical ranges. Since electron collision frequency is considered null when modelling the F-layer, a constant c = 10 MHz is adopted for the D-layer. This is only an attempt to study a varying D-layer, and may not be representative of an actual ionosphere. Fig. 8 shows the residuals of the data integrated over different lengths of time. It would be convenient if the chromatic ionospheric effects cancel over a long integration time; in such a case, the ionosphere can be ignored simply by temporal integration. However, by the way the data is analysed and fitted in this work (equation 25), there is no clear sign showing the cancellation of chromatic ionospheric effects over time for both constant and time-varying D-layer; the residual is still comparable to the global 21-cm signal in magnitude even after integration over five nights. The absence of this property suggests that an alternative approach to remove this significant distortion to the foregrounds should be proposed to accurately measure the global 21-cm signal.
One can also observe that the residuals decrease by about an order of magnitude in the case of a time-varying D-layer. It is due to a milder ionospheric condition compared to the theoretical values used to simulate a constant D-layer. It suggests that the significance of the ionospheric effects depends greatly on the current ionospheric condition.

Pattern of the Residuals
When simulated without a global 21-cm signal, the pattern of the residual after log-polynomial fitting (equation 24) is highly consistent in its shape in several different cases. This pattern is shown in Fig. 9 to be dominated by the D-layer, with the F-layer having a minor effect. Fig. 9 shows the residuals after log-polynomial fitting to the data simulated with six sets of randomised D-layer parameters and without the F-layer. The frequencies at which the local maxima and minima of the residuals are situated are consistent in all cases. The amplitude fluctuation with frequency is also consistent. The physics underlying this pattern and its mathematical composition is still unclear.
This characteristic motivates us to modify the fitting equation (equation 25) by adding a term that describes this pattern, so it becomes: where 0 is a constant and ( ) is the pattern as a function of fre- Figure 6. The plots show the residuals of the signal modelled with ionospheric effects, fitted by equation (24) and (25)   In this case, the amplitude is weakened by ∼ 15%. One can also observe a chromatic distortion in the signal, that it suffers a greater change at lower frequencies.  quency. However, the results (see Table 3) are only improved when the original fits are good in the first place; in the cases where the signals are already discernible by using equation (25), the residuals are significantly reduced by using the modified equation. This equation is applicable only to data for a single time step. When applying the same method to the temporally integrated cases, this method fails to generate a proper result. This suggests that different 0 's in each time step should be calculated separately and cannot be simplified by one 0 after the data are integrated over time. An alternative way to get round this issue is to subtract this term from the data at each time step before time integration. However, the Gaussian fitting and the foreground fitting should be done in one single process, and doing so, extracting the Gaussian signal at every time step, would therefore run Table 3. The table shows the root mean square (RMS) of the residuals after fitting and the amplitude | 0 |, central frequency , and width (standard deviation) of the fitted Gaussian signal. The reference (in red) is the original input signal used in the model. The same sets of data are fitted using two different equations, one is fitted with the additional D-layer pattern term. Some cases show improvements when the D-layer pattern is also being fitted while some yield a null fitted signal.

Residual
Fitted Gaussian Signal  counter the purpose of time integration, which is to let the chromatic ionospheric effect cancel itself out over time before signal extraction.

Modified Beam Chromaticity Correction and Spatially Varying Spectral Index
We are interested in knowing whether the modified chromaticity correction that bears a similar form to equation (17), used to correct the beam chromaticity, can also help mitigate the chromaticity imposed by the ionosphere. Technically, it is correcting for both the ionosphere and the beam, for the two are convolved for cases simulated with an ionosphere.
The difference between the data corrected by the two different Figure 10. Residual when the beam chromaticity is corrected by equation (19) (integrated over 5 night, 24 hours). Data corrected by the modified equation give a slightly lower residual than by the original equation (17).
equations is small, ∼ 0.1 K at lower frequencies. Although the magnitude of this difference is still comparable to the brightness of a global 21-cm signal, the difference does not decrease with frequency like a power law, but seems to fluctuate about 0 K. The residuals are similar in both cases, with the one corrected by the modified equation having a slightly smaller residual (Fig. 10), suggesting that the modified equation has more or less improved the result. However, more cases need to be tested in order to ascertain the effectiveness of the modified correction equation.
We have also studied the residual of a case simulated with a uniform spectral index = 2.5 and another with realistic spectral indices that vary across the sky. The map with realistic spectral indices is built by interpolating the GSM and the Haslam map (Anstey et al. 2020). As the adopted chromaticity correction (equation 17) only works under the assumption that the spectral index is uniform, it is expected that the correction would fail to effectively correct for the coupled chromatic disruption caused by the antenna beam itself and the ionospheric effects when the spectral index varies spatially.

Different Reference Frequency
We study the difference between two sets of data simulated with the same ionosphere, but whose beam/ionosphere chromaticity is corrected by equation (19) using different reference frequencies. When doing a magnitude comparison, they are also corrected by a factor calculated by equation (23) to compensate for the different reference temperature at different frequency caused by a chromatic beam.
The difference of the former case can be seen as null, which is expected, since the original correction factor derived by equation (17) should make the same correction irrespective of the reference frequency, as long as the spectral index is constant across the entire sky. The latter case does not have the same property, as the beam no longer integrates to unity but to different values at different frequencies after being stretched by the ionospheric refraction. The difference, ∼ 0.2 K at lower frequencies is significant enough to overwhelm the global 21-cm signal. The residuals in Fig. 11 shows obvious changes in the residual when a different reference frequency is adopted. This is due to the intrinsic beam chromaticity being coupled with the chromatic ionospheric effects, that it no longer integrates to unity.

Different Gaussian Signals
Apart from the signal motivated by Bowman et al. (2018), two different global 21-cm signals of lower amplitudes (∼ 0.1 K) and centred at other frequencies are also tested in the simulation. The parameters of these two signals lie within the possible range of high-redshift parameters suggested in Cohen et al. (2017). So far only two have Table 4. The table shows the root mean square (RMS) of the residuals after fitting and the amplitude | 0 |, central frequency , and width (standard deviation) of the fitted Gaussian signal. The reference (in red) is the original input signal used in the model. The residual increases with growing reference frequency. The overall fit shows an improvement as the reference frequency decreases, with the exception of width, which is shown to suffer greater reduction.

Residual
Fitted Gaussian Signal been done; more cases are left as further work. Lower amplitudes are chosen because, as shown earlier, a signal of ∼ 0.5 K could be extracted, albeit with a significant residual, we want to see if a weaker signal would be completely overwhelmed.
The results are shown in Fig. 12. The residuals of the fitted data generated in these simulations can be as large as 0.6 K at the lower frequency range, similar to the case where the signal is larger, overwhelming the amplitude of the added signals. What should be noted is that the fitted Gaussian signals shown in the plots apparently overlap with the local minima in the D-layer residual pattern mentioned previously, and it is quite probable that the strong fitted Gaussian signal is only a compensation for the deviation from the pattern at the less stable points, such as those at the extrema. The fit, therefore, should not be accepted as authentic. Moreover, the fitting used in this work is not sophisticated enough to properly extract the signal from a realistic foreground. Nevertheless, it is clear that the ionospheric effects have great potential to bury a weak global 21-cm signal of amplitude 0.6 K in an otherwise smooth foreground.

REACH Pipeline Analysis
We also tried to fit the simulated data using the current version of the REACH analysis pipeline, in which no feature to correct the ionospheric distortion has been added. In the pipeline, P C , a state-of-the-art nested sampling algorithm for high-dimensional parameter spaces using Bayesian inference (Anstey et al. 2020), is adopted. It tries to fit the data using the log likelihood given by Figure 12. The panels show the residuals of the signal modelled with ionospheric effects, fitted by equation (24) and (25)  , and dotted-magenta -the 21-cm signal added in the simulation. The residuals can be as large as 0.6 K at the lower frequency range, similar to the case where the signal is larger, overwhelming the amplitude of the added signals. The fitted Gaussian signals overlap with the local minima in the D-layer residual pattern (Fig. 9), and it is quite probable that the strong fitted Gaussian signal is only a compensation for the deviation from the pattern at the less stable points at the extrema.
where is a free parameter with a log uniform prior of 0.001-0.1 K. Fig. 13 shows the reconstructed Gaussian signal, in blue, and the residual, in red. It is apparent that the Gaussian signal cannot be recovered by the current version of the pipeline when the signal is contaminated by ionospheric effects, with a residual exceeding 4 K at lower frequencies, much higher than the intended signal detection at ∼ 0.5 K.

CONCLUSIONS
The works shown in this paper can be summarised by the following points: (i) The total chromatic ionospheric effect is dominated by the Dlayer absorption both in magnitude and chromatic distortion. The Flayer has a lower, yet significant, contribution. The residuals yielded by the simple log-polynomial fit seem to share a similar pattern, the local extrema occurring at the same frequencies in multiple cases when simulated without a signal.
(ii) Irrespective of the extreme ionospheric conditions, foreground temperature can differ within a ∼ 2500 K range, depending on the ionospheric condition. The relevant parameters include electron density, electron collision frequency, height and width in both layers.
(iii) The global 21-cm signal modelled as a Gaussian can have its amplitude weakened and its shape distorted after being subjected to ionospheric interference based on the parameters used in this analysis.
(iv) Preliminary results show that the chromaticity introduced by the ionosphere does not cancel out over time for a constant or a time-varying D-layer, with a time-varying F-layer, when using our simplified spatial model. It also suggests that the significance of the ionospheric effects greatly depend on the current ionospheric condition.
(v) When simulated with a uniform spectral index, the beam chromaticity function and its modified counterpart have a limited use to correcting the antenna temperature subjected to chromatic ionospheric interference. Moreover, the modified equation yields different results when adopting different reference frequencies. The adopted chromaticity corrections have little to no use in the cases where spatially varying spectral indices are adopted.
(vi) As precise removal of the foregrounds is critical in extracting the faint global 21-cm signal, the ionospheric model should to be included in data analysis pipelines to properly account for the chromaticity introduced by the ionosphere.
The ionospheric model adopted in this work makes strong assumptions on the complex ionosphere, excluding spatial azimuthal variations and the dynamic phenomena that might affect the result significantly. Based on our results, we conclude that better ionospheric models are essential to understanding the chromatic global 21-cm signal distortion in ground-based experiments.
A ionospheric model has not yet been included in the actual foreground model adopted in REACH. As a precise knowledge of the ionospheric condition is not currently achievable, the current ionospheric model can be improved by gaining further insights on the mathematical formulation of the ionospheric effects based on the observables as well as the simulations. With some final adjustments, it should eventually be deployed in the REACH analysis pipeline to enable a more accurate removal of the foregrounds.