A Tapered Gridded Estimator (TGE) for the Multi-Frequency Angular Power Spectrum (MAPS) and the Cosmological HI 21-cm Power Spectrum
Somnath Bharadwaj, Srijita Pal, Samir Choudhuri, Prasun Dutta
aa r X i v : . [ a s t r o - ph . C O ] D ec MNRAS , 000–000 (0000) Preprint 24 December 2018 Compiled using MNRAS L A TEX style file v3.0
A Tapered Gridded Estimator (TGE) for theMulti-Frequency Angular Power Spectrum (MAPS)and the Cosmological HI 21-cm Power Spectrum
Somnath Bharadwaj ⋆ , Srijita Pal , Samir Choudhuri and Prasun Dutta Department of Physics, & Centre for Theoretical Studies, IIT Kharagpur, Kharagpur 721 302, India National Centre For Radio Astrophysics, Post Bag 3, Ganeshkhind, Pune 411007, India Department of Physics, IIT (BHU), Varanasi 221005, India
ABSTRACT
In this work we present a new approach to estimate the power spectrum P ( k ) of redshifted HI 21-cm brightness temperature fluctuations. The MAPS C ℓ ( ν a , ν b ) completely quantifies the second order statistics of the sky signal un-der the assumption that the signal is statistically homogeneous and isotropicon the sky. Here we generalize an already existing visibility based estimatorfor C ℓ , namely TGE, to develop an estimator for C ℓ ( ν a , ν b ) . The 21-cm powerspectrum is the Fourier transform of C ℓ (∆ ν ) with respect to ∆ ν = | ν a − ν b | ,and we use this to estimate P ( k ). Using simulations of 150 MHz GMRT obser-vations, we find that this estimator is able to recover P ( k ) with an accuracy of5 −
20% over a reasonably large k range even when the data in 80% randomlychosen frequency channels is flagged. Key words: methods: statistical, data analysis - techniques: interferometric-cosmology: diffuse radiation
Measurements of the cosmological HI 21-cm power spectrum can be used to probe the largescale distribution of neutral hydrogen (HI) across a large redshift range from the Dark Agesto the Post-Reionization Era (e.g. Bharadwaj & Ali 2005; Furlanetto, Oh & Briggs. 2006; ⋆ Email:[email protected] (cid:13)
S. Bharadwaj et al.
Morales & Wyithe 2010; Prichard & Loeb 2012; Mellema et al. 2013). Being very faint innature, the 21-cm signal is buried in foregrounds which are four to five orders of magnitudelarger than the expected signal (Shaver et al. 1999; Santos et al. 2005; Ali, Bharadwaj & Chengalur2008; Bernardi et al. 2009; Ghosh et al. 2012; Iacobelli et al. 2013; Choudhuri et al. 2017a).There are several ongoing and future experiments, e.g. Donald C. Backer Precision Array toProbe the Epoch of Reionization (PAPER , Parsons et al. 2010), the Low Frequency Array(LOFAR , var Haarlem et al. 2013; Yatawatta et al. 2013), the Murchison Wide-field Array(MWA , Bowman et al. 2013; Tingay et al. 2013), the Giant Metrewave Radio Telescope(GMRT, Swarup et al. 1991) the Square Kilometer Array (SKA1 LOW , Koopmans et al.2015) and the Hydrogen Epoch of Reionization Array (HERA , DeBoer et al. 2017) whichare aiming to detect the 21-cm power spectrum from the Epoch of Reionization (EoR).The biggest challenge for a detection of the redshifted 21-cm signal are the foregroundswhich include point sources, the diffuse Galactic synchrotron emission, the free-free emissionfrom our Galaxy and external galaxies. Various techniques have been proposed to overcomethis issue. The foreground subtraction technique proposes to subtract a foreground modelfrom the visibility data or the image and use the residual data to detect the 21-cm powerspectrum (Jeli´c et al. 2008; Bowman et al. 2009; Paciga et al. 2011; Chapman et al. 2012;Trott et al. 2012; Paciga et al. 2013; Trott et al. 2016). Considering P ( k ⊥ , k k ), the cylin-drical power spectrum of the 21-cm brightness temperature fluctuations, the foregroundsare expected to be primarily confined to a wedge in the ( k ⊥ , k k ) plane. Here, k ⊥ ans k k refer to the components of the 3-dimensional wave vector k perpendicular and par-allel to the line of sight direction respectively. The foreground avoidance technique pro-poses to use the region outside this “Foreground Wedge” to estimate the 21-cm powerspectrum (Datta et al. 2010; Parsons et al. 2012; Vedantham et al. 2012; Pober et al. 2013;Thyagarajan et al. 2013; Parsons et al. 2014; Pober et al. 2014; Liu et al. 2014a,b; Dillon et al.2014, 2015; Ali et al. 2015).A large variety of estimators have been proposed and applied to measure the power spec-trum of the brightness temperature fluctuations using the visibility data measured in radiointerferometric observations. Image-based estimators (Seljak 1997; Paciga et al. 2013) have http://astro.berkeley.edu/dbacker/eor http://reionization.org/ MNRAS000
Morales & Wyithe 2010; Prichard & Loeb 2012; Mellema et al. 2013). Being very faint innature, the 21-cm signal is buried in foregrounds which are four to five orders of magnitudelarger than the expected signal (Shaver et al. 1999; Santos et al. 2005; Ali, Bharadwaj & Chengalur2008; Bernardi et al. 2009; Ghosh et al. 2012; Iacobelli et al. 2013; Choudhuri et al. 2017a).There are several ongoing and future experiments, e.g. Donald C. Backer Precision Array toProbe the Epoch of Reionization (PAPER , Parsons et al. 2010), the Low Frequency Array(LOFAR , var Haarlem et al. 2013; Yatawatta et al. 2013), the Murchison Wide-field Array(MWA , Bowman et al. 2013; Tingay et al. 2013), the Giant Metrewave Radio Telescope(GMRT, Swarup et al. 1991) the Square Kilometer Array (SKA1 LOW , Koopmans et al.2015) and the Hydrogen Epoch of Reionization Array (HERA , DeBoer et al. 2017) whichare aiming to detect the 21-cm power spectrum from the Epoch of Reionization (EoR).The biggest challenge for a detection of the redshifted 21-cm signal are the foregroundswhich include point sources, the diffuse Galactic synchrotron emission, the free-free emissionfrom our Galaxy and external galaxies. Various techniques have been proposed to overcomethis issue. The foreground subtraction technique proposes to subtract a foreground modelfrom the visibility data or the image and use the residual data to detect the 21-cm powerspectrum (Jeli´c et al. 2008; Bowman et al. 2009; Paciga et al. 2011; Chapman et al. 2012;Trott et al. 2012; Paciga et al. 2013; Trott et al. 2016). Considering P ( k ⊥ , k k ), the cylin-drical power spectrum of the 21-cm brightness temperature fluctuations, the foregroundsare expected to be primarily confined to a wedge in the ( k ⊥ , k k ) plane. Here, k ⊥ ans k k refer to the components of the 3-dimensional wave vector k perpendicular and par-allel to the line of sight direction respectively. The foreground avoidance technique pro-poses to use the region outside this “Foreground Wedge” to estimate the 21-cm powerspectrum (Datta et al. 2010; Parsons et al. 2012; Vedantham et al. 2012; Pober et al. 2013;Thyagarajan et al. 2013; Parsons et al. 2014; Pober et al. 2014; Liu et al. 2014a,b; Dillon et al.2014, 2015; Ali et al. 2015).A large variety of estimators have been proposed and applied to measure the power spec-trum of the brightness temperature fluctuations using the visibility data measured in radiointerferometric observations. Image-based estimators (Seljak 1997; Paciga et al. 2013) have http://astro.berkeley.edu/dbacker/eor http://reionization.org/ MNRAS000 , 000–000 (0000) GE P ( k ⊥ , k k ), including the 3D TGE (Paper II), consider aFourier transform of the measured visibilities V ( U , ν ) along the frequency axis ν to obtainthe visibilities V ( U , τ ) in delay space τ (Morales & Hewitt 2004). This is used to estimate P ( k ⊥ , k k ). A difficulty arises if the data is missing or flagged in a few frequency channels inwhich case the delay channel visibilities V ( U , τ ) and the estimated power spectrum P ( k ⊥ , k k )are both modified by a convolution with the Fourier transform of the frequency samplingfunction. Missing or flagged channels are quite common in any typical observation due toa variety of reasons including man made radio frequency interference (RFI). The CHIPSestimator developed by Trott et al. (2016) overcomes this problem by using Least-SquaresSpectral Analysis (LSSA) to evaluate V ( U , τ ). However this needs to be applied individuallyfor each baseline, and the entire process could be computationally expensive for large datavolumes. In this paper we propose an alternative approach to estimate P ( k ⊥ , k k ) which is MNRAS , 000–000 (0000)
S. Bharadwaj et al. able to handle the problem of missing or flagged data with relative ease. Another point tonote is that the earlier estimators all introduce a frequency filter which smoothly goes tozero at the two edges of the frequency band. This is introduced to avoid a discontinuity atthe edges of the band, however it results in the loss of some signal. Such a filter is not neededin the new estimator proposed here.The multi-frequency angular power spectrum C ℓ ( ν a , ν b ) (MAPS; Datta et al. 2007,Mondal et al.2018) completely quantifies the second order statistics of the sky signal under the assump-tion that the signal is statistically homogeneous and isotropic on the sky. This however doesnot assume that the signal is ergodic or statistically homogeneous along the frequency axis.We have C ℓ ( ν a , ν b ) = C ℓ (∆ ν ) where ∆ ν = | ν a − ν b | if we impose the additional conditionthat the signal is ergodic along frequency. The 3D 21-cm power spectrum P ( k ⊥ , k k ) is theFourier transform of C ℓ (∆ ν ). In the new approach presented here we first estimate C ℓ (∆ ν )and use the binned C ℓ (∆ ν ) to estimate P ( k ⊥ , k k ). Even if some channels are missing, it isquite possible that the frequency separations ∆ ν are all present in the data. In this case itis quite straight forward to evaluate P ( k ⊥ , k k ) through a Fourier transform of C ℓ (∆ ν ). Moresophisticated techniques like the LSSA can be used in case some ∆ ν are missing, howeverthis needs to be applied to the binned C ℓ (∆ ν ) and the task is not computationally expensive.The MAPS C ℓ (∆ ν ) has been used to quantify the statistical properties of the back-ground radiation in GMRT observations at 150 MHz (Ali, Bharadwaj & Chengalur 2008;Ghosh et al. 2012) and 610 MHz (Ghosh et al. 2011a,b). The HI signal contribution to themeasured C ℓ (∆ ν ) is expected to decorrelate rapidly when ∆ ν is increased whereas the fore-ground contribution is expected to remain correlated for large ∆ ν separations. This propertywas used (Ghosh et al. 2011b) to model and remove the foreground contribution and obtaina residual C ℓ (∆ ν ) which is consistent with noise. It was thereby possible to place an ob-servational limit on the HI 21-cm power spectrum at z ≈ .
3. The estimator used in theseearlier works individually correlates pairs of visibilities to estimate C ℓ (∆ ν ), a techniquewhich is computationally expensive. The 2D TGE (Paper II) presents an efficient techniqueto estimate the angular power spectrum C ℓ . In Section 2. of this paper we have generalizedthis earlier work to develop an estimator for the MAPS C ℓ ( ν a , ν b ). In Section 3. we presenthow P ( k ⊥ , k k ) is obtained from the estimated C ℓ (∆ ν ). Section 4. presents the Simulationswhich we have used to validate our estimator, Section 5. presents the Results and Section6. presents the Discussion and Conclusions. MNRAS , 000–000 (0000) GE The 2D TGE, presented in Paper II considers radio-interferometric observations at a singlefrequency ν and uses the measured visibilities V i to estimate the angular power spectrum C ℓ of the background radiation at the frequency ν . Here V i refers to the i -th visibilitymeasurement with a corresponding baseline U i . The measured visibilities can be expressedas V i = ∂B∂T ! Z d U ˜ a ( U i − U ) ∆ ˜ T ( U ) + N i . (1)Here, the first term is the sky signal which is the convolution of ˜ a ( U ) and ∆ ˜ T ( U ) wherethese are the Fourier transforms of the primary beam A ( θ ) and the temperature fluctuationsin the sky δT ( θ ) respectively, and B = 2 k B T /λ is the Planck function in the Rayleigh-Jeanslimit. The second term N i is the system noise contribution.In order to taper the sky response, the measured visibilities are convolved with a function˜ w ( U ) which is the Fourier transform of a window function W ( θ ) which falls off to a valueclose to zero well before the first null of the telescope’s primary beam pattern (Paper I).Further, in order to reduce the computation, the convolved visibilities are evaluated on agrid in uv space using V cg = X i ˜ w ( U g − U i ) V i . (2)where the ‘c’ in V cg refers to “convolved” and g refers to different grid points with corre-sponding baselines U g . The sky response of V cg is tapered with the window function W ( θ ).Here we have used W ( θ ) = e − θ /θ w where the value of θ w = 57 ′ is chosen so as to suppressthe contribution from the outer regions and sidelobes of the telescope’s primary beam pat-tern (Figure 1 of Choudhuri et al. 2016a). For comparison, the full width half maxima ofthe 150 MHz GMRT primary beam pattern may be estimated to be 1 . λ/D = 157 ′ where D = 45 m is the antenna diameter.The convolved gridded visibilities can be expressed as V cg = ∂B∂T ! Z d U ˜ K ( U g − U ) ∆ ˜ T ( U ) + X i ˜ w ( U g − U i ) N i , (3) MNRAS , 000–000 (0000)
S. Bharadwaj et al. where ˜ K ( U g − U ) = Z d U ′ ˜ w ( U g − U ′ ) B ( U ′ )˜ a (cid:16) U ′ − U (cid:17) (4)is an effective “gridding kernel”, andB( U ) = X i δ D ( U − U i ) (5)is the baseline sampling function of the measured visibilities.The 2D TGE estimator is defined asˆ E g = M − g | V cg | − X i | ˜ w ( U g − U i ) | | V i | ! . (6)with h ˆ E g i = C ℓg where ℓ g = 2 πU g , and h i denotes an ensemble average over multiple real-izations of the sky brightness temperature fluctuations which are recorded in the visibilities.The second term in the brackets ( ... ) in eq. (6) is introduced to subtract out the noise biascontribution which arises due to the correlation of a visibility with itself. M g is a normal-ization factor which we shall discuss later. Simulations show that the 2D TGE provides anunbiased estimate of the angular power spectrum C ℓ (Paper II) while effectively suppress-ing the contribution from the sidelobes and outer regions of the telescope’s primary beam(Choudhuri et al. 2017b). M g Calculation
As discussed in Paper II, the normalization constant M g can be written as, M g = V g − X i | ˜ w ( U g − U i ) | V (7)where, V g = ∂B∂T ! Z d U | ˜ K ( U i − U ) | . (8)and V = ∂B∂T ! Z d U | ˜ a ( U i − U ) | . (9)The values of M g (eq. 7) depend on the baseline distribution (eq. 5) and the form of thetapering function W ( θ ), and it is necessary to calculate M g at every grid point in the uv plane. Paper I presents an analytic approximation to estimate M g . While this has been foundto work very well in a situation where the baselines have a nearly uniform and dense uv coverage (Fig. 7 of Paper I), it leads to C ℓ being overestimated in a situation where we havea sparse and non-uniform uv coverage. Paper II presents a different method to estimate M g which has been found to work well even if the uv coverage is sparse and non-uniform . MNRAS000
As discussed in Paper II, the normalization constant M g can be written as, M g = V g − X i | ˜ w ( U g − U i ) | V (7)where, V g = ∂B∂T ! Z d U | ˜ K ( U i − U ) | . (8)and V = ∂B∂T ! Z d U | ˜ a ( U i − U ) | . (9)The values of M g (eq. 7) depend on the baseline distribution (eq. 5) and the form of thetapering function W ( θ ), and it is necessary to calculate M g at every grid point in the uv plane. Paper I presents an analytic approximation to estimate M g . While this has been foundto work very well in a situation where the baselines have a nearly uniform and dense uv coverage (Fig. 7 of Paper I), it leads to C ℓ being overestimated in a situation where we havea sparse and non-uniform uv coverage. Paper II presents a different method to estimate M g which has been found to work well even if the uv coverage is sparse and non-uniform . MNRAS000 , 000–000 (0000) GE M g is calculated for C ℓ estimationin eq. (6) . As discussed in Paper II, we proceed by constructing random realizations ofsimulated visibilities [ V i ] UAPS corresponding to a situation where the sky signal has an unitangular power spectrum (UAPS) C ℓ = 1. The simulated visibilities have exactly the samebaseline distribution as the actual observed visibilities. We then have (eq. 6) M g = h | V cg | − X i | ˜ w ( U g − U i ) | h| V i | ! i UPAS (10)which allows us to estimate M g . We average over N u independent realizations of the UPASto reduce the statistical uncertainty. The estimator ˆ E g provides an estimate of C ℓ at different grid points U g on the uv plane.We have binned the estimates in order to increase the signal to noise ratio and also reducethe data volume. The signal is assumed to be statistically isotropic on the sky whereby itis independent of the direction of U g . This allows us to average the C ℓ estimates within anannular region on the uv plane. We define the binned Tapered Gridded Estimator for bin a using ˆ E G ( a ) = P g w g ˆ E g P g w g . (11)where w g refers to the weight assigned to the contribution from any particular grid point.The choice w g = 1 assigns equal weightage to the value of C ℓ g estimated at each grid point,whereas w g = M g corresponds to a situation where the grid points which have a denserbaseline sampling (less system noise) would be given a larger weightage. The former wouldbe desireable if one wishes to optimize with respect to the cosmic variance whereas the latterwould be preferred to optimize with respect to the system noise contribution. The optimumchoice of w g to maximize the signal to noise ratio would depend on the window function andthe baseline distribution, and we plan to address this in future.The binned estimator has an expectation value¯ C ¯ ℓ a = P g w g C ℓ g P g w g (12)where ¯ C ¯ ℓ a is the average angular power spectrum at¯ ℓ a = P g w g ℓ g P g w g (13)which is the effective angular multipole for bin a . MNRAS , 000–000 (0000)
S. Bharadwaj et al.
The multi-frequency angular power spectrum C ℓ ( ν a , ν b ) (Datta et al. 2007) characterizes thejoint frequency and angular dependence of the statistical properties of the background skysignal. We decompose the brightness temperature fluctuations δT b ( ˆ n , ν ) in terms of sphericalharmonics Y m ℓ ( ˆ n ) using δT b ( ˆ n , ν ) = X ℓ,m a ℓ m ( ν ) Y m ℓ ( ˆ n ) (14)and define the multi-frequency angular power spectrum (hereafter MAPS) as C ℓ ( ν a , ν b ) = h a ℓ m ( ν a ) a ∗ ℓ m ( ν b ) i . (15)As discussed in Mondal et al. (2018), we expect C ℓ ( ν , ν ) to entirely quantify the secondorder statistics of the redshifted 21-cm signal.We now proceed to define a visibility based Tapered Gridded Estimator (TGE) for C ℓ ( ν a , ν b ). We generalize the analysis to consider visibility measurements V i ( ν a ) at mul-tiple frequency channels 1 a N c , each of width ∆ ν c , with N c channels that span abandwidth B bw . Here we allow for the possibility that several of the data are bad or missing.We assume that such data has been identified and flagged, and this information is storedusing a flagging variable F i ( ν a ) which has value 0 for the flagged data and value 1 otherwise.We then have V cg ( ν a ) = X i ˜ w ( U g − U i ) V i ( ν a ) F i ( ν a ) . (16)which allows us to define the Tapered Gridded Estimator (TGE) for C ℓ ( ν a , ν b ) asˆ E g ( ν a , ν b ) = M − g ( ν a , ν b ) R e V cg ( ν a ) V ∗ cg ( ν b ) − δ a,b X i F i ( ν a ) | ˜ w ( U g − U i ) | | V i ( ν a ) | ! . (17)where R e () denotes the real part, δ a,b is a Kronecker delta i.e. it is necessary to subtractthe noise bias only when the two frequencies are the same ( ν a = ν b ), and the noise in thevisibility measurements at two different frequencies ( ν a = ν b ) are uncorrelated.The TGE defined in eq. (17) provides an unbiased estimate of C ℓg ( ν a , ν b ) at the angularmultipole ℓ g = 2 πU g i.e. h ˆ E g ( ν a , ν b ) i = C ℓg ( ν a , ν b ) (18)We use this to define the binned Tapered Gridded Estimator for bin a ˆ E G [ a ]( ν a , ν b ) = P g w g ˆ E g ( ν a , ν b ) P g w g . (19)where w g refers to the weight assigned to the contribution from any particular grid point g . MNRAS , 000–000 (0000) GE M g ( ν a , ν b ) which roughlyaverages the visibility correlation V cg ( ν a ) V ∗ cg ( ν b ) across all the grid points which are sampledby the baseline distribution. The binned estimator has an expectation value¯ C ¯ ℓ a ( ν a , ν b ) = P g w g C ℓg ( ν a , ν b ) P g w g (20)where ¯ C ¯ ℓ a ( ν a , ν b ) is the bin averaged multi-frequency angular power spectrum (MAPS) at¯ ℓ a = P g w g ℓ g P g w g (21)which is the effective angular multipole for bin a .Paper II describes how we have estimated M g using UAPS simulations in the contextof observations at a single frequency. This has also been summarized in Section 2 of thispaper. Here we have extended the earlier analysis to simulate visibilities [ V i ( ν a )] UMAPS forwhich we have an unit multi-frequency angular power spectrum C ℓ ( ν a , ν b ) = 1. We alsoapply the same flagging variable F i ( ν a ) as the actual data to the simulated data. Using thesimulated visibilities [ V i ( ν a )] UMAPS and the actual flagging variable F i ( ν a ) in eq. (17), wehave an estimate of M g ( ν a , ν b ). We have used multiple realizations of the simulations toreduce the uncertainty in the estimated values of M g ( ν a , ν b ).We note that the estimator presented here does not take into account the fact that thebaselines U i = d i /λ (where d is the antenna spacing) and the primary beam pattern A ( θ , ν )both change with frequency and these are held fixed at the values corresponding to the cen-tral frequency ν c . While this may not have a very significant effect on the recovered 21-cmpower spectrum, it is very important for the foregrounds where this leads to the foregroundwedge (eg. Datta et al. 2010; Parsons et al. 2012; Vedantham et al. 2012). We note that thefrequency dependence of the baselines has been included in earlier versions of the MAPSestimator (Ali, Bharadwaj & Chengalur 2008; Ghosh et al. 2011a, 2012) which did not in-corporate gridding and tapering. It is possible to incorporate the frequency dependence ofthe baselines in the TGE by suitably scaling the baselines U i at the time of convolution andgridding (eq. 16), and we plan to address this in future work. P ( K ⊥ , K k )In order to estimate the 3D power spectrum P ( k ⊥ , k k ) we assume that the redshifted 21-cmsignal is statistically homogeneous (ergodic) along the line of sight (e.g. Mondal et al. 2018).We then have C ℓ ( ν a , ν b ) = C ℓ (∆ ν ) where ∆ ν = | ν b − ν a | i.e. the statistical properties of MNRAS , 000–000 (0000) S. Bharadwaj et al. the signal depends only on the frequency separation and not the individual frequencies. Inthe flat sky approximation, the power spectrum P ( k ⊥ , k k ) of the brightness temperaturefluctuations of the redshifted 21-cm signal is the Fourier transform of C ℓ (∆ ν ), and we have(Datta et al. 2007) P ( k ⊥ , k k ) = r r ′ Z ∞−∞ d (∆ ν ) e − ik k r ′ ∆ ν C ℓ (∆ ν ) (22)where k k and k ⊥ = ℓ/r are the components of k respectively parallel and perpendicular tothe line of sight, r is the comoving distance corresponding to ν c the central frequency ofour observations and r ′ (= dr/dν ) is evaluated at ν c . A brief derivation of eq. (22) is alsopresented in the Appendix of Mondal et al. (2018). In this paper we have used (eq. 22) toestimate P ( k ⊥ , k k ) from the MAPS C ℓ ( ν a , ν b ).First we impose the ergodic assumption on C ℓ ( ν a , ν b ) which has been estimated from thevisibility data using eq. (17) and binned using eq. (19,20 and 21). For a fixed ℓ and ∆ ν , weaverage over all the C ℓ ( ν a , ν b ) values for which | ν b − ν a | = ∆ ν to obtain C ℓ (∆ ν ). We thenhave C ℓ ( n ∆ ν c ) where − ( N c − n ( N c −
1) with C ℓ ( n ∆ ν c ) = C ℓ ( − n ∆ ν c ). We seethat C ℓ ( n ∆ ν c ) is a periodic function of n with period 2 N c −
2. We use the discrete Fouriertransform ¯ P ( k ⊥ , k k m ) = ( r r ′ ∆ ν c ) N c − X n = − N c +2 exp (cid:16) − ik k m r ′ n ∆ ν c (cid:17) C ℓ ( n ∆ ν c ) (23)with k k m = m × [ π/r ′ c ∆ ν c ( N c − P ( k ⊥ , k k m ) which is already binned in k ⊥ . Wehave further binned in k k m to obtain the Spherical Power Spectrum P ( k ), and the CylindricalPower Spectrum P ( k ⊥ , k k ). We have carried out simulation to validate the estimator presented here. We have simulated8 hours of 150 MHz Giant Meterwave Radio Telescope (Swarup et al. 1991) observationswith N c = 257 channels of width ∆ ν c = 62 . B bw ≈
16 MHz and integrationtime ∆ t = 16 s towards RA=10h46m00s and DEC=59 ◦ ′ ′′ . We note that the EoR 21-cmsignal is not expected to be ergodic over the 16 MHz bandwidth considered here due to theLight Cone effect (Mondal et al. 2018). However, we have not considered this effect here andassumed that the signal is ergodic. The sky signal, we assume, is entirely the redshifted HI 21-cm emission whose brightness temperature fluctuations are characterized by the 3D powerspectrum P m ( k ) = ( k/k ) n mK Mpc . For the purpose of this paper we have arbitrarily MNRAS000
16 MHz and integrationtime ∆ t = 16 s towards RA=10h46m00s and DEC=59 ◦ ′ ′′ . We note that the EoR 21-cmsignal is not expected to be ergodic over the 16 MHz bandwidth considered here due to theLight Cone effect (Mondal et al. 2018). However, we have not considered this effect here andassumed that the signal is ergodic. The sky signal, we assume, is entirely the redshifted HI 21-cm emission whose brightness temperature fluctuations are characterized by the 3D powerspectrum P m ( k ) = ( k/k ) n mK Mpc . For the purpose of this paper we have arbitrarily MNRAS000 , 000–000 (0000) GE PSfrag replacements∆ ν MHz ℓ = 2165 ℓ = 5255 Noise, 80% Flagging − . . . . . . . . C ℓ ( ∆ ν ) − m K ∆ ν MHz ℓ .
01 0 . Figure 1.
This shows C ℓ (∆ ν ) as a function of ∆ ν for two values of ℓ . The data points with 1 − σ error-bars are estimatedfrom 24 realizations of the simulations. Note that the ∆ ν = 0 points have been slightly shifted for convenience of plotting ona logarithmic scale. The lines show the theoretical predictions calculated by using the input model power spectrum P m ( k ) ineq. (24). chosen the values k = (1 . − / Mpc − and n = −
2. We have followed the procedureoutlined in Section 4 of Choudhuri et al. (2016b) to simulate visibilities V i ( ν a ) correspondingto different statistically independent realizations of the brightness temperature fluctuations.In addition to the sky signal, the visibilities also contain a system noise contribution. Wehave modelled the system noise contribution to the visibilities as Gaussian random variableswhose real and imaginary parts both have zero mean and variance σ N . For comparison wehave also estimated σ sky which is the same quantity for the simulated sky signal contribution.The ratio R = σ N /σ sky gives an estimate of the relative contribution of the system noisewith respect to the sky signal. In our simulations we have used R = 10 which corresponds toa situation where the noise contribution to an individual visibility is R = 10 times the skysignal contribution. We have generated 24 statistically independent realizations of both thesky signal and the system noise. The resulting 24 statistically independent realizations of thesimulated visibilities were used to estimate the mean and 1 − σ errors for the results presentedbelow. We have considered simulations both with and without flagging. For each baselinewe have generated random integers in the range 1 a N c and flagged the correspondingchannels. We have carried out simulations for various values of f FLAG (the fraction of flaggedchannels) in the range 0 f FLAG . U = d /λ and the primary beampattern A ( θ , ν ) have both been incorporated in the simulated visibilities. MNRAS , 000–000 (0000) S. Bharadwaj et al. PSfrag replacements C ℓ ( ∆ ν ) − m K ∆ ν M H z . ℓ Figure 2.
This shows C ℓ (∆ ν ) across the entire ℓ and ∆ ν range considered here for simulations with noise and 80% flagging.We have shifted the ∆ ν = 0 values to ∆ ν = 0 . The analysis here was restricted to baselines in the range 10 | U i | , uv plane was divided into 15 annular bins at equal logarithmic intervals for power spectrumestimation. This corresponds to the k ⊥ range 7 × − Mpc − to 2 .
03 Mpc − . Figure 1 showsthe binned power spectrum C ℓ (∆ ν ) at two different values of ℓ considering simulationswith noise and 80% flagging. For comparison we have also shown the theoretical predictioncorresponding to the input model power spectrum P m ( k ) calculated using C ℓ (∆ ν ) = 1 πr Z ∞ dk k cos( k k r ′ ∆ ν ) P ( k ⊥ , k k ) (24)which is the inverse of eq. (22). We see that the results from the simulations are in agreementwith the theoretical predictions. The results shown here are visually indistinguishable fromthe results from simulations with no noise and no flagging, or those with 20% ,
40% and 60%flagging, and we have not shown the other results here.We find that the value of C ℓ (∆ ν ) falls rapidly as ∆ ν is increased, and it has a value closeto zero for ∆ ν > C ℓ (∆ ν ) across the entire ℓ and ∆ ν range that we haveconsidered here. The results are visually indistinguishable even if we have no noise and noflagging (or less flagging), or if we evaluate C ℓ (∆ ν ) analytically using (eq. 24) and we havenot shown these here. We see that the value of C ℓ (∆ ν ) decrease as ℓ is increased. For a fixed ℓ , the value of C ℓ (∆ ν ) falls rapidly as ∆ ν is increased and it has a value close to zero at large MNRAS000
40% and 60%flagging, and we have not shown the other results here.We find that the value of C ℓ (∆ ν ) falls rapidly as ∆ ν is increased, and it has a value closeto zero for ∆ ν > C ℓ (∆ ν ) across the entire ℓ and ∆ ν range that we haveconsidered here. The results are visually indistinguishable even if we have no noise and noflagging (or less flagging), or if we evaluate C ℓ (∆ ν ) analytically using (eq. 24) and we havenot shown these here. We see that the value of C ℓ (∆ ν ) decrease as ℓ is increased. For a fixed ℓ , the value of C ℓ (∆ ν ) falls rapidly as ∆ ν is increased and it has a value close to zero at large MNRAS000 , 000–000 (0000) GE PSfrag replacements k k M p c − k ⊥ Mpc − P ( k ⊥ , k k ) m K M p c . . .
11 1 10 − − Figure 3.
The cylindrical power spectrum P ( k ⊥ , k k ) estimated from simulations with with noise and 80% flagging. ∆ ν . The decrease in the value of C ℓ (∆ ν ) with increasing ∆ ν is more rapid as we go to larger ℓ . The behaviour of C ℓ (∆ ν ) is directly manifested in the visibility correlation V ( U, ∆ ν ) = hV ∗ ( U , ν ) V ( U , ν + ∆ ν ) i with ℓ = 2 πU . This has been studied extensively in several earlierworks (Bharadwaj & Sethi 2001; Bharadwaj & Pandey 2003; Bharadwaj & Ali 2005), andwe do not discuss this any further here.Figure 3 shows the cylindrical power spectrum P ( k ⊥ , k k ) estimated by applying eq. (23)to the C ℓ (∆ ν ) obtained from the simulations with noise and 80% flagging (Figure 2). Weobtain estimates of P ( k ⊥ , k k ) in 15 bins of equal logarithmic spacing along k ⊥ each with N c = 257 values along k k . We have further binned P ( k ⊥ , k k ) into 16 bins of equal logarithmicspacing along k k to increase the signal to noise ratio and also for convenience of plotting. Theresults for the other cases which we have considered (lesser flagging, with/without noise)are very similar and they have not been shown separately. Note that the estimated powerspectrum turns out to be negative at a single pixel (top left corner of the figure) when wehave 80% flagging. In contrast, we obtain positive values at all the pixels when we considera smaller percentage of flagged data.The upper panel of Figure 4 shows the spherical power spectrum P ( k ) estimated fromthe simulations with no noise and no flagging, and also the simulations with noise and 80%flagging. Here the P ( k ⊥ , k k ) values were combined into 15 bins of equal logarithmic intervalin the k range 4 × − Mpc − to 3 Mpc − . We see that the results from these two sets ofsimulations are visually indistinguishable. The results for all the other cases considered hereare very similar and they have not been shown separately. The model power spectrum P m ( k ) MNRAS , 000–000 (0000) S. Bharadwaj et al.
PSfrag replacements P ( k ⊥ , k k ) mK Mpc . . − − P ( k ) m K M p c δ k Mpc − .
01 1 − . . . − . − . . . . . Figure 4.
The upper panel shows the estimated spherically-binned power spectrum P ( k ) and 1 − σ error-bars for simulationswith no noise and flagging and also with noise and 80% flagging. For comparison, the input model P m ( k ) is also shown by thesolid line. The bottom panel shows the fractional error δ = [ P ( k ) − P m ( k )] /P m ( k ) (data points) and the relative statisticalfluctuation σ/P m ( k ) (shaded regions). The values of σ are larger for simulations with noise and 80% flagging as compared tothose with no noise and no flagging. is also shown for comparison. We see that P ( k ) estimated from the simulations is below themodel predictions at k < .
02 Mpc − . Our estimator assumes that the convolution due to thetelescope’s primary beam pattern can be well approximated by a multiplicative factor whichwe have incorporated in M g (eq. 17). Earlier studies (Choudhuri et al. 2014) show that thisassumption does not hold at the small baselines (which also correspond to small k ⊥ ) thatprobe angular scales which are comparable to the angular extent of the telescope’s primarybeam pattern. The estimated power spectrum is in better agreement with the input model at k > .
02 Mpc − . We however notice that P ( k ) is somewhat overestimated at k > .
03 Mpc − ,but this difference goes down at larger k . The lower panel shows the fractional deviation δ = [ P ( k ) − P m ( k )] /P m ( k ) of the estimated power spectrum P ( k ) relative to the inputmodel P m ( k ), the shaded regions shows the 1 − σ errors σ/P m ( k ). We find that P ( k ) isoverestimated by 10 −
20% in the range 0 . k < . − , this falls to 5 −
15% in therange 0 . k < . − and the overestimate is less than 7 .
5% at k > . − . Thereis around ∼
1% difference in the estimated values when we have noise and 80% flaggingas compared to the situation when these are not incorporated. Further, the values of σ arelarger when we introduce noise and flagging, this is particularly more pronounced at large k . In all cases we find that the errors δ are less than the expected statistical fluctuations σ/P m ( k ). MNRAS000
1% difference in the estimated values when we have noise and 80% flaggingas compared to the situation when these are not incorporated. Further, the values of σ arelarger when we introduce noise and flagging, this is particularly more pronounced at large k . In all cases we find that the errors δ are less than the expected statistical fluctuations σ/P m ( k ). MNRAS000 , 000–000 (0000) GE The error-bars shown here are based on 24 independent realizations of the simulations.Choudhuri et al. (2014) and Choudhuri et al. (2016b) present analytical formulas for esti-mating the statistical errors, and it is possible to obtain similar formulas for C ℓ ( ν a , ν b ) andpropagate the resulting errors through the Fourier transform to predict errors for P ( k ). How-ever, simulations offer a more straight forward method to estimate the errors. It is possible touse the estimated power spectrum as an input for simulations, and use multiple realizationsof these simulations to estimate the error-bars for the estimated power spectrum.An earlier study (Choudhuri et al. 2014) used simulations to shows that the TGE over-estimates C ℓ due to the sparse and patchy uv coverage of the GMRT baseline distribution.This overestimate was found to come down if a more dense and uniform uv coverage wasconsidered instead. The estimator presented here overestimates P ( k ) by 5 −
20% across alarge portion of the k range, the exact cause for this is not known at present . We believethat this is a consequence of the sparse and patchy uv coverage of the GMRT baseline dis-tribution, and is not an inherent limitation of the estimator. We expect this effect to bemuch less severe for an array with a denser and more uniform uv coverage. Further studiesconsidering arrays with different uv coverage are needed to quantitatively establish this, andwe propose to address this in future work.The signal in the visibility measurements V ( U i , ν a ) at different baselines U i are not in-dependent due to the telescope’s primary beam pattern and the signal at baselines within D/λ are correlated (eq. (12) of Bharadwaj & Ali 2005) where D is the antenna diam-eter. Similarly, the visibility measurements V ( U i , ν a ) at different frequency channels ν a are not independent (Figure 9 of Bharadwaj & Ali 2005) and the signal remains corre-lated across different channels, the width of the correlation depending on the value of U i (Bharadwaj & Pandey 2003). The signal contained in the flagged data which is lost is alsocontained in the valid data which is available at our disposal for power spectrum estimation,and the estimator presented here is able to recover the power spectrum equally well even if80% of the data is flagged. While the estimated power spectrum is practically unchangedwith or without flagging, the statistical fluctuations σ are somewhat larger (particularly atlarge k ) when flagging is introduced. The entire analysis presented here is restricted to asituation where randomly chosen frequency channels were flagged. A variety of other situa-tions may occur in real life. For a given real data it would be best to first use the flagging MNRAS , 000–000 (0000) S. Bharadwaj et al. variables of the actual data in conjunction with simulations to verify if the estimator canreproduce the input model of the simulation. If needed, the discrete Fourier transform ofeq. (23) can be replaced by a more sophisticated spectral estimator. However, here it is nec-essary to apply this to the final binned data and not the individual baselines, and thereforethe problem is not computationally demanding. We propose to address these issues in moredetail in future work.In a recent paper Morales et al. (2018) has broadly classified the power spectrum esti-mators into two classes namely (1.) the delay spectrum or measured sky estimators, and(2.) the reconstructed sky estimators. The former class of estimators performs the Fouriertransform from ν to k k at a fixed antenna separation d which does not incorporate thefrequency dependence of the baseline. In contrast, the same Fourier transform is carried outat the baseline U corresponding to a fixed angular scale which effectively incorporates thevariation of baseline with frequency, however it uses a reconstructed sky model instead ofthe measured sky signal. The estimator presented here deals with the measured sky signal,it however differs from the usual delay spectrum estimators in that the signal is first corre-lated and then Fourier transformed. It is consequently possible to incorporate the frequencydependence of the baselines (as mentioned in Section 3). This has not been incorporated inthe present work, we plan to incorporate this and study its impact on foregrounds in futurework. REFERENCES
Ali S. S., Bharadwaj S.,& Chengalur J. N., 2008, MNRAS, 385, 2166AAli, Z. S., Parsons, A. R., Zheng, H., et al. 2015, ApJ, 809, 61Bernardi, G., de Bruyn, A. G., Brentjens, M. A., et al. 2009, A&A, 500, 965Bharadwaj , S. & Sethi , S. K., 2001 , JApA, 22 , 293Bharadwaj, S., & Pandey, S. K. 2003,JAPA , 24, 23.Bharadwaj S. , & Ali S. S. 2005, MNRAS, 356, 1519Bowman, J. D., Morales, M. F., & Hewitt, J. N. 2009, ApJ, 695, 183Bowman J. D. et al., 2013, PASA, 30, e031Chapman, E., Abdalla, F. B., Harker, G., et al. 2012, MNRAS, 423, 2518Choudhuri, S., Bharadwaj, S., Ghosh, A., & Ali, S. S., 2014, MNRAS, 445, 4351Choudhuri, S., Bharadwaj, S., Roy, N., Ghosh, A., & Ali, S. S., 2016a, MNRAS, 459, 151Choudhuri, S., Bharadwaj, S., Chatterjee, S., et al. 2016b, MNRAS, 463, 4093Choudhuri, S., Bharadwaj, S., Ali, S. S., et al. 2017a, MNRAS, 470, L11Choudhuri, S., Roy, N., Bharadwaj, S., et al. 2017b, New Astron., 57, 94Datta, K. K., Choudhury, T. R., & Bharadwaj, S. 2007, MNRAS, 378, 119Datta, A., Bowman, J. D., & Carilli, C. L. 2010, ApJ, 724, 526 MNRAS000
Ali S. S., Bharadwaj S.,& Chengalur J. N., 2008, MNRAS, 385, 2166AAli, Z. S., Parsons, A. R., Zheng, H., et al. 2015, ApJ, 809, 61Bernardi, G., de Bruyn, A. G., Brentjens, M. A., et al. 2009, A&A, 500, 965Bharadwaj , S. & Sethi , S. K., 2001 , JApA, 22 , 293Bharadwaj, S., & Pandey, S. K. 2003,JAPA , 24, 23.Bharadwaj S. , & Ali S. S. 2005, MNRAS, 356, 1519Bowman, J. D., Morales, M. F., & Hewitt, J. N. 2009, ApJ, 695, 183Bowman J. D. et al., 2013, PASA, 30, e031Chapman, E., Abdalla, F. B., Harker, G., et al. 2012, MNRAS, 423, 2518Choudhuri, S., Bharadwaj, S., Ghosh, A., & Ali, S. S., 2014, MNRAS, 445, 4351Choudhuri, S., Bharadwaj, S., Roy, N., Ghosh, A., & Ali, S. S., 2016a, MNRAS, 459, 151Choudhuri, S., Bharadwaj, S., Chatterjee, S., et al. 2016b, MNRAS, 463, 4093Choudhuri, S., Bharadwaj, S., Ali, S. S., et al. 2017a, MNRAS, 470, L11Choudhuri, S., Roy, N., Bharadwaj, S., et al. 2017b, New Astron., 57, 94Datta, K. K., Choudhury, T. R., & Bharadwaj, S. 2007, MNRAS, 378, 119Datta, A., Bowman, J. D., & Carilli, C. L. 2010, ApJ, 724, 526 MNRAS000 , 000–000 (0000) GE DeBoer, D. R., Parsons, A. R., Aguirre, J. E., et al. 2017, PASP, 129, 045001Dillon, J. S., Liu, A., Williams, C. L., et al. 2014, PRD, 89, 023002Dillon, J. S., Liu, A., Williams, C. L., et al. 2015, PRD, 91(12), 123011Furlanetto S. R., Oh S. P., Briggs F. H., 2006, Phys. Rep.,433, 181Ghosh, A., Bharadwaj, S., Ali, S. S., & Chengalur, J. N. 2011a, MNRAS, 411, 2426Ghosh, A., Bharadwaj, S., Ali, S. S., & Chengalur, J. N. 2011b, MNRAS, 418, 2584Ghosh, A., et al. 2012, MNRAS, 426, 3295Iacobelli, M., Haverkorn, M., Orr´u, E., et al. 2013, A&A, 558, A72Jeli´c, V., Zaroubi, S., Labropoulos, P., et al. 2008, MNRAS, 389, 1319Koopmans, L., Pritchard, J., Mellema, G., et al. 2015, Advancing Astrophysics with the Square Kilometre Array (AASKA14),1Liu, A., & Tegmark, M. 2012, MNRAS, 419, 3491Liu, A., Parsons, A. R., & Trott, C. M. 2014a, PRD, 90, 023018Liu, A., Parsons, A. R., & Trott, C. M. 2014b, PRD, 90, 023019Liu, A., Zhang, Y., & Parsons, A. R. 2016, ApJ, 833, 242McQuinn M., Zahn O., Zaldarriaga M., Hernquist L., Furlanetto S. R., 2006, ApJ, 653, 815Mellema, G., et al. 2013, Experimental Astronomy, 36, 235Mondal, R., Bharadwaj, S., & Datta, K. K. 2018, MNRAS, 474, 1390Morales, M. F., & Hewitt, J. 2004, ApJ, 615, 7Morales M. F., 2005, ApJ, 619, 678Morales M. F., Matejek M., 2009, MNRAS, 400, 1814Morales, M. F., & Wyithe, J. S. B. 2010, ARA&A, 48, 127Morales M. F., Beardsley A., Pober J., Barry N., Hazelton B., Jacobs D., Sullivan I., 2018, MNRAS, 2715Paciga G. et al., 2011, MNRAS, 413, 1174Paciga, G., Albert, J. G., Bandura, K., et al. 2013, MNRAS, 433, 639Parsons A. R. et al., 2010, AJ, 139, 1468Parsons, A. R., Pober, J. C., Aguirre, J. E., et al. 2012, ApJ, 756, 165Parsons, A. R., Liu, A., Aguirre, J. E., et al. 2014, ApJ, 788, 106Pen U.-L., et al., 2009, MNRAS, 399, 181Planck Collaboration, P. A. R. Ade et al., Planck 2015 results. XIII. Cosmological parameters, arXiv:1502.01589.Pober J. C. et al., 2013, ApJ, 768, L36Pober, J. C., Liu, A.,Dillon, J. S., et al. 2014, ApJ, 782, 66Pober, J. C., Hazelton, B. J., Beardsley, A. P., et al. 2016, arXiv:1601.06177Pritchard, J. R. and Loeb, A., 2012, Reports on Progress in Physics 75(8), 086901Santos, M.G., Cooray, A. & Knox, L. 2005, 625, 575Seljak, U. 1997, ApJ, 482, 6Shaver, P. A., Windhorst, R. A., Madau, P., & de Bruyn, A. G. 1999, A&A, 345, 380Swarup, G., Ananthakrishnan, S., Kapahi, V. K., Rao, A. P., Subrahmanya, C. R., and Kulkarni, V. K. 1991, CURRENTSCIENCE, 60, 95.Thyagarajan, N., Udaya Shankar, N., Subrahmanyan, R., et al. 2013, ApJ, 776, 6Thyagarajan, N., Jacobs, D. C., Bowman, J. D., et al. 2015, ApJ, 807, L28Tingay, S. et al. 2013, Publications of the Astronomical Society of Australia, 30, 7Trott, C. M., Wayth, R. B., & Tingay, S. J. 2012, ApJ, 757, 101Trott, C. M., Pindor, B., Procopio, P., et al. 2016, ApJ, 818, 139van Haarlem, M. P., Wise, M. W., Gunst, A. W., et al. 2013, A&A, 556, A2MNRAS , 000–000 (0000) S. Bharadwaj et al.
Vedantham, H., Udaya Shankar, N., & Subrahmanyan, R. 2012, ApJ, 745, 176Yatawatta, S. et al. 2013, Astronomy & Astrophysics, 550, 136 MNRAS000