21cm Epoch of Reionisation Power Spectrum with Closure Phase using the Murchison Widefield Array
Himanshu Tiwari Nithyanandan Thyagarajan Cathryn M. Trott Benjamin McKinley


§ ABSTRACT 
The radio interferometric closure phases can be a valuable tool for studying cosmological HI from the early Universe. Closure phases have the advantage of being immune to element-based gains and associated calibration errors. Thus, calibration and errors therein, which are often sources of systematics limiting standard visibility-based approaches, can be avoided altogether in closure phase analysis. In this work, we present the first results of the closure phase power spectrum of HI 21-cm fluctuations using the Murchison Widefield Array (MWA), with ∼ 12 hours of MWA-phase II observations centered around redshift, z≈ 6.79, during the Epoch of Reionisation. On analysing three redundant classes of baselines – 14 m, 24 m, and 28 m equilateral triads, our estimates of the 2σ (95% confidence interval) 21-cm power spectra are  ≲ (184)^2 pseudo  mK^2 at k_|| = 0.36 pseudo h  Mpc^-1 in the EoR1 field for the 14 m baseline triads, and ≲ (188)^2 pseudo  mK^2 at k_|| = 0.18 pseudo h  Mpc^-1 in the EoR0 field for the 24 m baseline triads. The “pseudo” units denote that the length scale and brightness temperature should be interpreted as close approximations. Our best estimates are still 3-4 orders high compared to the fiducial 21-cm power spectrum; however, our approach provides promising estimates of the power spectra even with a small amount of data. These data-limited estimates can be further improved if more datasets are included into the analysis. The evidence for excess noise has a possible origin in baseline-dependent systematics in the MWA data that will require careful baseline-based strategies to mitigate, even in standard visibility-based approaches.
 



§ INTRODUCTION
 

Epoch of Reionisation (EoR) is the period when the first stars and galaxies were formed in the early Universe between 15<z<5.3 and contributed to reionising the predominantly neutral intergalactic medium on cosmic scales <cit.>. It is also one of the least understood periods in the history of the Universe, mainly due to the lack of radiation influx from the first stars and galaxies, which are locally absorbed by the intervening medium. The hyperfine ground state of the atomic Hydrogen (HI) produces a weak transition of ∼ 1420 MHz, popularly known as the 21-cm line. It is considered a very promising probe of the EoR due to the abundance of Hydrogen in the early Universe. The intervening medium is largely transparent to the redshifted 21-cm line; therefore, it provides one of the best avenues to infer the astrophysical properties of the IGM and the cosmology of the early Universe. As the neutral IGM gets ionised, it weakens the strength of the 21-cm signal. 
One can interpret the stages of cosmic reionisation by estimating the depletion in the redshifted 21-cm signal through cosmic time <cit.>.To detect this forbidden transition from the early Universe, several radio instruments such as Murchison Widefield Array (MWA) <cit.>,  Hydrogen Epoch of Reionization Array (HERA) <cit.>, LOw-Frequency ARray (LoFAR) <cit.>, Giant metrewave Radio Telescope (GMRT) <cit.>, Precision Array for Probing the Epoch of Reionization (PAPER) <cit.>, Long Wavelength Array (LWA) <cit.>, Experiment to Detect the Global EoR Signature (EDGES) <cit.>, Shaped Antenna measurement of the background RAdio Spectrum (SARAS) <cit.>, Broad-band Instrument for Global HydrOgen ReioNization  Signal (BIGHORNS) <cit.>, Large Aperture Experiment to Detect the Dark Age (LEDA) <cit.>, Dark Ages Radio Explorer (DARE) <cit.>,  Sonda Cosmológica de las Islas para la Detección de
Hidrógeno Neutro (SCI-HI) <cit.>, Probing Radio Intensity at High-Z from Marion (PRIZM) <cit.> were built or are under construction. These instruments can either aim to detect the sky-averaged 21-cm signal spectrum (Global signal) or measure its spatial fluctuations. The former category of instruments can detect the overall IGM properties, whereas the latter can provide a detailed study of the three-dimensional topology of the EoR regime. 
The spatial signatures are quantified through statistical measures such as the power spectrum, which can probe the 21-cm signal strength as a function of cosmological length scales (k-modes). Alongside, the three-dimensional topology of the EoR can be studied via a two-dimensional 21-cm power spectrum <cit.>, which shows the variation of the 21-cm power spectrum along the line of sight and transverse axis.The significant challenges of detecting the 21-cm signal come from the foregrounds, ionospheric abnormalities, instrumental systematics, and Radio Frequency Interference (RFI), emitting in the same frequency range as the redshifted 21-cm signal from the EoR. In an ideal situation avoiding or minimising all the above factors, we still require to calibrate the instrument against the bright foregrounds that require calibration accuracy of ≳ 10^5 : 1 by the radio instruments to reach the required HI levels. 
Calibration accuracy is especially important for EoR observations using the MWA because of the presence of sharp periodic features in the bandpass produced by the polyphase filter bank used. The inability to accurately correct for these element-based bandpass structures significantly affects the power spectrum estimates <cit.>. The radio-interferometric closure phase has emerged as an alternative and independent approach to studying the EoR while addressing the calibration challenges. The main advantage of this approach is the immunity of closure phases to the errors associated with the direction-independent, antenna-based gains. Thus, calibration is not essential in this approach <cit.>. The closure phase in the context of the EoR was first investigated by <cit.> and further employed on the HERA data by <cit.>. The method has shown significant promise in avoiding serious calibration challenges, and with detailed forward modelling, one can ideally quantify the 21-cm power spectrum. This paper is the first attempt to utilise a closure phase approach on MWA-phase II observations. We followed the methods investigated by <cit.> and applied them to our datasets. This paper is organised as follows. In <ref>, we discuss the background of the closure phase. <ref> and <ref> of this paper explain the observations and forward modeling with simulations of the foregrounds, HI, and noise. In <ref>, we discuss the data processing and rectification, and finally, we present our results in <ref> and discuss them in <ref>.



§ BACKGROUND
 

In this section, we review the background of the interferometric closure phase in brief <cit.>. The measured visibility between two antenna factors at a given baseline (V_ij^ m) can be defined as the sum of true sky visibility and noise:

    V_ij^ m(ν) = 𝐠_i(ν) V_ij^ T(ν)𝐠_j^*(ν)  + V_ij^ N(ν);

where 𝐠_i(ν),  𝐠_j(ν) denote the element-based gain terms, {*} represents the complex conjugate, V_ij^ T(ν) is the true sky visibility, and V_ij^ N(ν) is the noise in the measurement.  The indices {ij} correspond to the antennae {a,b} forming a baseline. The true sky visibility can be further decomposed into the foregrounds and faint cosmological HI visibilities. 

    V_ij^ T(ν) = V_ij^ FG(ν) + V_ij^ HI (ν)

In general, the foreground ≥ 10^4 K orders of higher magnitude than the HI, and to reach the sensitivity limit of HI, the gains (𝐠_i's) are required to be precisely calibrated up to the HI levels. It presents challenges to the direct visibility-based HI power spectrum analysis as it requires accurate modelling of the foreground and mastering the calibration techniques.
In radio interferometry, the term `closure phase' is assigned to the phase derived from the product of N ≥ 3 closed loops of antenna visibilities <cit.>. When N=3, it is also sometimes referred to as the bispectrum phase in the literature, which can be defined as,


    ϕ_∇^ m(ν)    =  arg∏_ij=1^3 V_ij^ m(ν)  
       =   arg∏_ij=1^3 [ 𝐠_i(ν) V_ij^ T(ν)𝐠_j^*(ν)  + V_ij^ N(ν) ]  
       =  arg∏_ij=1^3 V_ij^ T(ν) + noise-like terms

where, {ij} runs through antenna-pairs {ab, bc, ca}, the gains of individual antenna elements (𝐠_i' s) gets eliminated in the closure phase; leaving only the true sky phase as the sole contributing factor in the closure phase. The closure phase delay spectrum technique also exploits the fact that the foregrounds predominantly obey a smooth spectral behavior, whereas the cosmological HI creates spectral fluctuations. Thus, in the Fourier delay domain, the foreground signal strength (power) gets restricted within the lower delay modes. In contrast, the HI power can be observed at the higher delay modes, creating the distinction between these two components in the Fourier domain, from which the faint HI can be detected. Because of its advantages in avoiding element-based calibration errors and processing simplicity, it promises to be an independent alternative to estimating the 21-cm power spectrum <cit.>. The delay spectrum of the closure phase can be estimated by taking the Fourier transform of the complex exponent of the closure phase with a window function along frequency,

    Ψ̃_∇(τ) = V_ eff∫ e^iϕ_∇(ν) W(ν) e^2π i ντ dν  ,

where τ represents delay, the Fourier dual of the sampling frequencies, and V_ eff is the effective visibility, which can be obtained through the model foreground visibilities or a calibrated visibility. In our work, we used the former estimated through foreground simulations, which are discussed in the next section. Note that we only take the amplitude of the V_ eff, which acts as a scaling factor in the delay spectrum. W(ν) is the spectral window function; we used a Blackman-Harris window <cit.> modified to obtain a dynamic range required to sufficiently suppress foreground contamination in the EoR window <cit.>,

    W(ν) = W_ BH(ν) ∗ W_ BH(ν)   ,
 
where, W_ BH(ν) is the Blackman-Harris window function, and {∗} represents the convolution operation. For a given triad, V_ eff is estimated from the sum of inverse variance visibilities weighted over the window,

    V_ eff^-2 = ∑_ij=1^3 ∫ V_ij(ν)^-2 W(ν) dν/∫ W(ν) dν  .

The inverse squaring ensures the appropriate normalisation of visibilities, V_ij, taking noise into account. From the delay spectrum, we estimate the delay cross power spectrum by taking the product of the two delay spectra and converting it into [pseudo mK^2 h^-3Mpc^3] units by assimilating the cosmological and antenna-related factors as,

    P_∇(k_||) = 2ℛ{Ψ̃_∇(τ)*Ψ̃_∇^'(τ)}×
    (1/Ω B_ eff)(D^2Δ D/B_ eff)
               (λ^2/2k_B)^2  
               [pseudo  mK^2h^-3 Mpc^3]   ,

where Ω is the antenna beam squared solid angle <cit.>, 
B_ eff is the effective bandwidth of the observation, D and Δ D are the cosmological comoving distance and comoving depth corresponding to the central frequency and the bandwidth, respectively. k_|| is the wavenumber along the line-of-sight <cit.>,

    k_|| = 2πτ B_eff/Δ D≈2πτν_r H_0 E(z)/c(1+z)^2  ,

where, ν_r is the redshifted 21-cm frequency, H_0 and E(z) are standard terms in cosmology. 

Ω B_ eff is related to the cosmological volume probed by the instrument and is defined as <cit.>

    Ω B_ eff   = ∬ |A(𝐬̂,ν)|^2   |W(ν)|^2  d^2 𝐬̂ dν  ,

where A(𝐬̂,ν) is the frequency-dependent, directional power pattern of the antenna pair towards the direction, 𝐬̂, and W(ν) is the spectral window function. However, we approximated by separating the integrals to obtain Ω as <cit.>

    Ω   = ∫ |A(𝐬̂,ν_r)|^2  d^2 𝐬̂

and effective bandwidth, B_eff, as <cit.>

    B_eff   = ϵ B = ∫_-B/2^+B/2 |W(ν)|^2 dν  ,

where, ϵ is the spectral window function's efficiency, and B=30.72 MHz is MWA's instantaneous bandwidth. For the MWA observations at the chosen band in this study, Ω≈ 0.076 Sr. For the modified Blackman-Harris window function adopted here, ϵ≈ 0.42 and hence, B_eff≈ 12.90 MHz. 

The “pseudo” in Equation (<ref>) is used to note that the power spectrum estimated via the closure phase method is an approximate representation of the visibility-based power spectrum <cit.>. Further, we used the power spectrum as defined in <cit.> with the scaling factor 2 instead of 2/3[The definition of V_ eff^2 already accounts for the factor of 3] to correct for the effective visibility estimates.



§ OBSERVATIONS
 

In this work, we used 493 zenith-pointed MWA-phase II compact observations from September 2016 under the MWA project . Each observation lasts 112 seconds in the MWA high-band frequency range of 167-197 MHz. Amongst these observations, 198 and 295 target the EoR0 (RA, Dec = 0h, -27deg) and EoR1 (RA, Dec = 4h, -30deg) fields, respectively. Figure <ref> shows the Local Sidereal Time (LST) and Julian Date (JD) of the observations. It amount to ≈ 6 and 9 hours of total observing time on the EoR0 and EoR1 fields, respectively. The MWA raw visibilities are stored in the  format. To convert them into measurement sets (MS) or uvfits, we used [<https://github.com/MWATelescope/Birli>] (an MWA-specific software that can perform data conversion, averaging in frequency and time, flagging, and other preprocessing steps). Using , we averaged the raw visibilities for 8 seconds at a frequency resolution of 40 kHz. 
Finally, we output the raw (uncalibrated) visibilities as standard . Note that we are required to keep all frequency channels for the closure phase analysis; thus, we avoided flagging channel-based RFI (e.g., DTV) and coarse band edge channels (around every 1.28 MHz), which are usually affected by the bandpass[The MWA’s signal processing chain contains filterbanks that yield 24 coarse channels of 1.28 MHz over the full 30.72 MHz band. The fine polyphase filterbank shape results in poor bandpass characteristics at the coarse channel edges <cit.>.]. 





§ SIMULATIONS
 

The closure phases are not linear in the visibilities; thus, forward modelling is key to understanding the data. We incorporated simulations of the foregrounds (FG), HI, and antenna noise to provide cross-validation and comparison with the data. Forward modelling can help identify the excess noise and systematic biases in the data and provide an idealistic estimate for comparison.





 §.§ foregrounds
 

The simulated FG are generated with the same parameters (i.e. matching LST, frequency and time resolution) as the observing data. In the first step, we generated the sky maps corresponding to each individual observation. We used the top 20,000 brightest radio sources from the PUMA catalogue <cit.> in the observing fields (EoR0, EoR1).
The catalogue includes point sources, Gaussians and shapelets. Note that our sky model does not account for the diffuse sky emission. Then, we generated the foreground sky visibilities and converted them to MWA-style . Initially, we experimented with various source counts (e.g., 15000, 25000, 45000) and their effect on the closure phase. We found the variation in closure phases saturates beyond 20000. Therefore, we settled for 20000 source counts in favour of faster computation. However, it should be noted that pinpointing the exact number of source counts where the closure phase saturates is challenging to find. The entire task of sky-map generation and foreground visibility estimation was accomplished using [<https://github.com/MWATelescope/mwa_hyperdrive>]. The sky visibilities are generated using fully embedded antenna element (FEE) beam with real-MWA observing scenario where the information of dead dipoles (if present during the observation) and antenna gains are incorporated in the simulation. Figure <ref> shows 20,000 sources around the EoR1 field for a given observation. For simplicity, only the single component of the sources (or point sources flux density)  is shown in the figure.




 §.§ Neutral hydrogen


Next, we estimate the HI visibilities as observed by MWA. In the limits of the cosmic and sample variance, the characteristic fluctuations in the HI-signal can be assumed to be the same across the sky; therefore, we can avoid simulating HI box multiple times; instead, we can use a single HI simulation box. The HI simulation was generated using 21cmFAST <cit.> with a simulation box size of 1.5 cGpc corresponding to 50^∘×  50^∘ in the sky at a redshift of 6. Then, we passed the simulated voxel data cube to [<https://github.com/JLBLine/WODEN>] <cit.> to generate the MWA-style visibilities of the HI and output it as . The HI visibilities were first generated for each 1.28 MHz coarse band separately with the matching frequency and time resolution of 40 kHz and 8 seconds of the foreground simulations and then manually stitched together to get the total 30.72 MHz bandwidth. The final HI visibilities were passed to the processing pipeline for further analysis.
The foreground and HI visibilities are added together. We computed the closure phase spectra of the foregrounds as well as of the HI imprinted on the foregrounds. Figure <ref> shows the smooth foreground spectra in the closure phase and the fluctuations (∼ 0.01 milliradian) introduced by the presence of HI. 




 §.§ Noise
 

The total noise consists of sky noise and receiver noise components. 


The receiver temperature for the MWA was assumed to be T_ rx=50 K <cit.>, while sky temperature follows a power law in our observing frequency range <cit.>,

    T_ sky = T_0(ν/ν_0)^α; T_0=180  K, ν_0=180  MHz, α=-2.5   ,


    T_ sys = T_ sky + T_ rx  .

From the system temperature, we estimated the system-equivalent flux density (SEFD) using <cit.>,

    SEFD = 2k_B T_ sys/A_ eff

and the RMS,

    σ (ν) =  SEFD/√(ΔνΔ t)  ,

where, k_B is Boltzmann's constant, A_ eff is the effective collecting area of the telescope, and Δν, Δ t are the frequency resolution and integration time, respectively.  The σ(ν) is used to generate the Gaussian random noise and converted into the complex noise visibilities with a normalisation factor of 1/√(2) in the real and imaginary parts. Finally, the noise visibilities were added to the corresponding foreground and HI visibilities to get the Model of the sky signal.











 §.§ Baseline-dependent gains

The eq. <ref> modifies to V_ij^ mod(ν) =  V_ij^ m(ν)𝐠_ij(ν); in scenarios incorporating baseline-dependent gains, where 𝐠_ij(ν) represents the baseline-dependent gain factor. We introduced the baseline-dependent gains using a simple uniform distribution in the 𝐠_ij(ν) phase with unity amplitude. The scaling factor introduced in the 𝐠_ij(ν) sampling is set to approximately match the RMS phase of the binned averaged closure phase of the DATA. We chose brute force method to find the optimal scaling factor for a given EoR field, with a the single scaling factor for given EoR field. Figure <ref> shows the binned averaged closure phase of DATA and Model with 𝐠_ij. The contribution due to the baseline dependent gains on the binned averaged closure phase about 0.05 rad in both EoR0, EoR1 fields, which is the RMS of the ratio between the Model with 𝐠_ij and Model closure phases.
From here onwards, we use two variants of models in the analysis, the first is a forward Model without baseline dependent gains and the second is a Model with 𝐠_ij.





§ DATA PROCESSING
 
The following section provides the basic data processing steps we incorporated into this analysis. The complete schematic flowchart of the data structure is shown in <ref>.




 §.§ DATA binning
 
In the first step, the repeated night-to-night observations are sorted based on the Local Sidereal Time (LST) and Julian Date (JD); see figure <ref>. We determined the 14 m, 24 m, and 28 m redundant baseline triads from both Hexagonal configurations of MWA (see an example in figure <ref> (right panel) for which the visibilities are estimated). A given triad {a, b, c} includes N_ vis = 3 visibilities which correspond to {V_ ab, V_ bc, V_ ca}.  The number of triads (N_ triads) varies depending on the baseline length. In our case, the N_ triads are 47, 32, 29 for 14 m, 24 m, and 28 m baselines, respectively. Please note that, when accurately measured, the 24 m baseline is 14√(()3)≈24.25 m; however, we chose the former for the simple denomination. On the other hand, the 14 m and 28 m baselines are nearly accurate for the antenna positional tolerance of MWA tiles. Each dual-polarisation observation was made for 112 seconds, which included N_ timestamps=14 each with 8 seconds of averaged data and a frequency resolution of 40 kHz, which provides a total of N_ channels=768 frequency channels with a bandwidth of 30.72 MHz. The entire observations can be restructured into;

    N_ obs≡{N_ LST , N_ JD, N_ timestamps, 
     
       N_ pol, N_ triads, N_ vis, N_ channels}




 §.§ RFI flagging
 

The MWA high band (167-197 MHz) lies in the digital television (DTV) broadcasting band; thus, we expect RFI to be present in our dataset  <cit.>, which in some cases can completely dominate the useful data from the observations. As mentioned in the previous sections, since our analysis required keeping all the frequency channels from our datasets, we did not perform any frequency channel-based RFI flagging in the data preprocessing step. Instead, we incorporated  <cit.>, which is designed to detect faint RFI in the MWA data, to either discard the entire frequency band or keep it based on the RFI occupancy of the dataset. Instead of assuming a persistent RFI along a frequency channel, we check for the RFI along the observation time (i.e., along the N_ timestamps axis).

The flagging was performed based on the RFI z-score of an observation. Note that the z-score was estimated at successive adjacent timestamps to measure if any faint or persistent RFI was present in the data across all timestamps (see figure  <ref>). We took a z-score threshold of 2.5, below which the data was considered good, and the channels where the z-score exceeded the threshold were considered RFI-affected. Then, we independently estimated the level of such RFI-affected channels along the frequencies at each timestamp and checked if the RFI occupancy at a given timestamp was more than 5%.  
As a first step in selecting good timestamps, we chose an RFI occupancy level of 5% as a threshold. We discarded the entire timestamp if the RFI occupancy exceeded this threshold. Figure <ref> presents RFI occupancy for the entire dataset. Since  z-scores are estimated relative to the successive adjacent timestamps, it might be difficult to quantify whether the RFI leakage between the adjacent timestamps (where one is good and another is bad) is there or not. Therefore, in the second step, we again flagged all the timestamps based on whether they had a bad neighboring timestamp. The flagging provides a masked array, further propagated through other data processing steps.


 §.§ Triad Filtering


The presence of faulty tiles or dipoles/antennae corrupts the voltages recorded by the correlator; therefore, we are required to cross-check the visibilities at each antenna triad. The easiest way was to perform a geometric median-based rectification on the closure phase. We performed a two-step median rectification on the closure phase. 
For a given observation, the data structure of the triad filtering can be shown as follows:

    ϕ_∇≡{N_ pol, N_ triads, N_ timestamps, N_ channels}

First, we estimated the median absolute deviation (MAD = Median (|X_i - X̅|)) of the closure phases against the mean along the N_ triads,

    ϕ_∇^ MAD≡{N_ pol, N_ triads, N_ timestamps, N_ channels}

and then we estimated the mean of the MAD (i.e. MAD_ mean) along the N_ channels axis. Finally, we estimated the MAD of the  MAD_ mean.

    μ{ϕ_∇^ MAD}≡{N_ pol, N_ triads, N_ timestamps, 1}

This step provides a single value of the MAD_ mean for a given triad at every timestamp. 

    MAD{ μ{ϕ_∇^ MAD}}≡{N_ pol, N_ triads, N_ timestamps, 1}

Finally, we masked the triads if the mean of the MAD is greater than 1σ≈ 1.4826 of MAD and considered the triad performing poorly at that given timestamp.





 §.§ Coherent Time Averaging
 
The coherent averaging gives us an estimate of a timescale up to which the sky signal can be assumed identical and averaged coherently to improve sensitivity. It can be estimated by measuring the variation in the sky signal with time for a fixed pointing. Indeed, these vary with instrument and frequency of observation since the beam sizes are different. To check this with MWA, we used a continuous drifted sky simulation (FG and HI) under ideal observing settings (i.e., unity antenna gains, equal antenna element elevation from the ground) for ≈ 0.5 hours while keeping a fixed zenith pointing. The sky moves about ≈ 7.5^∘ in 0.5 hours, which is less than the MWA beam size of ≈ 9^∘-7.5^∘ at the shortest (14 m) triad, thus justifying the simulation time range of 0.5 hour. 



We simply added the ideally simulated visibilities (FG and HI) and estimated the closure phase power spectrum as a function of time for higher delay (|τ|> 2 μ s). We used a fractional signal loss of 2% to measure the coherence threshold, a similar approach used by <cit.>.
The fractional loss in power is defined as,

    1-η = <|ψ_∇(t, τ)^2|> - |<ψ_∇(t, τ)>|^2/<|ψ_∇(t, τ)^2|>

The choice of |τ|> 2 μ s is to choose the timescale based on the loss of HI signal, which is where it would be significant, namely, the higher delay modes. Figure <ref> shows the fractional loss of coherent HI signal power as a function of averaging timescale for three triad configurations. We found the coherence averaging time for {∇: 14, 24, 28} triads to be approximately 408, 130, and 120 seconds, respectively. 





§ RESULTS
 

The bandpass structure of the MWA leaves significant systematics in the edge channels. Usually, these edge channels get flagged in the early data preprocessing step. However, we did not apply the flags since we require the full observing band for the delay spectrum analysis. The closure phase must be free of phase errors associated with the individual antenna elements. Thus, we expected the element-based bandpass structure to be removed in the closure phase. However, we still observed some residual bandpass structures in the closure phases; see figure <ref>. It can happen if there are some baseline-dependent gains present in the data in addition to antenna-dependent gains. We validated our hypothesis by developing a simple bandpass simulator, where we introduced an additional bandpass structure to MWA data that consisted of both element-based and baseline-dependent terms. The bandpass persists if baseline-dependent gains are present in the data, which otherwise gets completely removed if only antenna-based gains are present in the data (see appendix figure  <ref>). In total, about 128 frequency channels are affected by the bandpass, which accounts for about 16% of the total bandwidth. 




 §.§ Mitigation of baseline-dependent bandpass effects


We investigated two approaches to overcome the presence of baseline-dependent edge channel effects, the first being the Non-uniform Fast Fourier transform (NFFT), where we tried to avoid the bandpass-affected channels in the Fourier transform. The second Gaussian Process Regression (GPR) based data-inpainting to estimate the missing channel information at the location of the spikes in the bandpass spectrum.



  §.§.§ Gaussian Process Regression (GPR)
 

GPR <cit.> is a popular supervised machine learning method. GPR has previously been used in foreground subtraction <cit.>, and data inpainting <cit.> and characterisation <cit.> for 21-cm cosmology studies.
We follow a similar formalism to <cit.>. In the first step, we identified all the bad edge channels in the bandpass and flagged them in the closure phase. As mentioned before, edge channel contamination is present in the MWA data. The data at 40 kHz resolution has 768 frequency channels and 128 contaminated edge channels. We implemented GPR on the complex exponent of the closure phase exp(iϕ_∇ (ν)), with the real and imaginary components separately. The GPR implementation was rather simplistic since the complex phase varies between [0, 1]. We used the Matérn kernel as model covariance in our analysis. To optimise the kernel hyperparameters, we used the -based  <cit.> to find the minima of the objective function. The optimisation was reiterated over ten times to ensure the kernel hyperparameters' convergence. 
Note that for a given frequency range, we applied GPR to the entire N_ obs (see eq. <ref>) separately; therefore, the kernel Hyperparameters are also different for each closure phase frequency spectrum. 

Figure <ref> shows the closure phase of the data and GPR reconstruction. In the top panel, the data can be seen with spikes at regularly spaced edge channels of ≈ 1.28 MHz intervals. The interpolated values of the closure phase are plotted over the data. The bottom panel shows the difference in the RMS in the closure phase along the frequency axis. Note that, while estimating the relative difference in the data closure phase, we avoided the noisy edge channels, whereas, for the GPR case, we only included the relative difference near the edge channels. This will enable us to query whether the GPR closure phase has a similar variation across the frequency compared to the data. It can be seen that the RMS of the data is higher than the GPR values, which means that the GPR has performed quite well. We used the Python-based module [<https://github.com/SheffieldML/GPy>]  for the GPR implementation. 





  §.§.§ Non-uniform Fast Fourier transform (NFFT)
 

NFFT is a well-known method to get the Fourier transform of the data with missing samples <cit.>. To estimate the FFT of bandpass-affected data, first, we removed the 128 edge channels from the data and estimated the Ψ_∇(ν) (a Fourier conjugate of eq.<ref>), which is then supplied to the NFFT function to get Ψ_∇(τ). We used a Python-based  [<https://github.com/jakevdp/nfft>] package to develop the NFFT functions. 

The absolute values of the closure phase delay cross-power spectrum are shown in figure <ref>. It can be seen that the data is highly affected by the excess systematics in the power, evident by the periodic spikes. NFFT significantly reduces the bandpass systematics but does not eliminate it entirely. Since the GPR performed best between the two methods, we adopted only the GPR for the later analysis. From now on, we will be using `GPR-reconstructed DATA' as `DATA' for the entirety of the paper.





 §.§ Cross power spectrum estimation
 


We proceeded to the power spectrum estimation after eliminating a significant part of the baseline-dependent bandpass contamination using GPR inpainting. Based on the LST-JD of the observations (see figure. <ref>, we split the dataset equally along the JD axis (i.e., N_ JD=2). 
We binned the data along LST according to the coherent averaging time of MWA, which we already estimated for the redundant 14 m, 24 m, and 28 m baselines, resulting in 5, 14, and 17 LST bins, respectively.
The first level of weighted averaging (i.e., coherently averaging) was done to the bispectrum (visibility triple product) and effective foreground visibilities, (V_ eff), that lie within the same LST bin. The weights are the number of good observations left after being rectified by the RFI flagging and triad filtering within the given LST bin for a given polarisation and triad. 

In the next step, we estimated the delay spectra, Ψ_∇(τ), from the phase of the LST binned averaged closure spectra at each LST bin. It resulted in delayed spectra data structure, 

    { N_ LST, N_ JD, N_ pol, N_ triads, N_ delays}; N_ delays=N_ channels

Finally, we estimated the cross-power between the unique triad pairs of the delay spectra across the two JD-bins according to Equation (<ref>)
where,

    Ψ̃_∇(τ) . Ψ̃_∇^'(τ) = 1/N_ triadsC_2∑_i,j^N_ triadsΨ̃_∇(i, τ) . Ψ̃_∇(j, τ),    i>j,

where i, j runs over N_ triads (upper triangle of the {i,j} pairs) from the first and second JD bins. The normalising factor of 2 arises since phases only capture half the power in the fluctuations. After this operation, the data structure of the binned averaged cross P_∇(k_||) becomes 

    {N_ LST, N_ JDC_2, N_ pol, N_ triadsC_2, N_ delays}; N_ JDC_2=1.


We took the weighted mean (i.e., incoherently averaged) along the triad pairs, where the weights were propagated from the previous step (refer to data flowchart <ref> for details).
As the sky varies with LST, we applied the inverse variance weights along the LST axis and averaged them to get the final estimates of the power spectrum. The same operation was done for the imaginary part of the data to get an estimate of the level of systematics in the power spectrum. 


Ultimately, we incoherently averaged the two polarisations and downsampled the original delays to the effective bandwidth of the applied window function (i.e., ≈ 12.9 MHz). We used -based  to interpolate at the new downsampled delays. The final estimates of the closure phase power spectra are shown in figure <ref> and <ref> for the EoR0 and EoR1 fields, respectively. 
 


 §.§ Error estimation
 

The uncertainties on the power spectrum can be estimated in multiple ways. We primarily focused on estimating the uncertainties in two ways, the first being `noise+systematics' and the second only noise, where we tried to mitigate the systematics. 



  §.§.§ Noise+systematics
 

To estimate uncertainties in power, we increased the number of samples that go into the uncertainty estimation by splitting the JD-axis into four parts (i.e., N_ JD = 4). This led to the data (Ψ_∇(τ)) structure  being {N_ LST, N_ JD=4, N_ pol, N_ triads, N_ delays}. 
We similarly estimated the cross-power of the N_ triads along the unique pairs of N_ JD axis. This operation provided us with {N_ pol, N_ LST, N_ JDC_2=6, N_ triadsC_2, N_ delays} unique power spectra. The weighted mean power was estimated by along N_ triadsC_2 where the weights are coming from the number of good observations that went into the Ψ_∇ for a given triad.

Next, We estimated the weighted variance on the power using the standard error of the weighted mean provided in <cit.>,

    ( SEM_ wtd)^2 = n/(n-1)(∑ w_i)^2 [ ∑ (w_i X_i - w̅X̅_ wtd)^2 
    
        - 2X̅_ wtd∑ (w_i -w̅) (w_i X_i - w̅X̅_ wtd)
        + X̅_ wtd^2 ∑ (w_i -w̅)^2]   ,

where w_i are the weights, and n represents the weight count. Note that since the N_ JD = 4, we are required to normalise the variance. Figure <ref> illustrates the detailed data structure flow for the noise estimation.










  §.§.§ Noise
 

The uncertainties for the only noise case follow a similar procedure as the previous one, with the same JD split (i.e., N_ JD = 4), the cross-power is estimated, which led to the data structure {N_ pol, N_ LST, N_ JDC_2=6, N_ triadsC_2, N_ delays}. 
After this operation, we took the difference in the power spectra between the unique pairs of JD (i.e., the same sky signal). 

The only noise scenario can be understood assuming the cross-power between the two independent JD-bins correlates with the common signal and systematics across the triads. Assuming the power is coherent within the LST bin, the successive unique difference in the power within the LST bin eliminates the correlated power and systematics, leaving only the noise-like residuals. Also, since the differences eliminate the sky signal, the LST variation in the difference power is minimal; thus, we collapsed our datasets into a single axis and estimated the weighted standard deviation of the differenced power to get the final noise-like uncertainties. Finally, again, we took the weighted variance using eq. <ref>.







 §.§ Validation
 
We performed a two-sided KS test on the closure phase power spectra of the data and two model variants for the statistical comparison. The null hypothesis was rejected in all scenarios with the Model (without baseline dependent gains); however, it failed to reject the null-hypothesis in all scenarios when using the Model with g_ij. 
The former implies that the data and the model uncertainties were unlikely to be drawn from the same distribution, whereas the latter concludes the contrary. The test results are provided in table <ref>. The test results for the Model are not unexpected since we can see that the RMS floor levels of the data and Model are sufficiently different, which match the data and Model with g_ij. 
We modelled the excess RMS in the data as arising from baseline-dependent gain factors, although it is believed to have originated from the systematics and residual RFI. Note that we do not claim that the excess noise in the data is solely due to the baseline dependent systematics; however, if the argument is true, the baseline dependent gains introduced in our analysis suffice for the excess power in the data.






 §.§ 21-cm power spectrum

We estimated the final dimensionless 21-cm power spectrum from the closure phase power spectrum. The closure phase delay power spectrum can be written into a 21-cm power spectrum (“pseudo”) as follows:

    Δ_∇^2(k) = k^3P_∇(k_||)/2π^2 [pseudo  mK^2]

where, k^2 = k_⊥^2 + k_||^2,  with k_⊥ = 2π |b_∇|/λ D, where b_∇ is the baseline length of the triad, and D is the cosmological comoving distance. Note that the 21-cm power spectrum estimates from the closure phase power spectrum should not be interpreted as true but rather the approximate representation of the actual 21-cm power spectrum <cit.>. The power spectra converted to cosmological units for EoR0 and EoR1 fields are shown in Figures <ref> and <ref>, respectively.  



 

Assuming the convergence to normal distribution due to the central limit theorem, we estimated 2σ [95% confidence intervals (CI)] using the uncertainties since our sample size was sufficient (> 30). The upper limits on the 21-cm power spectrum [pseudo  mK^2] were then estimated 
{Δ_∇ UL^2 = (μ_Δ^2_∇± CI)  [pseudo  mK^2]}
for both the EoR0 and EoR1 fields are provided in the table <ref>, and the full table in the (<ref>).



§ DISCUSSION
 



We used the closure phase delay spectrum technique to obtain an independent estimate of the 21-cm power spectrum for the MWA phase II observations. These observations were centered on the EoR0 and EoR1 fields and were zenith-pointed, similar to the observing strategy of HERA. Our analysis revealed that MWA observations are possibly suffer from a baseline-dependent bandpass structure, which is especially noticeable in the edge channels. This bandpass structure results in structured bumps in the delay power spectrum (see figure <ref>), significantly contaminating the power spectrum. To address this issue, we explored two methods: Gaussian Process Regression (GPR) and Non-uniform Fast Fourier Transform (NFFT), to inpaint and mitigate the impact of the bandpass-affected edge channels on our power spectrum estimation. However, we decided against adopting the NFFT method because, although it reduced the magnitude of the bandpass, the bandpass remained evident in the NFFT delay spectra (see figure <ref>). Finally, we estimated the  21-cm power spectra using closure phase delay spectra. Additionally, we performed forward modelling in parallel with the observations to gain insights into the nature of the power spectrum under ideal observing conditions. The main findings of our analysis are summarised below.

    When we averaged closure phases across multiple timestamps within the same Local Sidereal Time (LST) bin, we noticed a significant residual bandpass structure, particularly noticeable in the edge channels (see Fig. <ref>). Since closure phases are unaffected by element-based bandpass variations, we concluded that these bandpass issues cannot be simplified into element-based terms and could instead be dependent on the baseline. To test this hypothesis, we simulated the same effect on foreground visibilities (see figure <ref>). We explored data-inpainting techniques to address these systematic errors to estimate closure phases in channels contaminated by baseline-dependent bandpass systematics (see figure <ref>). It is important to note that while these baseline-dependent issues are most noticeable in the edge channels, they could potentially affect all frequency channels since closure phases do not eliminate them. These issues also impact standard visibility-based power spectrum analysis methods. Understanding how the antenna layout contributes to such systematic errors is crucial for executing baseline-based mitigation strategies. Further investigation is needed to fully understand the extent to which these systematic errors affect MWA EoR power spectrum estimates. With a simple baseline dependent gains in the forward modelling, we aim to address the anomalies present in the DATA.
    
    On comparing the final closure phase power spectrum of the DATA and Model for the EoR0 field (see figure <ref>), we found that the peak power (at τ=0μ s) (i.e. ≈ 10^14 pseudo  mK^2h^-3 Mpc^3) of the DATA and Model only roughly matches for the 14 m triads, however the same for 24 m and 28 m triads differ significantly. During the initial closure phase estimation stage, we found EoR0 data having multiple phase wraps, which could be due to the presence of some systematics or residual RFI. It caused an overall shift in the peak power away from zero delays in the coherent averaging. This effect can be seen in the 28 m triad, which shows greater power next to the zeroth delay (see figure <ref>). The RMS floor level is between 1-2 orders of magnitude higher in the DATA (≈ 10^8 pseudo  mK^2h^-3 Mpc^3) compared to the Model (≈ 10^7 pseudo  mK^2h^-3 Mpc^3) and Model with g_ij (≈ 10^8 pseudo  mK^2h^-3 Mpc^3). This excess power in the data compared to the Model may arise from a smaller DATA sample size in the EoR0 field or systematics and residual RFI. Using a baseline-dependent gain factor in the simulation, we aimed to incorporate such systematics. We performed a 2-sided KS test on the DATA and Model at k_|| > 0.15  pseudo  h  Mpc^-1, which shows rejection of the null hypothesis for the likelihood of DATA and Model drawn from the same distribution at all baseline cases, which is expected since both differ significantly. In contrast, the KS-test setifies the null hypothesis when comparing DATA and Model with g_ij.In the closure phase power spectrum of the EoR1 field, the peak power of the DATA and Model (≈ 10^15 pseudo  mK^2h^-3 Mpc^3) match for all triads. The RMS floor level between the Model and DATA gets significantly better compared to the EoR0 field. They nearly match in all cases, except for the 14 m triads where the difference is approximately an order of magnitude higher in the DATA compared to the Model (see figure <ref>). It shows that we can improve our estimates of the power spectrum with an increased number of datasets. Thus, the analysis is data-limited for the amount of the data used. However, similar to the EoR0 case, the 2-sided KS test rejected the null hypothesis for all cases on Model, whereas fail to reject the null hypothesis in favour of the two distributions Model with g_ij might be drawn from the same distribution.

    Since our observations lie in the middle of the DTV broadcasting band, we further investigated for residual RFI in our data. We shifted our entire analysis to the lower frequency band (167-177 MHz), avoiding the central DTV-affected band (although, as shown in figure no. 2 in <cit.>, the DTV RFI nearly covers the entire EoR high band observations).
This analysis can help understand whether the residual RFI or other systematics are present in the data. Note that since we reduced our bandwidth by nearly three, we reduced our sensitivity by the same factor in the delay power spectrum. Therefore, the direct comparison of the mean power at k_|| > 0.15 pseudo h  Mpc^-1 might not be valid with the previous result.
 The closure phase power spectrum for the shifted spectrum is shown in figure <ref>, and <ref>. We can see significant improvements in the peak power of the DATA and Model compared to previous results. However, the overall RMS level was increased by an order of magnitude in all cases, possibly due to the lesser sensitivity (lower sample size). The DATA peak power at τ=0 μ s, especially in the EoR0 field, now matches the Model. Thus, we can justifiably argue that DTV RFI, which is expected to be prominent near 180 MHz, significantly contributed to the systematics present in the EoR0 DATA. On the other hand, EoR1 DATA seem relatively similar in both analyses, thus indicating that the systematics (such as persistent RFI) other than DTV RFI might be present in the data. Our findings are also confirmed when performing the KS-test on the DATA and Model, which shows non-rejection of the null hypothesis in all EoR1 field scenarios.
    We also compared the results with Model with g_ij which are shown in table <ref> shows KS-test outcomes of the Data and Model and Data and Model with g_ij on the shifted spectrum.
    
    We estimated the 2σ (95% confidence interval assuming Gaussianity from the convergence to Central Limit Theorem) on the 21-cm power spectrum for both the EoR0 and EoR1 fields. Our best upper limit on the 21-cm power spectrum of ≲ (184 pseudo  mK)^2 came from EoR1 field on 14 m triads at k=0.36 pseudo h Mpc^-1 with the noise only uncertainties. In the EoR0 field, our best estimate, ≲ (188 pseudo  mK)^2, came from the 24 m triads at k=0.18 pseudo h Mpc^-1 again using the noise-only uncertainty. However, as we discussed earlier, the systematics or residual RFI might have still affected these estimates, which we aim to address by introducing baseline-dependent gains in the modelling. It should be noted that, the exact nature of such baseline-dependent gains is not well understood. We have seen that, unlike Foregrounds, which usually gets restricted in lower delay modes, allowing faint HI signal to fluctuate visible at higher delay modes, the systematics equally affect all delay modes. Thus, for the scientific goal of observing milliradian-level sensitivity could be a significant challenge if such baseline dependent gains are present in the DATA. However, with extensive coherent averaging, the effect of such can be reduced. It should also be noted that the exact description of anomalies in the DATA can not be solely due to baseline-dependent systematics. Therefore, we state that if only the baseline-dependent systematics is present in the data, the level of noise introduced by the baseline-dependent gains justifies our DATA. Nevertheless, the level of fiducial HI and FG+HI powers are still lower than the data by 4-5 and 3-4 orders of magnitude, respectively; however, since our analysis is still data-limited, there is significant scope for improving upon the current estimates. Table <ref> provides the best 2σ estimates while table <ref> provides all estimates for each triad studied here.



§ SUMMARY

We present independent EoR 21-cm power spectrum results from the closure phase analysis of ≈ 12 hours of MWA phase-II data in the frequency range 167-197 MHz on three redundant baseline groups, namely, 14 m, 24 m, and 28 m baselines. Using the closure phase diagnostic, we found evidence for significant baseline-dependent systematics in the MWA data. Our best estimates of the 21-cm power spectrum at z=6.79 are ≲ (184 pseudo  mK)^2 at k=0.36 pseudo h Mpc^-1 in the EoR1 field using 14 m baseline triads and ≲ (188 pseudo  mK)^2 at k=0.18 pseudo h Mpc^-1 in the EoR0 field using 24 m baseline triads. Even with the limited amount of data analyzed, our closure phase method shows significant promise in independently constraining the 21-cm power spectrum during the EoR. Our results are still data-limited; hence, there is scope for further improvement by including more data in the analysis. Extensive sky modelling, such as accounting for the Galactic diffuse emission, is required before directly comparing the closure phase analysis with the standard visibility-based power spectrum. 



§ ACKNOWLEDGEMENTS
 
This project is supported by an ARC Future Fellowship under grant FT180100321. This research was partially supported by the Australian Research Council Centre of Excellence for All Sky Astrophysics in 3 Dimensions (ASTRO 3D) through project number CE170100013. The International Centre for Radio Astronomy Research (ICRAR) is a Joint Venture of Curtin University and The University of Western Australia, funded by the Western Australian State government. The MWA Phase II upgrade project was supported by the Australian Research Council LIEF grant LE160100031 and the Dunlap Institute for Astronomy and Astrophysics at the University of Toronto. This scientific work uses Inyarrimanha Ilgari Bundara, the CSIRO Murchison Radio-astronomy Observatory. We acknowledge the Wajarri Yamatji people as the traditional owners and native title holders of the Observatory site. Support for the operation of the MWA is provided by the Australian Government (NCRIS) under a contract to Curtin University administered by Astronomy Australia Limited. We acknowledge the Pawsey Supercomputing Centre, which is supported by the Western Australian and Australian Governments. Data were processed at the Pawsey Supercomputing Centre. 



§ DATA AVAILABILITY
 
This project was developed using Python and dependent libraries mentioned in the main text. Our data processing pipeline is publically available at [<https://github.com/himmng/Closure_phase_analysis>], along with the final processed datasets. The core data used in this work will be made available upon reasonable request. 













§ APPENDIX
 


 §.§ Inner vs Outer triads
 
The core of MWA Hexagons may be more affected by mutual coupling than the tiles near the edges. Therefore, we tried to perform an independent check on the power differences between the inner and outer tile triads to quantify the cross-talk level in the data. We chose 14m baseline triads to get a higher count and only estimated the closure phase power spectrum from them. We used first and second outer layer of the Hexagon configuration as the outer triads, and third and fourth tile layers as inner triads (see right panel of figure <ref>). In scenarios, where either two tiles from outer while the third from inner, or vice-versa, we considered those triads to be outer or inner triads, respectively. Figure <ref> shows the closure phase power spectrum of the inner and outer triads at 14m baseline lengths. We observed the relative percentage difference between the RMS power estimated for k_||>0.15  pseudo  h  Mpc^-1 between inner triads is higher than the outer triads are about 41%, 34% for EoR0 and EoR1 fields, respectively. However, the relative difference with the Model RMS (figure <ref>, <ref>) ≈ 6%, -25% for inner and outer triads in EoR0 field, and ≈ 186%, 114% for inner and outer triads in EoR1 field. These when comparing with the mean RMS value of the Model (≈ 10^8, 10^7 pseudo  mK^2h^-3 Mpc^3 for E0R0, EoR1 fields) consistent with each other.





 §.§ Avoiding DTV

The observations are centred around 180 MHz, which overlap with the Australian Digital Television (DTV) band; therefore, despite the extensive RFI flagging, some RFI could remain in the final processed data. Therefore, a useful test would be to shift the spectral window to lower frequencies, avoiding the DTV band in the analysis.  

To do this, we worked with the first ≈ 10 MHz band (167–177 MHz) of the total 30.72 MHz bandwidth and estimated the cross-power spectrum using the same procedure. However, the effective bandwidth for this analysis was reduced to ≈ 3.3 MHz, thus reducing the sensitivity by the same factor. 
The idea was to check whether the floor level of the power spectra (refer to figure <ref>, <ref>), which shows the excess power in the DATA, reduces and matches with the Model in the shifted window. It would suggest that the RFI could be localised in the central portion of the band, which is one of the major contributors to the excess power.
However, the power spectra of shifted window (see figure <ref> and <ref>) show similar behavior as that in figure <ref> and <ref>, indicating that the level of the RFI might be ubiquitous across the entire observing band along with the systematics in the data which are not being modelled in the simulations.








 §.§ Diagnosing bandpass systematics

We checked the closure phases for the baseline-dependent systematics, which persisted in the MWA data by modifying the antenna element-based gains (𝐠_i' s eq. <ref>).  
First, we created one set of MWA bandpass using random Gaussian between [0, 1] for each edge-channel frequency. We multiplied these by each visibilities in the correct parity pair order. It provides us with modified visibilities with randomised gains at the edge channels, which mimics the bandpass-affected visibilities. We used simulated flux densities for this procedure since they were produced ideally with unity antenna gains (although it does not affect having unity gains in the first place). 

Second, we used two scenarios to modify the bandpass further. In the first, the antenna-element-based gains were modified with new gains multiplied by the existing ones. This step verifies how the individual antenna-based gains vanish in the closure phase, illustrated in the top panel of figure <ref>. In the second scenario, instead of multiplicative gains from individual antenna elements (e.g. 𝐠_i, 𝐠_j), we multiplied an additional baseline-dependent term (𝐠_ij) that is not factorisable into element-based terms. It demonstrated that baseline-dependent gains do not cancel in the closure phase; see the same figure <ref>  bottom panel, where the residual difference between the residuals between the closure phase modified with 𝐠_ij and the original does not vanish.





 §.§ Incomplete modelling

Realistic vs. ideal beams or using fewer foreground sources can affect the final power estimates. Thus, we produced two test foreground simulations, one with the same 20,000 foreground sources but with ideal beam conditions where all dipole gains are active and set to unity (FGI_20k) and a second with only 5,000 foreground sources and real beam conditions, dead/missing dipoles incorporated in the beam evaluation (FGR_5k), which we compared against 20,000 foreground sources but with real beam conditions (FGR_20k).
Note that in the main results, we used 20,000 foreground sources with real beams, which are compared against the two test scenarios. The three cases' final closure phase power spectrum, foreground with real dipole gains with 20,000 sources,  foreground with unity dipole gains with 20,000 sources, and foreground with real dipole gains with 5,000 sources, are shown in figure <ref>. We obtained RMS power at ≥ 1.0  pseudo h Mpc^-1 to differentiate the two scenarios with the foreground with real dipole gains with 20,000 sources. The relative percentage error for unity dipole gains (20,000 sources) at 14m, 24m, and 28m baselines are 4%, 39%, 0.3%, and for real dipole gains (5,000 sources) at 14m, 24m, 28m baselines are 1.7%, 1.9%, 40%, respectively. Comparing with the RMS value of the Model from figure <ref>, <ref>) (≈ 10^8, 10^7 pseudo  mK^2h^-3 Mpc^3 for E0R0, EoR1 fields), these relative differences in both the scenarios are sufficiently less, thus, consistent with the Model.




 §.§ Data Structure Flowchart


As we proceed through the different averaging steps, the data structure is shown as a flowchart in figure <ref>.


 


