FlexKnot as a generalized model of the sky-averaged 21-cm signal at z ∼ 6–30 in the presence of systematics

ABSTRACT

Global 21-cm experiments are built to study the evolution of the Universe between the cosmic dawn and the epoch of reionization. FlexKnot is a function parametrized by freely moving knots stringed together by splines. Adopting the FlexKnot function as the signal model has the potential to separate the global 21-cm signal from the foregrounds and systematics while being capable of recovering the crucial features given by theoretical predictions. In this paper, we implement the FlexKnot method by integrating twice over a function of freely moving knots interpolated linearly. The function is also constrained at the lower frequencies corresponding to the dark ages by theoretical values. The FlexKnot model is tested in the framework of the realistic data analysis pipeline of the REACH global signal experiment using simulated antenna temperature data. We demonstrate that the FlexKnot model performs better than existing signal models, e.g. the Gaussian signal model, at reconstructing the shape of the true signals present in the simulated REACH data, especially for injected signals with complex structures. The capabilities of the FlexKnot signal model is also tested by introducing various systematics and simulated global signals of different types. These tests show that four to five knots are sufficient to recover the general shape of most realistic injected signals, with or without sinusoidal systematics. We show that true signals whose absorption trough is of amplitude between 120 and 450 mK can be well recovered with systematics up to about 50 mK.

1 INTRODUCTION

The cosmic epoch between recombination and reionization is still a missing piece in modern cosmology. This period of cosmic history sees the Universe through the dark ages, the first light, and then the onset of cosmic reionization. The hyperfine transition of neutral hydrogen of wavelength λ = 21 cm in its rest frame is an observable that can trace the local properties of gas in the intergalactic medium over time. A change in the brightness temperature relative to the radio background can occur when neutral hydrogen absorbs or emits at this wavelength, the level of which varies with redshift.

Existing observatories in an attempt to detect the 21-cm hydrogen line are designed following two main approaches. First, the interferometric instruments such as the SKA (Dewdney et al. 2009), HERA (DeBoer et al. 2017), PAPER (Parsons et al. 2010), MWA (Lonsdale et al. 2009), LOFAR (van Haarlem et al. 2013), and NeuFAR (Mertens, Semelin & Koopmans 2021) are designed to detect full spatially varying power spectrum of the 21-cm hydrogen line. Secondly, conceptually simpler but harder to calibrate approach adopts a wide-beam single antenna system to detect the spatially averaged, or sky-averaged, global 21-cm signal. These experiments include EDGES (Bowman et al. 2018), SARAS (Patra et al. 2013; Singh et al. 2022), BIGHORNS (Sokolowski et al. 2015), MIST (Monsalve et al. 2023), and REACH (de Lera Acedo 2019; de Lera Acedo et al. 2022).

The hyperfine transition of neutral hydrogen can be considered in terms of spin temperature, T_S, which describes the relative populations of the two energy levels. It can be defined as

where n₁ and n₀ are the relative populations of the higher and lower energy levels, g₁ and g₀ are the degeneracies of the higher and lower energy levels, respectively, and T_* is the transition energy divided by the Boltzmann constant.

The true shape of the global 21-cm signal during the epoch of reionization (EoR) is yet to be known. An example of the simulated signal by globalemu (Bevins et al. 2021b) is shown in Fig. 1: during the dark ages, gas is sufficiently dense for collisions to couple the spin temperature of the hydrogen gas to the kinetic temperature of the gas, T_S → T_K. Due to adiabatic cooling of gas, the kinetic temperature of the gas is lower than the background radiation temperature, T_K < T_γ, which is usually assumed to be the cosmic microwave background (CMB) temperature. This difference results in an absorption against the radio background, making T_S < T_γ. The collisional transitions become negligible as gas density decreases due to continuing expansion, setting T_S → T_γ. During cosmic dawn, the first luminous objects begin to emit radiation in the Lyman band. The spin temperature is then coupled to cold gas, T_S ∼ T_K < T_γ via Wouthuysen–Field effect (Wouthuysen 1952; Field 1958), resulting in another absorption. Fluctuations as well as total intensity in the Lyman α background no longer affect the 21-cm signal after Lyman α coupling saturates. The growing population of X-ray sources start heating the adiabatically cooling gas, driving T_S ∼ T_K > T_γ in general. In some models, this would turn absorption into emission. Ionizing photons begin to turn neutral hydrogen (H i) to ionized hydrogen (H ii), and then the signal eventually tends to zero (Shaver et al. 1999; Furlanetto 2016). Knowing the location and the shape of the global 21-cm signal absorption trough can enlighten us on the timing of the primordial star formation, the mass and star formation efficiency of the first star-forming haloes, and the luminosity of the first X-ray sources (Furlanetto, Oh & Briggs 2006; McQuinn & O’Leary 2012).

The EDGES group (Bowman et al. 2018) reported to have detected an absorption profile in the shape of a flattened Gaussian centred at 78 MHz with an amplitude of 500 mK in the sky-averaged 21-cm spectrum. Re-examination of the EDGES data analysis (Hills et al. 2018; Singh & Subrahmanyan 2019), however, has pointed out the non-physical parameters, the non-uniqueness of their solution, as well as potentially unaccounted for systematic structures in the data. Studies such as Sims & Pober (2020) and Bevins et al. (2021a) have shown that a damped sinusoidal systematic is strongly preferred in the EDGES data. Further observations such as SARAS3 also reported non-detection of the EDGES profile with 95.3  per cent confidence (Singh et al. 2022).

Existing methods in data analysis and interpretation such as using the shape of a Gaussian or a flattened Gaussian as a signal model cannot describe both absorption and emission peaks at the same time. On the other hand, physically modelled signals (Mesinger, Furlanetto & Cen 2011) such as those of 21cmGEM (Cohen et al. 2020) (available at https://people.ast.cam.ac.uk/~afialkov/Publications.html), and globalemu (Bevins et al. 2021b) have limited signal trough without enhanced radio backgrounds. In this paper, we adopt and test a signal model that has the potential of describing a variety of shapes while still having the ability to separate the signal from the foregrounds as well as the systematics. This is realized by adopting a function parametrized by freely moving knots stringed together by splines called the FlexKnot model. The FlexKnot method is introduced by Vázquez et al. (2012a) to reconstruct the primordial power spectrum from the Planck data, and then adopted by Millea & Bouchet (2018) for parametrizing the reionization history. Heimersheim et al. (2022) used it to parametrize the EoR history and the potential of hypothetical high-redshift Fast Radio Bursts to constrain cosmic reionization.

Recently, Heimersheim et al. (2024) have also adopted the FlexKnot method as a global signal model. They used FlexKnot to separate the non-foreground component using EDGES low-band data. There are several differences between our implementation and the approach of Heimersheim et al. (2024). In contrast to the work we present here, they do not distinguish between the cosmological signal and systematics. They modelled the foregrounds using a polynomial, while we simulate the foregrounds with a physically motivated method. Moreover, instead of adopting Piecewise Cubic Hermite Interpolating Polynomial (pchip) to perform interpolation like in aforementioned works, our model starts by parametrizing the second derivative of the global 21-cm signal via linearly interpolated splines, and the signal is recovered by integrating it twice. We test the FlexKnot signal model by implementing it in the REACH data analysis pipeline (Anstey, de Lera Acedo & Handley 2021) where the physically motivated foreground and the signal are jointly fitted using the Bayesian nested sampling algorithm PolyChord (Handley, Hobson & Lasenby 2015). We also test the FlexKnot model on cases in which uncalibrated sinusoidal systematic structures are present. Unaccounted-for systematics would potentially arise in practice, due to e.g. ground plane artefact (Bradley et al. 2019) and calibration issues (Sims & Pober 2020), so it is crucial to know how susceptible the signal model is to the presence of systematics.

This paper is organized as follows. In Section 2, we detail how the function is parametrized and introduce the additional constraints applied to the signal model. Section 3 briefly describes how the foreground is modelled in the REACH data analysis pipeline as well as Bayesian inference on which the pipeline is based. The results are presented in Section 4, which covers comparison between different signal models, signal recovery in the presence of sinusoidal systematics in the simulated antenna temperature data, and lastly the optimal number of knots. Section 5 concludes the work.

2 CONSTRUCTING THE FLEXKNOT FUNCTION

The FlexKnot model used in this paper is based on Vázquez et al. (2012a), who proposed the model to reconstruct the primordial power spectrum from the Planck data. It has also been adopted in several other works for its higher degree of flexibility (Bridges et al. 2009; Vázquez et al. 2012b, 2013; Abazajian et al. 2014; Aslanyan et al. 2014; Hee et al. 2016; Planck Collaboration XX 2016; Hee et al. 2017; Millea & Bouchet 2018; Olamaie et al. 2018; Planck Collaboration I 2020a; Planck Collaboration X 2020b; Heimersheim et al. 2022; Escamilla & Vazquez 2023). In this paper, it is formed by a number of knots interpolated by splines. It has the advantage of describing functions of all possible shapes, given sufficient number of knots; this attribute is desirable because the exact shape of the global 21-cm signal is yet to be known. Normally, there is a deep absorption trough at z ∼ 10–30 and an emission at z ∼ 6–10; but this is model-dependent. To model the global 21-cm signal, piecewise cubic spline interpolation seems to be the preferable choice. In this work, the FlexKnot function is constructed by integrating twice over a function of freely moving knots interpolated by piecewise linear splines:

which is the equivalent to a function made up of piecewise cubic splines, and works around the difficulties in interpreting the posteriors and setting priors when directly employing cubic spline interpolation with nested sampling (Handley et al. 2019). The positions of the knots are determined by Bayesian nested sampling. Fig. 2 shows an example of a signal recovered by the two different methods: the upper three panels show how the FlexKnot function is implemented as a signal model in three steps, and the last panel shows the direct cubic spline interpolation done via pchip. To recover a signal alike to a Gaussian, the number of knots required in the two different methods is different. The reconstructed Gaussian-like signal in Fig. 2 requires two more knots for the second derivative parametrization method to recover a Gaussian signal with reasonable accuracy than the direct cubic spline interpolation method.

The FlexKnot signal model is then implemented in the REACH data analysis pipeline to test its functionality. The REACH data analysis pipeline, based on Bayesian nested sampling, is detailed in Section 3. In the cases of only a few knots, the FlexKnot function modelled via parametrizing the second derivative is suitable for describing shapes with obvious troughs or crests, such as the shape of the global 21-cm signal, which, typically, is predicted to have one bigger trough and an ensuing emission in the range z ∼ 6–30. To describe a smooth shape like the power-law function accurately, it would require a much higher number of knots distributed throughout the argument range, as the second derivative of a power law is non-linear except when the order is 2 or 3. This characteristic can also potentially help separate the global 21-cm signal from the smooth foregrounds. Fig. 3 shows an example of how the FlexKnot model of different N_knot aims to recover a power-law function:

In this case, it requires nine knots for it to be properly recovered, and a lower number of knots leads to inaccurate results. The optimal number of knots to use in each case would also depend on the step size in numerical integration.

2.1 Parameters and priors

Each knot K_i is defined by two parameters: the position in the observing frequency and the second derivative in temperature. There is an additional highest frequency knot K_{N + 1} fixed at 0 K after which the integrated value remains zero, as it is known that reionization will eventually eliminate the 21-cm signal. This additional highest frequency knot is excluded from the total number of knots N_knot. The total number of parameters n_signal of the FlexKnot signal model is then

As the FlexKnot signal model is built in the mindset of the Bayesian nested sampling technique, priors need to be set for the parameters. A sorted prior is adopted for the positional parameter of the knots so that the positional parameters would always be in ascending order (see Appendix  A). The prior range of the positional parameters is [50, 200] MHz (sorted uniform prior), which covers the observational frequency range of the REACH experiment.

For the second derivative parameters, the prior is set to be [−0.02, 0.02] (uniform prior); the prior limits should be adjusted depending on the interval between the frequency points, as the derivative describes the rate of change between two discrete points in this model. In our tests, the interval between the frequency points is set to 1 MHz. The second derivative of the lowest frequency point (the first knot K₁), however, is set to be negative only [−0.02, 0] (uniform prior) to prevent the function from turning uncharacteristically high at the low frequencies corresponding to higher redshifts, which could potentially occur when the chromatic foregrounds in the form of power-law function are present in the data. The additional highest frequency knot K_{N + 1} does not have the second derivation parameter; its twice integrated value is fixed at 0 K, with both its first and second derivative being also null.

2.2 Theoretical dark ages primer

Unlike during the EoR, the theoretical global 21-cm signal is better understood during the dark ages (Hogan & Rees 1979; Scott & Rees 1990; Mondal & Barkana 2023). To take advantage of that, in our model (Fig. 4), a set of theoretical values at frequencies corresponding to the dark ages (z ∼ 30–40) lying outside the observational frequency range of REACH is implemented as a primer, or the first segment that initializes the ensuing function, to not only provide a reasonable and reliable constraint, but also further prevent potential uncharacteristic positive surge at the lower frequencies that sometimes occurs, as the case shown in Fig. 5. The low-frequency end of the function is then interpolated between the last point of the theoretical dark ages primer and the first knot. The value of the theoretical dark ages signal is generated by the globalemu emulator (Bevins et al. 2021b) modified for the dark ages where different parameters are used. The global 21-cm signal during the dark ages is characterized by the following astrophysical parameters: baryon density, Ω_b = 0.01495074, matter density, Ω_m = 0.29395689, curvature density, Ω_k = 0.05966342, and reduced Planck value, h = 0.75665809 km s⁻¹ Mpc⁻¹.

3 DATA ANALYSIS PIPELINE

To test the FlexKnot signal model, the REACH data analysis pipeline (Anstey et al. 2021) is adopted to extract the injected global 21-cm signal from the simulated antenna temperature data (Section 3.2). Simulated data that include the foregrounds, the injected global signal, systematics, and noise would be generated and then fed into the pipeline for it to separate the signal from other elements. The pipeline is built with the purpose to avoid degeneracy with potential systematics by adopting physically motivated foreground modelling, where the foreground model and the instrument model are jointly fitted with the signal model using the Bayesian inference technique. The Bayesian nested sampling algorithm PolyChord is opted for its ability to estimate Bayesian evidence efficiently in the case of high model dimensionalities, which is inherent to this implementation.

3.1 Bayesian inference

Bayesian inference is a statistical modelling technique that has its merit in parameter estimation and model comparison. A model |$\mathcal {M}$| parametrized by |$\theta _\mathcal {M}$| can be used to calculate the probability of observing the data |$\mathcal {D}$| by updating previous knowledge of the parameters, the prior |$\mathrm{P}(\theta _\mathcal {M}|\mathcal {M})$|⁠. This can be done by applying Bayes’ theorem,

where |$\mathcal {P}$| is the posterior distribution, |$\mathcal {L}$| is the likelihood, the probability of the data given a model and the set of parameters describing the model, |$\mathcal {\pi }$| is the prior distribution of the parameters, and |$\mathcal {Z}$| is the Bayesian evidence or marginal likelihood, which gives the probability of observing the data |$\mathcal {D}$| given the model |$\mathcal {M}$| (Sivia & Skilling 2006). One can achieve marginalization by integrating over the prior distribution:

To compare different models, one can derive the probability of a model given the data by applying Bayes’ theorem on the Bayesian evidence:

where |$\mathrm{P}(\mathcal {D})$| is a normalization factor independent of the model. As such, one may compare two competing models |$\mathcal {M}_1$| and |$\mathcal {M}_2$| by taking the ratio of the two evidences weighted by |$\mathrm{P}(\mathcal {M})$|⁠, or by taking the logarithmic Bayes factor:

under the assumption of uniform weighting |$\mathrm{P}(\mathcal {M}_1) = \mathrm{P}(\mathcal {M}_2)$|⁠. A positive |$\Delta \log (\mathcal {Z})$| indicates the preference of model |$\mathcal {M}_2$| over |$\mathcal {M}_1$| with betting odds of |$e^{\Delta \log (\mathcal {Z})}:1$|⁠. In the context of this work, we would first take the |$\Delta \log (\mathcal {Z})$| between the model with a signal:

and the one without:

to first make sure the model with a signal is indeed statistically favourable before proceeding to compare the |$\Delta \log (\mathcal {Z})$| yielded by the different signal models. |$\Delta \log (\mathcal {Z})$| is also referred to as the log evidence in this work.

3.2 Simulated antenna temperature data

The simulated antenna temperature data generated to test the capability of the FlexKnot signal model to extract the global 21-cm signal in the REACH data analysis pipeline include the following four components: the foregrounds, the sky-averaged 21-cm signal, the systematics, and the antenna temperature noise:

3.2.1 Foreground

The spatially varying spectral index map is simulated by the sky model based on instances of the 2008 Global Sky Model (de Oliveira-Costa et al. 2008) at 408 and 230 MHz:

with which the sky model can be generated:

It is then convolved with the beam pattern of a conical log spiral antenna (Dyson 1965) of beam pattern D(θ, ϕ, ν) to generate the foreground component:

3.2.2 Sky-averaged 21-cm signal

The Gaussian signal, the flattened Gaussian signal, and the signals emulated by globalemu are the three different types of signal that have been injected as the true signal in the simulated antenna temperature data in this work.

A Gaussian is generally written as

where A₂₁ is the absorption amplitude, ν₀ is the centre frequency, and σ is the standard deviation. The explored parameter space of the injected signal in the form of Gaussian is as follows:

1. amplitude: {155, 255, 455}  mK

2. centre frequency: {85, 105, 125}  MHz

3. standard deviation: {10, 15, 20}  MHz

Bowman et al. (2018) defined the flattened Gaussian function by the equation:

with

where w is the full width at half-maximum, and τ is a flattening factor. The unsmooth flat bottom of a flattened Gaussian is what distinguishes it from the two other types of signal. The parameters used to characterize the flattened Gaussian signal EDGES claimed to have detected have the following values: |$A_{21}\, = 0.52$| K, |$\nu _0\, = 78.3$| MHz, |$w\, = 20.7$| MHz, and |$\tau \, = 6.5$|⁠.

globalemu (Bevins et al. 2021b) is a sky-averaged 21-cm signal emulator, and the signal during cosmic dawn and the EoR is characterized by the following astrophysical parameters: the star formation efficiency, f_*, the minimal virial circular velocity, V_c, the X-ray efficiency, f_X, the CMB optical depth, τ, the slope of the spectral energy density, α, the low energy cut-off of the X-ray spectral energy density, ν_min, and the mean free path of ionizing photons, R_mfp.

Table 1 lists the parameters used to generate the signals adopted in this work. Some signals are emulated with enhanced radio background, where the parameter α is replaced by f_radio, for the deeper absorption troughs that is easier to recover. The models with an enhanced radio background are built to explain the EDGES result (Sharma 2018; Fialkov & Barkana 2019), the detected signal of which has an amplitude incompatible with the standard 21-cm models. For the purpose of this paper, only signals with absorption trough deeper than 15 mK are used as the true signal.

3.2.3 Systematics

In the cases where systematics are present, the simulated antenna temperature data include a sinusoid with an amplitude of {20, 50, 100 mK}. The phase of the sinusoids is shifted with respect to the centre frequency of the injected signal or the frequency at which the signal is at its lowest, so the phase difference between the signal and the sinusoid is the same across all cases:

where A_sys is the amplitude of the sinusoid, P_sys is the period, and ϕ_sys is the phase.

3.2.4 Antenna temperature noise

The antenna temperature noise is generated by the generalized normal distribution noise model where the noise at each frequency channel is randomly assigned based on the probability density distribution. It is a reasonable alternative to the radiometric noise and is easier to generate. The scale parameter of the generalized normal distribution noise model can be expressed as

where β_n is the shape parameter, σ_n is the standard deviation, and Γ denotes the gamma function. In this work, the shape parameter is set to β_n = 2, and standard deviation is σ = 25 mK unless otherwise specified.

3.3 Foreground modelling and fits

In the REACH data analysis pipeline, a foreground model, |$\mathcal {M}_\mathrm{f}$|⁠, and a signal model, |$\mathcal {M}_\mathrm{21}$|⁠, are jointly fitted; with the noise model, |$\mathcal {M}_{\mathrm{noise}}$|⁠, together they work as a single model |$\mathcal {M}$| to be analysed by the Bayesian algorithm (equation 9). The foreground modelling follows the same process as how the foreground component of the simulated data is generated (Section 3.2.1). Moreover, in order to take chromatic distortions into account when modelling the foregrounds, the sky is first divided into N regions of similar spectral indices, and in each region, a distinct uniform spectral index parameter is assigned for the purpose of scaling the base map (Anstey et al. 2021). The antenna temperature can then be modelled by convolving it with a beam pattern. Chromatic distortions due to beam chromaticity and non-uniform spectral index can be characterized by fitting this physically motivated foreground function. The foreground model, |$\mathcal {M}_\mathrm{fg}$|⁠, is then fit to the data jointly with the signal model, |$\mathcal {M}_\mathrm{21}$|⁠, via Bayesian inference with a likelihood function of

where i is the index of each frequency bin, under the assumption of a simple model of uniform uncorrelated Gaussian noise α_n (equation 19) across the frequency band. |$T^{*}_\mathrm{fg}$| and |$T^{*}_\mathrm{21}$| are the antenna temperature data given by the models |$\mathcal {M}_\mathrm{fg}$| and |$\mathcal {M}_\mathrm{21}$|⁠, respectively. The prior given to the spectral index parameters β_i is [2.45844,3.14556] (uniform prior), which is the full range of spectral indices in the map, and the Gaussian noise parameter [10⁻⁴, 10¹] K (logarithmically uniform prior). The PolyChord settings used in this paper are listed in Table 2.

4 RESULTS

4.1 Model comparison

In theory, the FlexKnot signal model should have the versatility to fit injected signals of different types of shapes, and for the cubic spline twice integrated from linear spline method, shapes that are smooth without rigid turns are favoured. Signals of different types of shapes, namely, the Gaussian signal, the flattened Gaussian signal, and the globalemu signal, are injected in the data to test the model’s capabilities and make comparison with the commonly used Gaussian signal model. Apart from |$\Delta \log (\mathcal {Z})$| (equation 8), signal root-mean-square error (RMSE), the RMSE between the injected signal, T₂₁, and the reconstructed signal, |$T^{*}_{\mathrm{21}}$|⁠, is also calculate to compare between different models:

A lower signal RMSE means a better fit.

4.1.1 Gaussian

In the tested cases where the injected signal is a Gaussian (the tested parameter space is listed in Section 3.2.2), the FlexKnot model usually does not recover the signal better than the Gaussian model itself. In most of the cases (Table 3 is the list of all tested cases), the FlexKnot model nevertheless recovers the signals with significant confidence in terms of |$\Delta \log (\mathcal {Z})$| and with the signal RMSE being only slightly higher than that of the Gaussian fit. Tested cases with an injected Gaussian signal yield very similar results, and Fig. 6 shows one of the cases where the true signal is a Gaussian centred at 85 MHz with an amplitude of 155 mK and a standard deviation of 10 MHz. The dashed lines in the figure are the corresponding reference values yielded by the Gaussian signal model, and the solid lines are that of the FlexKnot model of different N_knot. The signal RMSEs are lower than when using the Gaussian signal model, indicating better signal fits. The signal RMSE also shows that the presence of systematic decreases the quality of fit progressively with the level of its amplitude in all cases, with the case of 50 and 100 mK having significantly higher values, indicating the unreliability of the signal recovery. The impact of systematics is discussed more extensively in Section 4.2.

4.1.2 GLOBALEMU

Within the tested parameter space (listed in Table 1), the FlexKnot model performs better than the Gaussian model in the cases of injected globalemu signal (Table 4); it is able to recover various absorption histories of the globalemu signal yielded by different set of astrophysical parameters. The parameters used to generate the globalemu signals tested in the work are shown in Table 1.

One of the results is shown in Fig. 7. The high |$\Delta \log (\mathcal {Z})$| indicates the significant confidence in the results of all the explored cases. The signal RMSE shows a significant difference between the FlexKnot model and the Gaussian model, the former outperforming the latter. The signal RMSE may seem high for an injected signal at the absorption level of A ∼ 145 mK. It is due to the continuous significant emission at the higher frequencies in this particular case, which works against the FlexKnot model whose values after the highest frequency knot are set to zero. This is shown in Fig. 8. The discrepancies at the high-frequency range contribute significantly to the signal RMSE. If the higher frequency (emission) range is excluded, the signal RMSE becomes much lower.

4.1.3 Flattened Gaussian

Implementing flattened Gaussian-shaped global signals is motivated by Bowman et al. (2018). The shape of a flattened Gaussian is different from that of a Gaussian or a globalemu in its flat bottom, making the overall shape resemble a piecewise linear spline more than a piecewise cubic spline. In the cases of injected flattened Gaussian signals, the centre frequency is well recovered for both models. The sharp corners at the bottom of the flattened Gaussian already rule out the possibility for the Gaussian model to recover it accurately. On the other hand, as shown in Fig. 9, the FlexKnot model is able to accurately capture the width; it recovers the shape well up until the bottom of the signal, where the Gaussian is flattened. The sharp corners of the signal, as well as the flat bottom, however, are not. In the presence of the chromatic foregrounds, it often fails to capture the feature using the necessary four knots, but instead uses fewer knots to describe a smoother shape. When there are more than sufficient knots, the model would begin to fit structures such as the foreground chromaticity. The reconstructed signal plot at the bottom of Fig. 9, in which the yellow vertical lines indicate the mean value of each knot location yielded by the fit when N_knots = 7, shows that the locations of most knots are not too far off, but it is one knot short of the four knots necessary to form a flat bottom. This is partly due to the high noise level σ = 25 mK; the recovery of the flattened bottom could sometimes be improved by lowering the noise level in T_data. An alternative is to simply adopt linear spline interpolation instead of cubic spline in the FlexKnot signal model.

4.2 Signal recovery with systematics

In order to explore the capability of the FlexKnot signal model in the presence of systematics, simulated antenna temperature data with different levels of sinusoids (at an amplitude of {0, 20, 50, 100 mK}) and injected signals have been tested. Each case is tested using a range of N_knots, of which the lowest signal RMSE is recorded. To compare the fit of the signals at different absorption levels, the signal RMSE is divided by the amplitude of the true signal, yielding a dimensionless value. The results are shown in Fig. 10, where the cases are grouped by the amplitude, the centre frequency, and the standard deviation of the injected signal.

In all cases, as expected, greater systematics in T_data would yield higher signal RMSE, or poorer fits. The first panel is grouped by the amplitude of the true signal, showing that the stronger the signal absorption level, the less susceptible it is to sinusoidal systematics, and thus the better the signal recovery is.

The second panel, grouped by the centre frequency of the true signal, shows that when there is no systematic present, signal recovery is slightly better when the true signal is centred at a higher frequency. The foreground brightness decreases exponentially with frequency, and thus the signal presence becomes more distinctive at higher frequencies, the likely reason for the better results. It no longer holds, however, in the presence of sinusoidal systematics, which, depending on its phase, can distort the shape of the signal and cause a shift in the location of the signal trough.

There is no apparent difference between the groups by the standard deviation. It probably is due to the period of the added sinusoidal systematics, which is larger than the entirety of the explored parameter space covering only {10, 15, 20} MHz.

Overall, the results show that antenna temperature data with sinusoidal systematics of an amplitude higher than 50 mK could not be recovered with sufficient accuracy using the FlexKnot signal model.

4.3 Optimal number of knots

In this section, we investigate the optimal number of knots that could represent the global 21-cm signal. This is done by running the algorithm using a range of number of knots, N_knot = 4–10, and examine the fitted signals by log evidence and signal RMSE.

4.3.1 Data without systematics

The optimal number of knots in general is the same with or without the presence of foregrounds and noises for both Gaussian and globalemu signals. In the cases without sinusoidal systematics in T_data, |$\Delta \log (\mathcal {Z})$| mostly peaks at N_knot ∼ 4–6 and in general decreases with the increasing N_knot, meaning a higher N_knot does not necessarily contribute significant improvement to the signal recovery. This is true for both injected Gaussian and globalemu signals, and it means that N_knot ∼ 4–6 is sufficient to recover signals in the form of a Gaussian or globalemu. Fig. 11 shows the signal RMSE to N_knot of multiple different injected globalemu signals. The reconstructed signal plots as well as signal RMSE also show that more knots does not necessarily improve the signal recovery, and instead introduces unwanted structures to the signal fit. Interestingly, despite the different implementations and interpolation methods, our conclusion with regards to the optimal number of knots for data without systematics is similar to Heimersheim et al. (2024), in which it is shown that the evidence peaks at N_knot = 6 using EDGES low-band data.

On the other hand, the optimal number of knots is different with or without foregrounds and noises for a flattened Gaussian signal. When the data consist of a flattened Gaussian signal only, the total number of knots required to resolve the flattened bottom would be N_knot ∼ 7–8. When the chromatic foregrounds and noises are present, however, it is often the case that the flat bottom becomes unrecoverable, regressing the signal fit to a smoother bottom, the degree of which depends on the noise level. When the flat bottom is masked by the noise, the log evidence peaks at N_knot ∼ 5–6, similar to a Gaussian or a globalemu signal.

4.3.2 Data with systematics

When the systematics are at the level where the signal can still be recovered with significant confidence and low signal RMSE, mostly when A_sys < 50 mK, the optimal N_knots is about the same as the cases without systematics in T_data. Additional knots beyond that tend to start fitting the systematics mixed with the foreground residuals instead of improving the overall signal recovery. |$\Delta \log (\mathcal {Z})$| does not improve with the increasing N_knots either.

In the cases where there are large sinusoidal systematics (A_sys > 50 mK) in the data that the signal RMSE is high, one sometimes observes a significant increase in |$\Delta \log (\mathcal {Z})$| with increasing N_knot. This is still a result of the FlexKnot signal model starting to fit the sinusoidal systematics alongside the foreground residuals instead of the true signal; it is confidently fitting the other structures such as part of the foregrounds or the systematics other than the true signal in the given data.

5 CONCLUSIONS

In this work we have explored the capability of a physics-agnostic FlexKnot signal model implemented in the data analysis pipeline of the REACH experiments to recover the cosmological global 21-cm signal from simulated data. We examined the performance of the signal model with and without systematics, as well as for several different implementations of the injected global signal.

Our FlexKnot signal model is characterized by a cubic function integrated twice from a function of freely moving knots interpolated by piecewise linear splines, a different implementation of the FlexKnot signal model from previous works (e.g. Heimersheim et al. 2024). It is conditioned by the dark ages primer outside the REACH observing frequency range (50–200 MHz, which covers the period between cosmic dawn and reionization) at the lower frequency end. The second derivative prior of the first knot is constrained to be negative to prevent the potential surge of the signal, which would suggest an emission at the lower frequency end corresponding to higher redshifts that is against established theoretical predictions. An additional highest frequency knot is also included to set the ensuing values at higher frequencies to zero.

We implement the FlexKnot signal model in the REACH data analysis pipeline, where the foreground and the signal are jointly fitted, to test how well the signal model manages to recover the injected signal. For comparison, the same tests are performed using the Gaussian signal model. Three types of global 21-cm signal have been used as the true signal: the Gaussian signal, the flattened Gaussian signal, and physically motivated signals emulated by globalemu. The FlexKnot signal model can recover all tested signals with confidence and with reasonably low signal RMSE, and outperforms the Gaussian signal model when the true signal is not in the form of a Gaussian. Its ability to resolve narrow features that are not smooth, such as the flat bottom of a flattened Gaussian, however, has been shown to be weak.

In the presence of sinusoidal systematics as a proxy for unaccounted for systematics e.g. cable reflections, the FlexKnot signal model has the tendency to start fitting the systematics, as the number of knots increases. This can be avoided by setting fewer number of knots. For the tested signals, which all have absorption level higher than 120 mK, the FlexKnot signal model can reliably recover the signal with confidence and reasonably low signal RMSE when the sinusoid as systematic has an amplitude lower than 50 mK. The deeper the absorption trough is, the less susceptible it is to the presence of systematics. The accuracy in finding the centre frequency, or the minimum point, of the true signal is also affected by sinusoidal systematics.

The optimal number of knots to recover the absorption history of the physical global signal is N_knot ∼ 4–6. When there are more than sufficient knots, the FlexKnot signal model would begin to fit structures that are unwanted such as the foreground chromaticity, yielding poorer fits. Heimersheim et al. (2024) arrived at a similar conclusion, showing that the evidence peaks at N_knot = 6 using EDGES low-band data despite the different implementations and interpolation methods. It is about the same for a signal in the form of a flattened Gaussian, as the FlexKnot signal model often struggles to fit the flat bottom and eventually recovers a signal alike to the other two types of signal that are smooth throughout. The optimal number of knots remains the same even in the presence of sinusoidal systematics, and adding more knots would result in it starting to fit other unwanted structures such as the systematics and the foregrounds.

Overall, we have shown that the FlexKnot signal model could be a better alternative to the existing models such as a Gaussian fit, having the potential to recover the absorption history of various different types that are smooth with reasonable confidence and accuracy in the likelihood analysis; it is capable of separating the signal as well as recovering the features given by the theoretical prediction. It also performs well in the framework of the REACH data analysis pipeline, where the foreground and the signal are jointly fitted.

ACKNOWLEDGEMENTS

We thank Will Handley and Stefan Heimersheim for helpful discussions. ES is supported by Cambridge Trust and Taiwan Ministry of Education. DA is supported by STFC and EdLA is supported by STFC Ernest Rutherford Fellowship. AF is supported by the Royal Society University Research Fellowship.

DATA AVAILABILITY

The data underlying this article will be shared on reasonable request to the corresponding author.

References

APPENDIX A: SORTED PRIOR OF THE POSITIONAL PARAMETERS

Sorted prior is applied for the positional parameters in our model, so that the higher order knot will always be at a higher frequency than the lower one. The two corner plots in Fig. A1 show the prior and posterior distributions of a case with and another without systematics. The prior distribution in scattered blue dots lying only in the upper triangle is a characteristic of sorted prior. The overall narrower 2σ posterior distributions in the case without systematics indicates that the position parameters are being constrained.