Large-scale 21-cm Cosmology with LOFAR and AARTFAAC
Gehlot, Bharat Kumar
IMPORTANT NOTE: You are advised to consult the publisher's version (publisher's PDF) if you wish to cite from it. Please check the document version below.
Document Version
Publisher's PDF, also known as Version of record
Publication date: 2019
Link to publication in University of Groningen/UMCG research database
Citation for published version (APA):
Gehlot, B. K. (2019). Large-scale 21-cm Cosmology with LOFAR and AARTFAAC. University of Groningen.
Copyright
Other than for strictly personal use, it is not permitted to download or to forward/distribute the text or part of it without the consent of the author(s) and/or copyright holder(s), unless the work is under an open content license (like Creative Commons).
Take-down policy
If you believe that this document breaches copyright please contact us providing details, and we will remove access to the work immediately and investigate your claim.
Downloaded from the University of Groningen/UMCG research database (Pure): http://www.rug.nl/research/portal. For technical reasons the number of authors shown on this cover page is limited to 10 maximum.
4
974
AARTFAAC Cosmic Explorer: Pilot observations of
the 21-cm power spectrum in the EDGES absorption
trough
B. K. Gehlot, F. G. Mertens, L. V. E. Koopmans, A. R. Offringa, A. Shulevski, M. A. Brentjens, M. Kuiack, M. Mevius, V. N. Pandey, A. Rowlinson, A. M. Sardarabadi, H. K. Vedantham, R. A. M. J. Wijers, S. Yatawatta and S. Zaroubi
4
Abstract
Recently, the EDGES collaboration (Bowman et al. 2018) reported a tentative detection of the 21-cm absorption feature seen against the Cosmic Microwave Background (CMB) in the sky-averaged spectrum. The feature is several times stronger than that predicted by traditional astrophysical models. If genuine, it can lead to much stronger fluctuations on the 21-cm signal on degree scales (up to 1 Kelvin in rms). This will allow these fluctuations to be detected in integration times nearly fifty times shorter compared to previous predictions. To detect or set limits on the 21-cm signal fluctuations, we commenced the “AARTFAAC Cosmic Explorer” (ACE) program to measure or set limits on the power spectrum of the 21-cm fluctuations in the redshift range z = 17.9− 18.7 (∆ν = 72− 75 MHz) corresponding to the deepest part of the EDGES absorption feature. We use the Amsterdam-ASTRON Radio Transient Facility and Analysis Centre (AARTFAAC), which is a LOFAR wide-field imager. In this chapter, we present the first results from 6 h of data taken in ‘semi drift-scan’ mode early in the ACE program. We observe an excess factor of∼ 2.4 in the ratio of the difference of the Stokes I and V gridded visibilities (at 4 s level times scales) when averaged over 6 h. Finally, we present a first exploratory 21-cm power spectrum from ACE data analysed using the LOFAR-EoR data calibration and imaging pipeline and Gaussian Process Regression foreground removal. After 6 h of incoherent averaging, a systematics limited 2σ upper limit of ∆2
21< (71 K)2 at k = 0.145 h cMpc−1 is reached in the redshift range
4
4.1: Introduction 99
4.1
Introduction
Observations of the redshifted 21-cm signal of neutral hydrogen from the Cos-mic Dawn and Epoch of Reionization hold the potential to revolutionise our understanding of how these first stars and galaxies formed and the nature of their ionizing radiation (Madau, Meiksin & Rees 1997; Shaver et al. 1999; Furlanetto, Oh & Briggs 2006; Pritchard & Loeb 2012; Mesinger, Furlanetto & Cen 2011; Zaroubi 2013). During the Cosmic Dawn (CD) (12≲ z ≲ 30), the first luminous objects formed in the dark and neutral Universe (Pritchard & Furlanetto 2007). X-ray and ultraviolet radiation from these first stars heated and ionized neutral hydrogen (HI) in the surrounding Inter-Galactic Medium (IGM) during the Epoch of Reionization (EoR). This process continued (span-ning the redshift range 6 ≲ z ≲ 12) until hydrogen in the IGM became fully ionized (Madau, Meiksin & Rees 1997).
In recent years, a large number of observational efforts got underway to ob-serve this faint 21-cm signal from the CD and EoR. Radio interferometers such as the LOw Frequency ARray1(LOFAR; van Haarlem et al. 2013), the
Giant Meterwave Radio Telescope2(GMRT; Paciga et al. 2011), the Murchison Widefield Array3(MWA; Tingay et al. 2013; Bowman et al. 2013), the Precision
Array for Probing the Epoch of Reionization4(PAPER; Parsons et al. 2010)
as well as the next generation instruments such as the Hydrogen Epoch of Reionization Array5(HERA; DeBoer et al. 2017), the Long Wavelength Array 6,7(LWA; Greenhill et al. 2012), the NENUFAR8(New Extension in Nançay
Upgrading loFAR; Zarka et al. 2012), and the upcoming Square Kilometre Array9(SKA; Mellema et al. 2013; Koopmans et al. 2015) are working towards
measuring the spatial brightness temperature fluctuations in the cosmological 21-cm signal.
In parallel, single-element radiometers such as the Experiment to Detect the Global Epoch of Reionization Signature (EDGES; Bowman et al. 2018), the Large-aperture Experiment to Detect the Dark Ages (LEDA; Bernardi et al. 2016), the Shaped Antenna measurement of the background RAdio Spectrum 2 (SARAS 2; Singh et al. 2017), the Sonda Cosmológica de las Islas para la Detección de Hidrógeno Neutro (SCI-HI; Voytek et al. 2014), the
Prob-1 http://www.lofar.org/ 2 http://gmrt.ncra.tifr.res.in/ 3 http://www.mwatelescope.org/ 4 http://eor.berkeley.edu/ 5 http://reionization.org/ 6 http://www.tauceti.caltech.edu/LWA/ 7 http://lwa.phys.unm.edu/epo.html 8 https://nenufar.obs-nancay.fr/ 9 http://skatelescope.org/
4
ing Radio Intensity at high z from Marion (PRIZM; Philip et al. 2018), and the Netherlands-China Low frequency Explorer10,11 (NCLE) are seeking to measure the global sky-averaged 21-cm signal as a function of redshift. Recently, a deep spectral feature centred at 78 MHz was reported by the EDGES collaboration (Bowman et al. 2018). The feature was presented as the long sought-after 21-cm absorption feature seen against the CMB during the CD at z ∼ 17. The location of this putative absorption trough is con-sistent with redshift predictions from theoretical models and simulations of the Cosmic Dawn (Furlanetto, Oh & Briggs 2006; Pritchard & Loeb 2010b; Mesinger, Ferrara & Spiegel 2013; Cohen et al. 2017). However, the depth of the feature is ∆T21 ∼ 0.5 K (99% confidence level), which is 2 − 3 times
stronger and considerably wider (∆ν ∼ 19 MHz) than that predicted by the most optimistic astrophysical models (e.g. Pritchard & Loeb 2010b; Fialkov & Barkana 2014; Fialkov & Loeb 2016; Cohen et al. 2017). Moreover, the observed feature is flat-bottomed instead of a smooth Gaussian-like shape. Several “exotic” theoretical models have already been proposed which might explain the depth of the feature, such as a considerably colder IGM due to interaction between baryons and dark matter particles causing a lower spin-temperature and therefore a deeper absorption feature (e.g. Barkana 2018; Fialkov, Barkana & Cohen 2018), or a stronger radiation background against which the absorption is taking place (e.g. Feng & Holder 2018; Ewall-Wice et al. 2018; Dowell & Taylor 2018). Although the 21-cm signal is expected to be stronger at these redshifts, the foreground emission is several times brighter at these frequencies compared to EoR 21-cm signal observations at 150 MHz (Bernardi et al. 2009, 2010). Moreover, ionospheric effects are amplified at lower frequencies (Gehlot et al. 2018a; de Gasperin et al. 2018), making the measurement of the signal equally (or even more) challenging than in EoR experiments. As of now, Ewall-Wice et al. (2016) reported a systematics-limited power spectrum upper limit of ∆2
21 < (104mK)2 on co-moving scales
k ≲ 0.5 h cMpc−1 of the 21-cm signal brightness temperature in the redshift range 12≲ z ≲ 18 using MWA. This overlaps with the low redshift edge of the 21-cm absorption feature (Bowman et al. 2018). Gehlot et al. (2018b) recently provided a 2σ upper limit of ∆2
21< (1.4× 104mK)2 on the 21-cm signal power
spectrum at k = 0.038 h cMpc−1using the LOFAR-Low Band Antenna (LBA) system in the redshift range 19.8 ≲ z ≲ 25.2, which corresponds to the high redshift edge of the absorption feature.
Although concerns have been raised about the validity of the detection of the absorption feature in terms of foreground modelling and instrumental effects
10 https://www.ru.nl/astrophysics/research/radboud-radio-lab-0/projects/
netherlands-china-low-frequency-explorerncle/
4
4.1: Introduction 101
(Hills et al. 2018), if the detection is confirmed, the strength of the 21-cm absorption feature can also cause a large increase in the 21-cm brightness temperature fluctuations in the redshift range z = 17− 19 (Barkana 2018; Fialkov, Barkana & Cohen 2018). This redshift range corresponds to the deepest part of the absorption profile and enables detection of the 21-cm signal brightness temperature fluctuations on degree angular scales in this redshift range within a much shorter integration time (∼ 50 times shorter) compared to what was previously expected.
Motivated by this, we have commenced a large scale program called “AART-FAAC Cosmic Explorer” (ACE) to measure or limit the power spectrum of the brightness temperature fluctuations of the 21-cm signal from z∼ 18 using the LOFAR Amsterdam-ASTRON Radio Transients Facility And Analysis Centre (AARTFAAC) wide-field imager (Prasad et al. 2016). AARTFAAC correlates up to 576 individual receiver elements (LBA dipoles or High Band Antenna tiles) in the core of LOFAR, thereby providing a much larger field of view and more sensitivity on large angular-scales compared to regular LOFAR ob-servations. The redshift range targeted by ACE is z = 17.9− 18.6 which corresponds to 72.39− 75.05 MHz frequency range. The ACE programme12
aims for a 1000 h deep integration of a large part of the northern sky to mea-sure the power spectrum. In this work, we present first results using a small ACE pilot dataset and successfully demonstrate the end-to-end application of the LOFAR-EoR data processing pipeline to AARTFAAC data. Readers may refer to Yatawatta et al. (2013); Patil et al. (2017); Gehlot et al. (2018b) for an overview of the LOFAR-EoR data processing pipeline, and a description of HBA and LBA data processing.
The chapter is organised as follows: section 4.2 briefly describes the AART-FAAC wide-field imager, the observation setup of the ongoing ACE obser-vations and the current status, and the basic preprocessing steps for the raw data e.g. flagging and averaging. The calibration and imaging strategy for the ACE data is described in section 4.3. In section 4.4, we estimate and discuss the noise in Stokes I and V for the ACE data using various statistical meth-ods. In section 4.5, we describe the power spectrum estimation methodology, foreground removal process and show various power spectrum results. We summarise the work and discuss future outlook in section 4.6. We use ΛCDM cosmology throughout the analyses with cosmological parameters consistent with Planck (Planck Collaboration et al. 2016b).
4
Figure 4.1 – The left panel shows the LBA-dipole layout of the LOFAR AARTFAAC
array in its 12 station configuration (A12). The x− y axes represent the distance towards the East and the North direction, respectively. The middle panel shows the snapshot uv-coverage towards zenith at 74 MHz. The right panel shows the radial uv-density 2π|u|d|u|dN .
4.2 Observations and preprocessing
We used the LOFAR-AARTFAAC wide field imager to observe the northern sky in the frequency range 72.4− 75 MHz. The sky was observed in ’semi drift-scan’ mode13, and the observed snapshot data was processed using the
LOFAR-EoR data processing pipeline (Yatawatta et al. 2013; Patil et al. 2017; Gehlot et al. 2018a,b). The observational setup and the preprocessing steps are briefly described in following subsections.
4.2.1 The LOFAR AARTFAAC wide-field imager
The Amsterdam-ASTRON Radio Transient Facility and Analysis Centre (AART-FAAC) is a LOFAR based all-sky radio transient monitor (Prasad et al. 2016; Kuiack et al. 2018). It piggybacks on ongoing LOFAR observations and taps the digital signal streams of individual antenna elements from six or twelve core stations depending on the requirements. AARTFAAC operates in two modes viz. A6 where the six innermost stations (also called the “superterp”) of the LOFAR core are used, and A12 where twelve stations of the LOFAR core are used. The A6 mode consists of 288 dual-polarization receivers (e.g. Low Band Antenna (LBA) dipoles or High Band Antenna (HBA) tiles) within a 300 m di-ameter circle and the A12 mode consists of 576 such receivers spread across 1.2
13 AARTFAAC does not beam-form to track; it only points to the instantaneous zenith
direction per integration time in drift-scan mode. During preprocessing, the drift-scan data recorded during a long observational run are split into 15 min observation blocks. The phase centre for each block follows a constant declination and passes through the zenith mid-way during the corresponding observation block. This constitutes the ’semi drift-scan’ mode.
4
4.2: Observations and preprocessing 103
km (van Haarlem et al. 2013). Figure 4.1 shows the LBA-dipole layout, the in-stantaneous snapshot uv-coverage and the radial uv-density (dN /2π|u|d|u|). The latter is relatively flat between u = 2− 100λ except for a dip around
u∼ 20λ which is due to slightly patchy LBA-dipole layout in the “superterp”.
The inner region of the uv-plane is densely sampled and provides the same
uv-coverage at different times corresponding to a different part of the sky.
The array is co-planar at centimetre level within 0− 60λ, which is beneficial for wide-field imaging. In addition to this, the baselines up to 1.2 km support intermediate resolution imaging which helps to improve calibration and better captures compact structure in the sky. Each of the inner twelve LOFAR core stations consist of 96 LBA dipoles14(only 48 out of 96 dipoles can be used at a time) and 48 HBA tiles15. At a given time, AARTFAAC can only observe in
either LBA or HBA mode depending upon the ongoing LOFAR observation. The digitized signal from the corresponding receiver elements is tapped and transported to the AARTFAAC correlator (located at the Centre for Infor-mation Technology (CIT) in Groningen, Netherlands) prior to beam-forming. Due to network limitations, only 16 sub-bands can be correlated for 16-bit mode limited by the current network capacity. Each sub-band is 195.3 kHz wide and consists of up to 64 channels providing a maximum frequency reso-lution of 3 kHz, with currently maximum instantaneous system bandwidth of 3.1 MHz. The correlator subsystem is a GPU based correlator, which produces correlations (XX,XY,YX,YY) for all dipoles pairs for every frequency channel with 1 s integration. The correlator has 1152 input streams with 576 signal streams per polarization. The output correlations can either be dumped as raw correlations on storage disks on the AARTFAAC storage/compute clus-ter or can be routed to the AARTFAAC real-time calibration and imaging pipeline for transient detection. AARTFAAC can only observe in drift-scan mode. However, phase-tracking can be applied to raw data during or after preprocessing. The raw data from the AARTFAAC storage/compute cluster can be streamed via a fast network (1 Gbit/s) to the LOFAR-EoR processing cluster at the CIT. The raw data can be converted to standard Measurement Set (MS) format using a custom software package aartfaac2ms16 which can
also apply offline phase tracking. Readers may refer to Prasad et al. (2016) for further information about AARTFAAC system design and capabilities, and van Haarlem et al. (2013) for observing capabilities of LOFAR.
14 LBA dipoles are dual-polarization (X-Y) dipoles optimised to operate between 30-80
MHz
15 HBA tiles consist of 16 dual-polarization dipoles arranged in a 4× 4 grid, which are
analogue beam-formed to produce a single tile beam. HBA tiles are optimised to operate between 110-240 MHz. Every core station is also split into two HBA sub-stations, each consisting of 24 tiles, operating independently.
4
Table 4.1 – Observational and correlator setting details.
Parameter value
Telescope LOFAR AARTFAAC
Observation cycle and ID Cycle 10, LT10_006 Antenna configuration A12
Number of receivers 576 (LBA dipoles) Observation start time (UTC) May 25, 2018; 19:02:00 Duration of observation 6 h
Number of pointings 24 (1 pointing per 15 min) Phase centre (mean) Zenith
Minimum frequency 72.36 MHz
Maximum frequency 75.09 MHz
Target bandwidth 2.73 MHz
Outrigger sub-bands 68.36 MHz and 78.90 MHz Primary Beam FWHM 120◦ at 74 MHz
Field of View 11000 deg2 at 74 MHz
Polarisation Linear X-Y
Time, frequency resolution:
Raw Data 1 s, 65.1 kHz
After flagging and averaging 4 s, 65.1 kHz
4.2.2 ACE observational setup and status
We use the A12 mode of AARTFAAC to observe the Northern sky in the drift-scan mode with the mean phase centre at the zenith. We use 14 contiguous sub-bands (2.73 MHz bandwidth) to observe the 72.36− 75.09 MHz, targeting the redshift range z = 17.9− 18.6. We place the two remaining sub-bands
∼ 4 MHz away from the targeted band centre on either side of the band, to
aid in assessing the wide-band systematics and calibration quality as well as help foreground modelling and subtraction. We choose three channels per sub-band (with 65.1 kHz resolution) and 1 s correlator integration. High spectral and time resolution provides improved RFI excision and a better handle on frequency transform (discussed in section 4.5). So far we have recorded 125 h of driftscan observations mostly spanning the night time LSTs during May -September 2018, with a typical span of 4− 12 h per observation. From these observations, we chose an early observing run of 6 h duration (observed on May 25, 2018) for the current analysis. Table 4.1 summarises the observational and correlator setting details.
4
4.3: Calibration and Imaging 105
4.2.3 Data preprocessing
The first step of data processing is to apply tracking to the drift-scan obser-vations. For instruments with very large Field of Views (FoV; FWHM∼ 120◦ in our case), phase referencing to a single stationary point in the sky during long observations limits the portion of the sky which is visible. This is not an optimal strategy for long-duration observations. Therefore, instead of fixing the phase reference to a single stationary point for the entire observation, we choose to re-phase every 15 min observation block. The phase centre for each observation block is a constant declination point (on a great circle through zenith) which passes through zenith mid-way during the 15 min observation. We refer these re-phased 15 min observation blocks as ‘time-slices’ throughout the chapter.
The next step is RFI flagging, which is performed on the highest resolution data to minimise information loss. We use AOFlagger (Offringa et al. 2010; Offringa, van de Gronde & Roerdink 2012) to perform RFI excision on raw data and also flag all visibilities that include non-working LBA dipoles (∼ 15%). The remaining data is averaged to a resolution of 4 s and 65.1 kHz and subsequently divided into 15 min time-slices for every individual phase centre. Each time-slice is separately written into Measurement Set (MS) format, which are stored permanently on the LOFAR-EoR processing cluster. The data volume in MS format is around 150 GB for 15 min time-slices, and ∼ 3.4 TB for 6 h observations, respectively. The aartfaac2ms package performs the re-phasing and flagging tasks and returns the phased and flagged data in MS format.
4.3
Calibration and Imaging
Visibilities measured by AARTFAAC are corrupted by the errors caused due to instrumental imperfections such as complex receiver gain, primary-beam and global band-pass, as well as environmental effects, for example, due to the ionosphere. Calibration of AARTFAAC refers to the estimation of these errors and correcting the observed visibilities to obtain a reliable estimate of the true sky visibilities. The errors that corrupt the visibilities can be classified into two broad categories: Direction Independent (DI) errors and Direction Dependent (DD) errors. DI errors are independent of the direction of the incoming signal from the sky and comprise of complex receiver gain and frequency band-pass, as well as a global ionospheric phase. On the contrary, DD errors change with sky direction, e.g. as a result of the antenna voltage pattern, ionospheric phase fluctuations, and Faraday rotation (Hamaker, Bregman & Sault 1996b; Sault, Hamaker & Bregman 1996; Smirnov 2011a,b).
4
4.3.1 Direction Independent calibration
Direction Independent (DI) Calibration involves estimation of complex gains (full-Jones) per dipole (represented by a complex 2× 2 Jones matrix for two linear polarizations). We use NDPPP17to calibrate the raw visibilities and sub-sequently apply the gain solutions obtained in the calibration to the visibili-ties. Unlike SAGECal-CO (Yatawatta 2015, 2016; Yatawatta, Diblen & Spreeuw 2017; Yatawatta 2018) that we used previously in Gehlot et al. (2018b) for LBA-beam-formed data, NDPPP employs the primary beam model for individ-ual LBA dipoles. We use Cas A, Cyg A, and Vir A (the three brightest sources in the northern sky) to calibrate the visibilities. Their sky-model consists of 241 components (Delta functions and Gaussians). The models of these sources were obtained using LOFAR-LBA observations and the source fluxes in the model, within a few percent, are consistent with the Very Large Array (VLA) observations at 74 MHz (Cohen et al. 2007; Kassim et al. 2007). We use a power-law with a spectral index of −0.8 to represent the source spectra. We choose a calibration solution-time interval of 12 s for each 65.1 kHz channel to account for DI (or beam-averaged) instrumental and ionospheric effects while maintaining a reasonable signal-to-noise ratio (≳ 30) over the calibration in-terval. During calibration, we only exclude the baselines|u| < 10λ in order to avoid the large-scale diffuse Galactic emission biasing the calibration solutions. We apply a LOFAR-LBA dipole beam model18 during the model prediction
step to adjust the flux scale. Absolute flux scale can be obtained by applying the beam model before the imaging step.
4.3.2 Direction Dependent calibration
The two brightest sources Cas A and Cyg A dominate the visibilities and su-perpose significant PSF side-lobes over the field. It is crucial to subtract these sources to reduce the confusion due to these side-lobes. We use SAGECal-CO for DD-calibration and subtraction of these sources. We use a calibration so-lution interval of 24 s and 65.1 kHz respectively, with 5 ADMM iterations and a regularization parameter of ρ = 104 (Yatawatta 2016; Yatawatta, Diblen &
Spreeuw 2017; Yatawatta 2018) to enforce smoothness across frequency direc-tion. We also exclude baselines |u| < 10λ in this step for the same reason as in DI-calibration. We only use Cas A and Cyg A model in this step and use two directions corresponding to these sources. The directional gain solutions obtained are used to subtract these sources. Note that the sky-model building
17 https://www.astron.nl/lofarwiki/doku.php?id=public:user_software:
documentation:ndppp
18 Current LOFAR-LBA dipole beam models are based on Electro-Magnetic (EM)
simu-lations of the LOFAR-LBA dipoles (private communication with LOFAR Radio Obser-vatory).
4
4.3: Calibration and Imaging 107
Figure 4.2 – Intermediate resolution dirty continuum mosaics (72.39− 75.06 MHz)
combining the images of 23 DI-calibrated time-slices. The top panel shows the mo-saic in the Galactic coordinate system projected onto a Molleweide projection. The bottom panel shows the mosaic in the Equatorial coordinate system (same projection as above) with the NCP at the centre. The mosaics correspond to the LST range 11.25−17.00 h and averaged over all channels. The images were trimmed to a slightly smaller size (100◦ in pixel coordinates) to avoid projection effects near the horizon. The colour-scale is linear and is in units of Kelvins. The dotted curves represent the parallels and meridians corresponding to DEC and RA with 30◦ separation respec-tively. The solid curves represent RA=0h and DEC= 00◦. In the top panel, the grid poles on the left and right side of the mosaic correspond to the NCP and SCP, respectively.
4
for AARTFAAC is in an early stage and is not yet sufficiently complete for a DD-calibration step with more independent directions. Therefore, we only use the DD-calibrated visibilities to create a mosaic of the northern sky (discussed next). For rest of analyses in this chapter, we use the DI-calibrated visibilities where Cas A and Cyg A are still present.
4.3.3 Imaging
The visibilities after DI and DD-calibration are imaged with WSClean pack-age (Offringa et al. 2014; Offringa & Smirnov 2017). WSClean is a CPU-based wide-field interferometric imaging software that uses the w-stacking algorithm. We use a ‘Kaiser-Bessel’ kernel (Kaiser & Schafer 1980) (which is an approx-imation of the Prolate Spheroidal Wave Function (Jackson et al. 1991)) for gridding with kernel-size of 7 pixels and an oversampling factor of 63. The cell-size used in gridding is calculated from the image size after padding. We use 10− 120λ baseline range (which provide a resolution of ∼ 28′ ), with ‘natural’ weighting scheme to produce Stokes I, V and PSF images for all channels and time-slices over the full visible sky for further analysis. Another set of high-resolution Stokes I images for every channel and time-slice were produced using 0− 300λ baselines with ‘uniform’ weighting scheme, and are only used to make a mosaic map of the visible sky. Note that we do not use (multiscale) CLEAN in the imaging step. Figure 4.2 shows the latter Stokes
I continuum mosaic of the northern sky created by re-projecting the Stokes I
image for each DD-calibrated time-slice to HEALPix projection19. We observe that subtraction of the bright sources Cas A and Cyg A using DD-calibration is imperfect and still leaves significant residuals. The poor subtraction of these sources is possibly due to the incomplete sky-model used for DD-calibration.
4.4 Thermal noise statistics
To derive the thermal noise on the Stokes I and V visibilities, we use a visibility differencing method (Gehlot et al. 2018b). The DI-calibrated visibilities are divided into even and odd samplings at 4 s time resolution (whereas we used 12 s resolution in Gehlot et al. 2018b). At this time cadence, the sky and the PSF do not vary and cancel out in the difference, also, the system is expected to be stable. Thus, the difference between these even and odd time samplings provides an estimate of the thermal noise per visibility. However, because AARTFAAC has a large number of samples in the uv-plane (∼ 1.6 × 105
per channel and integration time for A12), it is computationally expensive to estimate the difference for every even and odd visibility pair. Therefore,
4
4.4: Thermal noise statistics 109
we utilised an approach where the even and odd visibility samples are first gridded using WSClean. We use the same gridding kernel and parameters as described in section 4.3.3. Note that WSClean also applies the w-correction. The obtained gridded visibility cubes are then used to calculate the difference. This approach is significantly faster than the former approach and the gridded visibilities are only a small fraction of the data volume of the raw visibilities. The azimuthally averaged spatial power spectrum (per channel) of the Stokes
I and V difference can be then written as
P∆tI(|u|, ν) ≡ ⟨∆tI˜⟩
2 = 1
2⟨(˜Ieven− ˜Iodd)
2⟩ (4.1)
where ˜Ieven and ˜Iodd are the Fourier transforms of the even and odd image
cubes respectively, u = (u, v) is the baseline vector (in units of wavelength) in the uv-plane and|u| =√u2+ v2 and ν is the frequency. For each time-slice,
we use Stokes I and V image-cubes of even and odd visibility samplings at 4 s cadence to estimate P∆tI and P∆tV respectively. We use a baseline bin-size of δ|u| = 1λ for azimuthal averaging. We also average the difference visibili-ties ∆tI corresponding to the 23 time-slices of 15 min each and estimate the˜
time-averaged noise spectra P⟨∆tI⟩ and P⟨∆tV⟩ from these averaged difference visibilities.
Figure 4.3 shows P∆t (for a single time-slice at LST 11.33 h on May 25, 2018; top row), and P⟨∆t⟩(averaged over 11.25−17.25 h LST range on the same date; bottom row) for both Stokes I and V . We observe that both Stokes I and
V noise spectra are constant along the frequency direction for baselines|u| ≲
50λ except for a few outliers possibly caused by the residual RFI. The noise estimated with the approach we employed here also depends on the number of visibilities per uv-cell, if the noise on the visibilities is uncorrelated. This effect is clearly visible along the baseline axis. The vertical structure around
|u| ∼ 20λ is due to relatively sparse uv-density compared to the neighbouring
bins, and both spectra do not vary significantly within |u| = 20 − 50λ. The power also increases for |u| > 50λ which again have sparser uv-density. This behaviour correlates well with the uv density variations shown in the right panel of figure 4.1. The power in both Stokes I and V on baselines |u| ≲ 50λ reduces by a factor of∼ 25 when averaged for 6 h, as expected for incoherent noise. However, the power in Stokes I in most cells on |u| > 50λ does not appear to average down compared to Stokes V , which is surprising. We suspect that gridding errors (or some other instrumental effect) cause this effect. The sky itself cannot cause this effect as the transition is too sharp at (∼ 50−60λ). However, the level of leakage might be proportional to the sky brightness. We are investigating this further, and leave this analysis for the future.
4
Figure 4.3 – The spatial power spectrum of the difference between gridded visibility
cubes (even - odd) at 4 s cadence. The top row shows the Stokes I and V noise spectra P∆tI(|u|, ν) (left panel) and P∆tV(|u|, ν) (right panel) for a single time-slice (at LST 11.33 h on May 25, 2018). The bottom row shows P⟨∆tI⟩and P⟨∆tV⟩for the time-slices averaged over 6 h (11.25− 17.25 h LST range).
spectra are statistically identical, assuming there is a negligible polarised com-ponent from the sky which can cause Stokes I and V noise to deviate slightly from each other (see e.g. van Straten 2009). However, the ratio P∆tI/P∆tV shown in figure 4.4 (top left panel) is significantly larger than unity and is also flat in frequency except for some outliers. The ratio between|u| = 20 − 50λ is relatively constant and has a mean and a median value of ∼ 2.70 and ∼ 2.65 respectively, which correspond to an excess of∼ 1.6 in Stokes I rms noise com-pared to Stokes V . Both the mean and the median values are consistent with each other within a few percents. Although the noise in Stokes I and V both decrease after averaging for 6 h, their ratio remains approximately the same on baselines |u| = 20 − 50λ and the distribution of these ratio values has a median value of∼ 2.4. This suggests that this excess is incoherent in time and does not average down for longer integration times. This excess is similar in behaviour to the physical excess noise we observed in Gehlot et al. (2018b) but is almost two times higher than the physical excess we observed in the beam-formed data. We also perbeam-formed a simple Kolmogorov-Smirnov (K-S) test to compare the underlying distribution of the ratios corresponding to single
4
4.5: Power Spectrum Analysis 111
Figure 4.4 – The left column shows the ratios P∆tI/P∆tV (top panel), and P⟨∆tI⟩/P⟨∆tV⟩ (bottom panel) of the spectra in figure 4.3. The right column shows normalized histograms of the distribution of the ratio values within the range |u| = 20 − 50λ. The red and blue vertical lines represent the mean and the median values of the corresponding distribution.
time-slice and 6 h average. The p-value obtained from the K-S test is∼ 10−19, suggesting that the two distributions are unlikely to be drawn from the same underlying probability distribution. We also observe that Stokes I noise for baselines > 50λ does not average down in time compared to the Stokes V noise. We are currently investigating this but we suspect gridding errors, and flagging differences might cause this effect (Offringa et al. submitted). The effect of visibility migration from one uv-cell to another is negligible at the 4 s level, since, typical time-scale for baseline migration are of order ∼ 2 min at|u| = 60λ. As the noise in Stokes I and V is lowest and their ratio is flat within |u| = 20 − 50λ range, we only analyse this baseline range throughout the power spectra and foreground removal in the remainder of this chapter.
4.5
Power Spectrum Analysis
Cylindrically averaged power spectra in (k⊥, k∥) space are the most commonly used statistical tool to study the challenges associated with foreground contam-ination and systematic biases (Bowman, Morales & Hewitt 2009; Vedantham, Udaya Shankar & Subrahmanyan 2012). The wave mode k⊥ represents the scale of the brightness temperature fluctuations in the plane perpendicular to the line of sight and the wave mode k∥ represents the scale of the fluctuations
4
along the line of sight. Foregrounds, the ionospheric effects and systematic biases which are smooth in frequency reside within a region often called “the wedge”. The cylindrically averaged power spectrum P (k⊥, k∥) is defined as (Parsons et al. 2012; Thyagarajan et al. 2015a):
P (k⊥, k∥) = X
2Y
ΩAB⟨|˜V(u, τ)|
2⟩, (4.2)
where ˜V(u, τ) is the Fourier transform of the visibilities along the frequency direction, B is the bandwidth of the visibility cube and ΩAis the primary beam solid angle. X and Y are the cosmological conversion factors to convert angle and frequency to transverse co-moving distance (D(z)) and the co-moving depth along the line of sight (∆D), respectively. The wave numbers k⊥ and
k∥ are related to the baseline vector (u = (u, v) in units of wavelength) and the delay (τ ) as (Morales & Hewitt 2004):
k⊥ = 2π (|u|)
D(z) , k∥ =
2πν21H0E(z)
c(1 + z)2 τ , (4.3)
where ν21 is the rest frame frequency of the 21-cm spin flip transition of
neu-tral hydrogen, z is the redshift corresponding to the observation frequency,
c is the speed of light, H0 is the Hubble constant at z = 0 and E(z) ≡
[
ΩM(1 + z)3+ Ωk(1 + z)2+ ΩΛ
]1/2
is a function of the standard cosmological parameters . The spherically averaged dimensionless power spectrum (∆2(k))
can be determined from P (k⊥, k∥) as
∆2(k) = k
3
2π2P (k), (4.4)
where k =√k2
⊥+ k∥2 and P (k) is the spherically averaged power spectrum.
4.5.1 Cylindrically averaged power spectrum estimation
One approach is to estimate the power spectrum from the gridded visibilities obtained after Fourier transforming the images along the spatial and frequency axes. This approach is utilised in LOFAR power spectra analyses where phase tracking is used (Patil et al. 2016, 2017; Gehlot et al. 2018b). In the case of AARTFAAC, however, tracking is applied only for short time periods, and each time-slice needs to be combined later in a coherent manner to improve the signal-to-noise ratio before power spectrum estimation. Moreover, AART-FAAC is a wide-field imager and imaging large FoVs for long time integrations with traditional gridding algorithms is non-trivial and produces sub-optimal results. Drift-scan instruments such as PAPER and HERA use ‘fringe rate
4
4.5: Power Spectrum Analysis 113
filter’ (Parsons & Backer 2009; Parsons et al. 2016), which is similar to grid-ding in uv-plane, to coherently average the visibilities for longer integrations. However, this approach is only useful for per-baseline power spectrum mea-surement and this methodology leads to loss of the cosmological signal (Cheng et al. 2018). Furthermore, it is computationally expensive for large AART-FAAC datasets. Another popular method to produce images with coherent sky averaging is m-mode analysis (Shaw et al. 2014, 2015; Eastwood et al. 2018), which is used to make sky maps using drift-scan observations. This approach is based on spherical harmonics and takes advantage of the fact that sky re-peats after every 23h56m (a sidereal day). This approach has been successfully used to create sky-maps with the Long Wavelength Array (see Eastwood et al. 2018 for more details). However, an m-mode analysis inversion is not possible with current AARTFAAC system as it requires full sidereal day observations. Because of limited storage space on the AARTFAAC storage/compute cluster, such observations with required resolution settings (see section 4.2) cannot be accommodated on the cluster. Also, the daytime observations are unfavourable at LBA frequencies due to RFI and erratic solar/ionospheric effects.
We are currently developing a hybrid strategy to combine the “semi drift-scan” time-slices in a coherent manner to estimate the power spectrum of the aver-aged sky. In this approach, the gridded visibilities for each time-slice are trans-formed to spherical harmonics basis (represented with spherical Fourier-Bessel basis functions) using Spherical Harmonics Maximum Likelihood (SpH-ML) in-version technique (see Ghosh, Mertens & Koopmans 2018). The visibilities in spherical harmonics basis can be combined coherently by rotating the fitted spherical harmonics modes for each time-slice to a common reference. The coefficients of the fitted spherical harmonics basis functions provide a direct estimate of the angular power spectrum Cℓ (which can be converted to a
cosmological power spectrum after proper scaling with cosmological parame-ters). Moreover, these spherical harmonics modes provide full-sky maps for modelling large-scale diffuse emission in mitigating the foreground contamina-tion with Gaussian Process Regression (GPR; Mertens, Ghosh & Koopmans 2018). Although, this method is reasonably fast for small time-slices and fur-ther speed up can be achieved by using a limited number of ℓ-modes in the inversion, it is computationally much more expensive than w-projection algo-rithm and m-mode analysis. Therefore, we currently determine power spectra from the gridded visibilities for every time-slice after w-correction, ignoring sky curvature.
Cylindrical power spectra
We use the DI-calibrated image cubes to first determine the spatial power spectra P (k⊥, ν) for every frequency channel. We do not use DD-calibrated
4
Figure 4.5 – The spatial power spectra per frequency channel P (k⊥, ν) for Stokes I (left column) and V (middle column) DI-calibrated gridded visibilities before fore-ground removal for|u| = 20 − 50λ. The right column shows corresponding full-sky Stokes I images (but produced with|u| = 0 − 120λ range). The black squares show the 50◦× 50◦ region used to determine the power spectra. A smaller region is chosen to avoid the Milky Way, however, it also reduces the power spectrum sensitivity as a large portion of the sky is removed. Different rows correspond to different LSTs in the range 11.25 - 17.25 h separated by 1 h each, with the top row being the lowest LST. The ripples at∼ 200 kHz level in Stokes V spectra for the last two rows are due to edge-effects of the poly-phase filter.
4
4.5: Power Spectrum Analysis 115
Figure 4.6 – The cylindrically averaged power spectra P (k⊥, k∥) for Stokes I (left column) and V (right column) before foreground removal, corresponding to the time-slices shown in figure 4.5. The right column shows their ratio PI(k⊥, k∥)/PV(k⊥, k∥). The horizontal lines at k∥∼ 2.0 in Stokes I and V correspond to the 200 kHz ripple due to the poly-phase filter. The gray lines correspond to the instrumental horizon delay as a function of baseline.
4
image cubes because the sky-model building for AARTFAAC is in an early stage and is not yet sufficiently complete for a reliable DD-calibration step. The DD-calibration can distort the desired 21-cm signal if the sky-model is still incomplete (Mouri Sardarabadi & Koopmans 2018). We are working on improving the sky-model to include more sources in future analyses. The image cubes used here are produced using ‘natural’ weighting scheme and the gridding parameters described in section 4.3.3. The inner 50◦ × 50◦ region of these DI-calibrated image cubes is used for the power spectrum estimation (see e.g. Offringa et al. submitted) to limit primary beam effects and avoid the Milky Way. However, it also leads to a relatively lower power spectrum sensitivity. We use the baselines|u| = 20 − 50λ and a bin size of δu = 1.15λ (corresponding to 50◦ FoV) for the power spectrum estimation and azimuthal averaging.
Figure 4.5 shows P (k⊥, ν) for DI-calibrated Stokes I (left column) and V
(middle column) visibilities for different time-slices (different rows), within the LST range 11.25− 17.25 h. The first time-slice (first row) corresponds to the LST around the time of sunset (11.25h), whereas the other four time-slices were observed after sunset. The emission in Stokes I is dominated by smooth foregrounds, and as the time passes, the overall power in Stokes I increases. This is predominantly caused by Cas A and Cyg A, which are near the horizon at the beginning of the observation and move to higher elevations during the observation. Moreover, at the same time, the Galactic plane also rises in the sky. A combination of these bright sources dominate the foreground emission and leads to the increased power in Stokes I. The Stokes V power spectrum for the first time-slice have noisier channels compared to the time-slices in the second and third row, suggesting that the observations when the Sun is above the horizon might have worse noise level due to the more erratic ionosphere or the Sun itself. It is also possible that the RFI environment becomes more benign during the night time. However, the time-slices during LST 15.0h to 16.5h (last two rows in figure 4.5) show a prominent ripple of∼ 200 kHz width across frequency direction. This ripple is caused by the frequency band-pass structure at the sub-band level due to the edge-effects of the poly-phase filter and is present in all time-slices. It is amplified when Cas A, Cyg A and the galactic plane are near/inside the analysis region possibly due to polarization leakage from Stokes I to V .
To obtain the cylindrically averaged power spectra, P (k⊥, k∥), the visibili-ties are Fourier transformed along frequency direction. Before the frequency transform, we flag some channels affected by residual RFI to avoid any ar-tifacts in the power spectrum. We use a Discrete Fourier Transform (DFT) along frequency. We apply a ‘Blackman’ spectral window20 to the data prior 20 The ’Blackman’ window is defined as:
4
4.5: Power Spectrum Analysis 117
to the DFT to reduce the side-lobe leakage along the delay axis. Figure 4.6 shows the cylindrically averaged Stokes I and V power spectra, PI(k⊥, k∥) and
PV(k⊥, k∥) respectively, for the same time-slices as shown in figure 4.5. We
observe that the Stokes I and V power spectra for the first time-slice have relatively high power compared to later times, as a result of the noisier chan-nels we observed in figure 4.5. The “wedge” is also clearly visible in Stokes
I, and it extends beyond the horizon (which lies below k∥ ∼ 0.5 h cMpc−1). Moreover, Stokes V also shows a fainter “wedge” on the smallest k∥ modes, which might be the bright sources e.g. Cas A, Cyg A, Galactic Centre etc. leaking from Stokes I to V due to instrumental polarization leakage. Gehlot et al. 2018a observed strong polarization leakage in LBA due to Cas A. We also observe an artefact of the 200 kHz band-pass ripple due to the poly-phase filter in Stokes I at k∥∼ 2.0 in the time-slices after sunset, whereas Stokes V starts to show the ripple artefact when bright sources are at higher elevations in the sky. The ratio PI(k⊥, k∥)/PV(k⊥, k∥), above the “wedge” is higher for
the time-slice before sunset and becomes lower at later times. This suggests that the night-time observations are preferable over day-time observations as the power in Stokes I and V outside the “wedge” are comparable during night time observations.
4.5.2 Foreground Subtraction with GPR
Subtraction of the bright foreground emission a crucial step in 21-cm signal experiments. The intrinsic foreground emission has two dominant components being diffuse Galactic synchrotron and thermal emission, and extra-galactic sources (e.g. radio galaxies, supernova remnants etc.) (Di Matteo et al. 2002; Zaldarriaga, Furlanetto & Hernquist 2004; Bernardi et al. 2009; Ghosh et al. 2012). In addition to this, the instrument also imparts a spectral structure on the data called instrumental mode-mixing due to its frequency response (Datta, Bowman & Carilli 2010; Morales et al. 2012; Trott, Wayth & Tingay 2012; Vedantham, Udaya Shankar & Subrahmanyan 2012; Hazelton, Morales & Sullivan 2013). On the other hand, the 21-cm signal varies rapidly with fre-quency. Gaussian Process Regression (GPR) (Rasmussen & Williams 2005) exploits this distinct spectral behaviour of the intrinsic foregrounds, instru-mental mode-mixing, and the 21-cm signal to separate them from each other. GPR models these different components with Gaussian Processes (GPs), using different covariance functions representing the spectral correlation functions of the different components. Readers may refer to Mertens, Ghosh & Koopmans (2018) for an overview of GPR and its application for foreground removal and
W (n) = 0.42− 0.5 cos [ 2πn (M− 1) ] + 0.08 cos [ 4πn (M− 1) ] , where 0≤ n ≤ M − 1 (see e.g. Blackman & Tukey 1958).
4
Figure 4.7 – Cylindrically averaged Stokes I (left column) and V (middle column)
power spectra (P (k⊥, k∥)) after foreground removal, corresponding to the time-slices shown in figure 4.5. The right column shows their ratio PI(k⊥, k∥)/PV(k⊥, k∥).
4
4.5: Power Spectrum Analysis 119
signal separation. The GPR has been demonstrated to work successfully on both LOFAR-HBA and LBA data in Mertens, Ghosh & Koopmans (2018) and Gehlot et al. (2018b) respectively.
In this pilot analysis, we use DI-calibrated visibilities for foreground removal. We select the 42 channels of 65.1 kHz (totalling 2.73 MHz bandwidth) to perform the foreground removal. It is difficult to differentiate between the smooth foregrounds emission and the mode-mixing component if the band-width is only several MHz. On these short baselines, the wedge only occupies a small part of the power spectrum space (see figure 4.6). Therefore, we use a single RBF (Radial Basis Function) covariance function with a coherence scale with a prior of 0.8− 10 MHz to model the foreground component (both in-trinsic and mode-mixing). The 21-cm signal is modelled using a Matern class kernel with a coherence scale with a prior of 0.05− 1.5 MHz. The coherence scales are optimised by maximising the Bayesian evidence (Mertens, Ghosh & Koopmans 2018). We use the Stokes V noise, estimated in section 4.4, as a prior for the noise covariance in GPR. To remove the 200 kHz ripple from Stokes I and V visibilities, we use Principal Component Analysis (PCA) in conjunction with GPR. First, GPR is used to remove the foregrounds. Next, PCA is run on the residuals, and the first principal component is subtracted from the original visibilities before the foregrounds were removed. Finally, GPR is performed on the residuals after subtraction of the first principal com-ponent to remove the remaining foreground emission. GPR is run separately on each time-slice since the foreground and mode-mixing components change as a function of time.
Figure 4.7 shows the cylindrically averaged Stokes I and V power spectra after foreground and band-pass structure removal. We observe that the foreground emission that originally dominated the lowest k∥ modes in Stokes I and V power spectra (see figure 4.6), as well as the 200 kHz ripple, are effectively removed. Both Stokes I and V appear noise-like after foreground removal. Their ratio is approximately constant over most k⊥ and k∥ modes with a median value∼ 3 (except the last time-slice), consistent with the ratio of noise spectra estimates in section 4.4. The behaviour of this excess remains puzzling because it decorrelates on a 4 s time-scale and does not show any frequency coherence. We also observed a similar but lower excess in the beam-formed LBA data (Gehlot et al. 2018b). We are currently investigating this effect in more detail and we leave this analysis for the future.
4.5.3 Spherically averaged power spectrum estimation
We determine the spherically averaged dimensionless Stokes I and V power spectra (∆2(k)) from the residual visibilities after foreground removal for each
4
Figure 4.8 – The spherically averaged Stokes I (blue curves) and V (red curves)power spectra for each time-slice before (left panel) and after (right panel) foregroundremoval. Different shades of colour represent different time-slices with the darkest shade being the first and lightest shade being the last time-slice.
time-slice. Figure 4.8 shows ∆2
I(k) and ∆2V(k) before and after foreground
removal for each time-slice. The power in Stokes I on the smallest k modes is dominated by diffuse foregrounds and it decreases by more than 3 orders of magnitude after foreground removal. The sky rotation also affects the shape of the power spectrum. We observe that the power increases on small k modes with time as the Galactic plane, Cas A and Cyg A move to higher elevations. The power in Stokes V is presumably dominated by polarization leakage from I to V , which is removed by GPR because it is spectrally smooth. Both ∆2
I(k)
and ∆2
V(k) after foreground removal appear power-law like for most
time-slices, which is expected if P (k) is approximately constant on all k modes. Surprisingly, the power in residual Stokes I and V is lower for the time-slices where Cas A and Cyg A are at higher elevations. This could be due to DI-calibration being less erroneous when the calibrator sources are closer to zenith thus resulting in lower noise in the images.
We also compare the residual power spectra with the power spectra of Stokes
I and V noise determined from visibility differencing in section 4.4. In order
to have a better comparison and visualisation, we average the residual Stokes
I, V and their noise spectra using inverse variance weighting (which is an
optimal method for averaging the power spectra), where we use the squared error on each power spectra as corresponding weights. The averaged power spectra⟨∆2
I⟩t,⟨∆2V⟩t,⟨∆2I,n⟩tand ⟨∆2V,n⟩tare shown in figure 4.9. We observe
4
4.5: Power Spectrum Analysis 121
Figure 4.9 – Optimally time averaged Stokes I and V power spectra ⟨∆2
I⟩t and ⟨∆2
V⟩t (solid blue and red lines respectively), and their corresponding noise spectra ⟨∆2
I,n⟩tand⟨∆2V,n⟩t(dashed grey and black lines respectively).
Stokes V , which is similar to the ratio that we observed before. Similarly, the Stokes I noise spectrum is also higher than the Stokes V noise by a factor of
∼ 2. Also, the Stokes V power spectrum after foreground removal is consistent
with Stokes I noise power spectrum. Noise bias correction
The excess that we observe, needs to be accounted for while removing the noise bias from the residual Stokes I power spectra. Since the origin of the excess power is not clear, we employ two different strategies to remove this noise bias from Stokes I and estimate noise bias corrected Stokes I power spectrum ∆2I−∆2N: (1) We scale the Stokes V power spectrum (after foreground removal) with the median of the ratio of Stokes I and V noise estimates for each time-slice individually. The scaled Stokes V is then used to remove the bias in Stokes
I; (2) We assume that ∆2I,n represents the noise bias in Stokes I, and use it to correct for the noise bias in Stokes I for each time-slice. Each strategy has its advantages. Whereas the first strategy assumes that potential gain errors might amplify incoherent noise represented by Stokes V , the second strategy assumes that an unknown excess power exists (as shown in e.g. Gehlot et al.
4
Figure 4.10 – Noise bias corrected power spectra ∆2
I− ∆2N for two strategies. The dashed grey line shows the Stokes V noise power spectrum ∆2
V,n. The errorbars represent the 2σ errors on the power spectra.
2018b and in section 4.4) in Stokes I over Stokes V . However, both strategies neglect the effect of calibration errors leaking power to shorter baselines. ∆2
I−
∆2N obtained for the two strategies are averaged incoherently in time using the inverse variance weighting scheme. We are currently investigating more optimal strategies to coherently combine different time-slices to obtain a single estimate of Stokes I and V power spectra, but for now, we adhere to the incoherent averaging.
Figure 4.10 shows the noise bias corrected power spectra ∆2
I− ∆2N for the two
strategies. We observe that both power spectra are still higher than the noise power spectrum by factors of few. However, the errors are much lower for both cases. With strategy (1), a 2σ upper limit of ∆2
21< (71 K)2 is reached at
k = 0.145 h cMpc−1 with the 1σ error of ∼ (7.9 K)2. With strategy (2), a 2σ
upper limit of ∆2
21 < (86.7 K)2 is reached at k = 0.145 h cMpc−1 with the 1σ
error of∼ (7.6 K)2. Although, the upper limits that still include a considerable
bias are currently not so interesting, the error on the power spectrum after averaging is quite small. We plan to incorporate coherent averaging of the residual visibilities soon, which will improve these upper limits by a large factor.
4
4.6: Summary and future work 123
4.6
Summary and future work
In this work, we described the AARTFAAC Cosmic Explorer (ACE) program motivated by the recently reported detection of the deep absorption feature in sky averaged spectrum of the 21-cm signal during Cosmic Dawn by the EDGES collaboration (Bowman et al. 2018). Main results of the chapter are summarised below:
1. We demonstrate the successful end-to-end application of the LOFAR-EoR data processing pipeline to ACE data, starting from preprocessing and calibration to power spectrum estimation after foreground removal. 2. We observe that the ratio of noise estimates in Stokes I and V determined by differencing even and odd visibility samples with a cadence of 4 s and averaged over 15 min has a median value of ∼ 2.8. This excess remains approximately constant when averaged over 6 h integration, even though the noise decreases in both Stokes I and V . This excess has a peculiar behaviour and is of unknown origin. A similar but lower excess has been found in beam-formed LBA data (Gehlot et al. 2018b). We are currently investigating this in detail and leave the analysis for the future.
3. The Stokes V power during the observation period before and around the time sunset is higher than the power after the sunset in several channels. It is possibly caused by ionospheric activity around sunset or a more benign RFI environment after sunset.
4. With noise bias removal strategy (1), 2σ upper limit of ∆2
21 < (71 K)2
is reached at k = 0.145 h cMpc−1 with the 1σ error of∼ (7.9 K)2. With
strategy (2), we reach a 2σ upper limit of ∆2
21 < (86.7 K)2 at k =
0.145 h cMpc−1 with the 1σ error of ∼ (7.6 K)2 in the redshift range
z = 17.9− 18.7.
4.6.1 Future outlook and forecast
Although we demonstrate the successful application of the LOFAR-EoR data processing pipeline on ACE data, most of the steps used in the analysis are fairly rudimentary. In the future, we plan to improve the processing and analysis by improving several aspects such as:
1. Improve DI and add DD calibration by including a detailed sky-model and enforce frequency smoothness of the gain solutions as a constraint in the calibration. A further improvement in DD-calibration can still be achieved by increasing ADMM iteration and by better initialisation of gain solutions.
4
2. Improve the imaging strategy and develop optimal methods to coherently combine the visibilities in time to achieve a lower noise floor.
3. Study the effect of polarization leakage and the ionosphere in ACE ob-servations and mitigate them if required.
4. Investigate the origin of the excess noise in Stokes I compared to Stokes
V and avoid this to determine a more accurate estimator of thermal
noise.
The above-mentioned improvements and avoidance of excess noise will allow us to remove the noise bias properly, leading to a power spectrum sensitivity of ∆221< (0.5 K)2. This is still above the thermal noise level but is already below several predictions by Fialkov, Barkana & Cohen (2018). Further improve-ments in calibration and foreground removal will improve the power spectrum sensitivity.
