Science, Faculty of
Physics and Astronomy, Department of
DSpace
UBCV
Pinsonneault-Marotte, Tristan
Doctor of Philosophy - PhD
University of British Columbia
The Canadian Hydrogen Intensity Mapping Experiment (CHIME) is a radio telescope that we built to map the large-scale structure of the Universe between redshifts 0.8 and 2.5, when dark energy is expected to begin the transition from a decelerating to an accelerating phase in the expansion of the Universe. CHIME was designed to perform an intensity mapping survey using the 21 cm line of neutral hydrogen, a novel method that has the potential to enable enormous surveys of the distant Universe, but also significant observational challenges to overcome. In this thesis, I describe contributions I made to the CHIME data acquisition system and calibration effort, culminating in a detection of cosmological 21 cm emission in cross-correlation with measurements of the Lyman-𝜶 forest. The large data rate from the CHIME correlator is processed in real time by a high-performance digital pipeline, the development of which I participated in extensively. A few specific processing tasks where I led the design and implementation are highlighted in this work. In order to detect the 21 cm signal amidst the much brighter foreground emission from nearby sources, a very precise instrumental calibration is required. Calibrating the telescope's beam is a particular concern. One of many approaches being pursued for CHIME is the holographic observation of bright celestial sources in concert with a second radio telescope. I describe work I did to derive beam measurements from such observations and their analysis, including a scheme for calibrating the polarised beam response. I report the detection of 21 cm emission at an average redshift z = 2.3 in the cross-correlation of CHIME maps with measurements of the Lyman-𝜶 forest from the eBOSS. Data collected by CHIME over 88 days in the 400-500 MHz frequency band (1.8 < z < 2.5) were formed into maps of the sky and high-pass delay filtered to suppress the foreground power. Line-of-sight spectra to the eBOSS background quasar locations were extracted from the CHIME maps and combined with the Lyman-𝜶 forest flux transmission spectra to estimate the 21 cm-Lyman-𝜶 cross-correlation function. Fitting a simulation-derived template to this measurement results in a detection of 9-𝝈 significance.
A Detection of Cosmological 21 cm Emission fromCHIME in Cross-correlation with the eBOSSLyman-α ForestbyTristan Pinsonneault-MarotteMSc Physics, University of British Columbia, 2018BSc Honours Physics, McGill University, 2016a thesis submitted in partial fulfillmentof the requirements for the degree ofDoctor of Philosophyinthe faculty of graduate and postdoctoralstudies(Physics)The University of British Columbia(Vancouver)October 2023© Tristan Pinsonneault-Marotte, 2023The following individuals certify that they have read, and recommend to the Facultyof Graduate and Postdoctoral Studies for acceptance, the thesis entitled:A Detection of Cosmological 21 cm Emission from CHIME inCross-correlation with the eBOSS Lyman-α Forestsubmitted by Tristan Pinsonneault-Marotte in partial fulfillment of the require-ments for the degree of Doctor of Philosophy in Physics.Examining Committee:Gary Hinshaw, Professor, Physics and Astronomy, UBCSupervisorMark Halpern, Professor, Physics and Astronomy, UBCSupervisory Committee MemberDouglas Scott, Professor, Physics and Astronomy, UBCUniversity ExaminerEd Grant, Professor, Chemistry, UBCUniversity ExaminerAdditional Supervisory Committee Members:Scott Oser, Professor, Physics and Astronomy, UBCSupervisory Committee MemberDavid Morrissey, Adjunct Assistant Professor, Physics and Astronomy, University ofVictoria; Research Scientist, TRIUMFSupervisory Committee MemberiiAbstractThe Canadian Hydrogen Intensity Mapping Experiment (CHIME) is a radio telescopethat we built to map the large-scale structure of the Universe between redshifts 0.8and 2.5, when dark energy is expected to begin the transition from a deceleratingto an accelerating phase in the expansion of the Universe. CHIME was designed toperform an intensity mapping survey using the 21 cm line of neutral hydrogen, a novelmethod that has the potential to enable enormous surveys of the distant Universe,but also significant observational challenges to overcome. In this thesis, I describecontributions I made to the CHIME data acquisition system and calibration effort,culminating in a detection of cosmological 21 cm emission in cross-correlation withmeasurements of the Lyman-α forest.The large data rate from the CHIME correlator is processed in real time by ahigh-performance digital pipeline, the development of which I participated in exten-sively. A few specific processing tasks where I led the design and implementation arehighlighted in this work.In order to detect the 21 cm signal amidst the much brighter foreground emissionfrom nearby sources, a very precise instrumental calibration is required. Calibratingthe telescope’s beam is a particular concern. One of many approaches being pursuedfor CHIME is the holographic observation of bright celestial sources in concert witha second radio telescope. I describe work I did to derive beam measurements fromsuch observations and their analysis, including a scheme for calibrating the polarisedbeam response.I report the detection of 21 cm emission at an average redshift z¯ = 2.3 in thecross-correlation of CHIME maps with measurements of the Lyman-α forest from theeBOSS. Data collected by CHIME over 88 days in the 400-500 MHz frequency band(1.8 < z < 2.5) were formed into maps of the sky and high-pass delay filtered to sup-press the foreground power. Line-of-sight spectra to the eBOSS background quasarlocations were extracted from the CHIME maps and combined with the Lyman-α for-est flux transmission spectra to estimate the 21 cm-Lyman-α cross-correlation func-iiition. Fitting a simulation-derived template to this measurement results in a detectionof 9σ significance.ivLay SummaryWe have known for about a century that the Universe is expanding, but more recentlyit was discovered that the expansion is accelerating. To explain this within ourunderstanding of gravity, a previously unsuspected type of energy is hypothesizedto exist and make up ∼ 70% of the energy in the Universe today. This gap in ourunderstanding of fundamental physics is referred to as dark energy.The Canadian Hydrogen Intensity Mapping Experiment is a large radio telescopethat we built to map the distribution of matter in the distant reaches of the Universeby measuring light emitted by hydrogen, and thus track the expansion to learn aboutdark energy. The design of CHIME is novel, and making it work will require over-coming many challenges. In this thesis I describe my contributions to building thesystem that records data from CHIME, efforts towards its calibration, and finally apreliminary detection of hydrogen located nearly 11 billion light-years away.vPrefaceThe Canadian Hydrogen Intensity Mapping Experiment (CHIME) was designed andbuilt by a team of researchers spanning many institutions over many years, andmuch of the work that I performed thus was in collaboration with others or relies onprevious contributions. In particular, the analysis of CHIME data is implementedin Python code that was collaboratively developed and is reused extensively withinthe collaboration. A fraction of this code is publicly released (Shaw et al., 2020c,d,a;Renard et al., 2021), where my individual contributions are explicitly recorded in thecode versioning.An overview of the experiment was published in CHIME Collaboration et al.(2022a) (“the overview paper” in what follows), of which I am a co-author, and Ihave reproduced some text and figures from that paper when describing CHIMEin this work. The first detection of the 21 cm signal by CHIME was published inCHIME Collaboration et al. (2022b) (“the stacking paper” hereafter), where I amalso a co-author. The analysis that led to that result includes a number of steps thatare shared with the analysis I report in Chapter 4, and where appropriate I haveincluded and adapted material from that paper to describe those steps. Chapter 4is the subject of a paper currently in preparation that will soon be submitted forpublication. Instances where material from these published works was reused in thisthesis are listed in the remainder of the preface. Otherwise, all text and figures, aswell as the work they report, are from me.Figure 1.2 was assembled from the publications Planck Collaboration et al. (2014)and Planck Collaboration et al. (2020), and Figure 1.3 was reproduced from publicmaterials provided by eBOSS. The introduction to hydrogen intensity mapping andCHIME in sections 1.1.4 and 1.3 uses some snippets of text from the overview paper.Figures 1.5, 1.6 and 1.7 and their captions are reproduced entirely from that paper.Figure 2.1 was adapted from a figure published in the overview paper. Figure 3.1is a photograph taken by a member of the CHIME collaboration and made availablefor public use. A sentence and references describing holography in Section 3.1 werevireused from the overview paper. Figure 3.2 and its caption are reproduced fromthe overview paper. Figure 3.3 was generated by Alex Reda for a paper that is inpreparation.Figure 4.3 was reproduced from the overview paper, and its caption adaptedfrom there. The discussion in Section 4.4.3 was adapted from the stacking paper.Figure 4.5 was reproduced from the stacking paper. Portions of Section 4.5.2 wereadapted from the stacking paper, as were Sections 4.5.3 and 4.5.4. The descriptionof the cosmological and CHIME signal simulations in Sections 4.6.1 and 4.6.2 areadapted from the stacking paper. The analysis and code that are described in thesections mentioned so far in this paragraph are the work of other members of thecollaboration, which I have used and adapted for my analysis. The remainder of theanalysis and code implementation presented in Chapter 4 are my work. Figure 4.13was created by me based on a similar figure published in the stacking paper.viiTable of ContentsAbstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iiiLay Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vPreface . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viTable of Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viiiList of Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xiList of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xiiList of Abbreviations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xvAcknowledgements . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xvii1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11.1 Cosmology . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21.1.1 The Geometry of Spacetime . . . . . . . . . . . . . . . . . . . 21.1.2 The Accelerating Universe and ΛCDM . . . . . . . . . . . . . 61.1.3 BAO Standard Ruler . . . . . . . . . . . . . . . . . . . . . . . 111.1.4 Hydrogen Intensity Mapping . . . . . . . . . . . . . . . . . . 131.2 Radio Interferometers . . . . . . . . . . . . . . . . . . . . . . . . . . 171.2.1 Beyond the Aperture . . . . . . . . . . . . . . . . . . . . . . . 171.2.2 Visibilities . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 201.3 CHIME . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 211.3.1 Instrument . . . . . . . . . . . . . . . . . . . . . . . . . . . . 241.3.2 Real-time Calibration and Flagging . . . . . . . . . . . . . . . 292 Receiver Pipeline Development . . . . . . . . . . . . . . . . . . . . . 32viii2.1 Final Data Compression Stage . . . . . . . . . . . . . . . . . . . . . . 332.1.1 Floating Point Truncation . . . . . . . . . . . . . . . . . . . . 332.1.2 Software Implementation . . . . . . . . . . . . . . . . . . . . . 372.1.3 Compression Performance . . . . . . . . . . . . . . . . . . . . 372.2 Pulsar Gating . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 382.2.1 Gating Scheme . . . . . . . . . . . . . . . . . . . . . . . . . . 382.2.2 Timing Model . . . . . . . . . . . . . . . . . . . . . . . . . . . 402.2.3 Software Implementation . . . . . . . . . . . . . . . . . . . . . 412.3 Real-time Maps . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 412.4 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 433 Beam Calibration . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 463.1 Holography Data Processing . . . . . . . . . . . . . . . . . . . . . . . 463.1.1 Fringestopping and Regridding . . . . . . . . . . . . . . . . . 503.1.2 Decorrelation Correction . . . . . . . . . . . . . . . . . . . . . 503.1.3 Data Packing . . . . . . . . . . . . . . . . . . . . . . . . . . . 513.2 Analysis of Unexpected Slopes in Holography Amplitudes . . . . . . . 523.2.1 Characterisation of the Slopes . . . . . . . . . . . . . . . . . . 523.2.2 Expected Effects . . . . . . . . . . . . . . . . . . . . . . . . . 563.2.3 Other Checks . . . . . . . . . . . . . . . . . . . . . . . . . . . 613.3 Pulsar Holography . . . . . . . . . . . . . . . . . . . . . . . . . . . . 653.3.1 Disentangling the CHIME and Galt Gains . . . . . . . . . . . 663.3.2 Separate Gain Fit . . . . . . . . . . . . . . . . . . . . . . . . 693.3.3 Including Polarisation . . . . . . . . . . . . . . . . . . . . . . 703.3.4 Calibrating the Galt Jones Matrix . . . . . . . . . . . . . . . 713.4 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 744 Cross-correlation of CHIME and Lyman-α forest Measurements . 764.1 Motivation/Background . . . . . . . . . . . . . . . . . . . . . . . . . 774.2 Data . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 774.2.1 eBOSS . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 774.2.2 CHIME . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 794.3 Cross-correlation Method . . . . . . . . . . . . . . . . . . . . . . . . 814.3.1 Frequency Remapping . . . . . . . . . . . . . . . . . . . . . . 824.3.2 Correlation Function Estimation . . . . . . . . . . . . . . . . 834.4 Foreground Delay Filter . . . . . . . . . . . . . . . . . . . . . . . . . 834.4.1 Delay Power Spectra . . . . . . . . . . . . . . . . . . . . . . . 84ix4.4.2 Foreground Wedge . . . . . . . . . . . . . . . . . . . . . . . . 864.4.3 DAYENU Filter . . . . . . . . . . . . . . . . . . . . . . . . . . . 874.5 Data Processing Pipeline . . . . . . . . . . . . . . . . . . . . . . . . . 894.5.1 Daily Processing . . . . . . . . . . . . . . . . . . . . . . . . . 904.5.2 Sidereal Average . . . . . . . . . . . . . . . . . . . . . . . . . 934.5.3 Point-Source Subtraction . . . . . . . . . . . . . . . . . . . . 944.5.4 Map-Making . . . . . . . . . . . . . . . . . . . . . . . . . . . 944.5.5 Foreground Filter . . . . . . . . . . . . . . . . . . . . . . . . . 1004.5.6 Extraction of Line-of-sight Spectra . . . . . . . . . . . . . . . 1024.5.7 Correlation Function . . . . . . . . . . . . . . . . . . . . . . . 1054.6 Signal Model and Simulations . . . . . . . . . . . . . . . . . . . . . . 1054.6.1 Large-scale Structure Simulation . . . . . . . . . . . . . . . . 1054.6.2 CHIME Maps Simulation . . . . . . . . . . . . . . . . . . . . 1064.6.3 Lyman-α Forest Simulation . . . . . . . . . . . . . . . . . . . 1104.6.4 Cross-correlation Template . . . . . . . . . . . . . . . . . . . 1124.7 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1144.7.1 Detection Significance . . . . . . . . . . . . . . . . . . . . . . 1154.7.2 Comparison to fast fourier transform (FFT) Cross-correlationMethod . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1164.8 Validation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1184.8.1 Lyman-α Forest Permutations . . . . . . . . . . . . . . . . . . 1184.8.2 Even/Odd Days . . . . . . . . . . . . . . . . . . . . . . . . . 1184.9 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1225 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 124Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 127A Reverse Lanczos Interpolation . . . . . . . . . . . . . . . . . . . . . 133B Spectral Incoherence . . . . . . . . . . . . . . . . . . . . . . . . . . . 135C Spatial Aliasing in Interferometric Measurements . . . . . . . . . . 137xList of TablesTable 4.1 Detection significance of CHIME-eBOSS correlation . . . . . . . . 116xiList of FiguresFigure 1.1 Evolution of the energy densities over the age of the Universe . . 8Figure 1.2 BAO as measured by Planck . . . . . . . . . . . . . . . . . . . . . 12Figure 1.3 Summary of the eBOSS survey . . . . . . . . . . . . . . . . . . . 14Figure 1.4 Interferometry diagram . . . . . . . . . . . . . . . . . . . . . . . . 19Figure 1.5 BAO scale with CHIME . . . . . . . . . . . . . . . . . . . . . . . 23Figure 1.6 Aerial photograph of CHIME . . . . . . . . . . . . . . . . . . . . 25Figure 1.7 Diagram of CHIME signal chain . . . . . . . . . . . . . . . . . . . 26Figure 1.8 CHIME primary beam extent . . . . . . . . . . . . . . . . . . . . 27Figure 2.1 Schematic of the CHIME receiver pipeline . . . . . . . . . . . . . 34Figure 2.2 Histogram of truncation errors . . . . . . . . . . . . . . . . . . . . 36Figure 2.3 Ringmap live-viewer screenshot . . . . . . . . . . . . . . . . . . . 44Figure 3.1 Photo of the CHIME and Galt telescopes . . . . . . . . . . . . . . 47Figure 3.2 Amplitude and phase of the CHIME beam from a holographic ob-servation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48Figure 3.3 Holographic beam measurements from a few bright point sources 49Figure 3.4 Magnitude of decorrelation as a function of delay in the CHIMEcorrelator . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51Figure 3.5 Slopes in the amplitude of uncalibrated holography on meridian . 53Figure 3.6 Slopes in the amplitude of calibrated holography on meridian . . 54Figure 3.7 Slope fits as a function of frequency for Tau A . . . . . . . . . . . 55Figure 3.8 Slope fits as a function of frequency for Cyg A . . . . . . . . . . . 56Figure 3.9 Slope fits as a function of frequency for Cas A . . . . . . . . . . . 57Figure 3.10 First SVD mode for Tau A . . . . . . . . . . . . . . . . . . . . . . 58Figure 3.11 Second SVD mode for Tau A . . . . . . . . . . . . . . . . . . . . 59xiiFigure 3.12 Magnitude of the first and second modes of the Tau A SVD ex-panded into feed-frequency space. Within each cylinder, feed indexincreases from south to north. . . . . . . . . . . . . . . . . . . . . 60Figure 3.13 Residuals from the Tau A amplitude data to a two-mode recon-struction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61Figure 3.14 Single SVD mode for an independent decomposition of each cylin-der for Tau A . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62Figure 3.15 Decorrelation effect on meridian as a function of baseline . . . . . 63Figure 3.16 Source extent effect as a function of baseline . . . . . . . . . . . . 64Figure 3.17 Dynamic spectrum for a gated holographic observation of the pul-sar B0329+54 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 66Figure 3.18 Beam measurement derived from gated holography of the pulsarB0329+54 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 69Figure 4.1 Histogram of redshifts of Extended Baryon Oscillation Spectro-scopic Survey (eBOSS) Lyman-α forest samples and their overlapwith the CHIME data. . . . . . . . . . . . . . . . . . . . . . . . . 78Figure 4.2 Binned map of background quasars in the eBOSS survey . . . . . 79Figure 4.3 RFI environment at the CHIME site . . . . . . . . . . . . . . . . 80Figure 4.4 Delay power spectrum for CHIME visibilities . . . . . . . . . . . 87Figure 4.5 Diagram of the CHIME sidereal averaging pipeline . . . . . . . . 91Figure 4.6 CHIME beam model and data in m-space . . . . . . . . . . . . . 98Figure 4.7 Slices through the CHIME beam model . . . . . . . . . . . . . . . 100Figure 4.8 Delay power spectrum of a CHIME map . . . . . . . . . . . . . . 101Figure 4.9 CHIME maps at subsequent stages of filtering . . . . . . . . . . . 103Figure 4.10 Sample of corresponding eBOSS and CHIME spectra . . . . . . . 104Figure 4.11 Modelled parameter values described in this section as a functionof redshift. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 107Figure 4.12 Steps for generating simulated Lyman-α spectra . . . . . . . . . . 112Figure 4.13 Diagram of the 21 cm and Lyman-α simulation pipeline . . . . . . 113Figure 4.14 Cross-correlation templates from simulations . . . . . . . . . . . . 114Figure 4.15 Cross-correlation function of CHIME and eBOSS measurements . 115Figure 4.16 Cross-correlation functions evaluated using two methods . . . . . 117Figure 4.17 Correlation function evaluated on the difference of even and odddays . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 119xiiiFigure 4.18 Variance in the correlation functions of the full stack, even-odddifference, and a noise realisation . . . . . . . . . . . . . . . . . . 120Figure 4.19 Histograms of pixels in the filtered maps for the lower band . . . 121Figure 4.20 Histograms of pixels in the filtered maps for the upper band . . . 122Figure C.1 Diagram of spatial aliasing in interferometry . . . . . . . . . . . . 138xivList of AbbreviationsADC analogue-to-digital converterBAO baryon acoustic oscillationsCHIME Canadian Hydrogen Intensity Mapping ExperimentCMB cosmic microwave backgroundCPU central processing unitDFT discrete fourier transformDPSS discrete prolate spheroidal sequenceDRAO Dominion Radio Astrophysical ObservatoryeBOSS Extended Baryon Oscillation Spectroscopic SurveyEW east-westFFT fast fourier transformFGPA fluctuating Gunn-Peterson approximationFLRW Friedmann-Lemaître-Robertson-WalkerFPGA field programmable gate arrayFRB fast radio burstGPU graphical processing unitGR general relativityHA hour angleHI neutral hydrogenIGM inter-galactic mediumΛCDM Λ cold dark matterLSB least-significant bitLSS large-scale structureNCP north celestial poleNGC north galactic capNS north-southxvPFB polyphase filterbankQFT quantum field theoryRA right ascensionRF radio frequencyRFI radio frequency interferenceSDSS Sloan Digital Sky SurveySGC south galactic capSVD singular value decompositionUBC University of British ColumbiaVLA Very Large ArrayxviAcknowledgementsI am very grateful to my supervisor Gary for guiding me along the path of scientificinquiry, and through the intricacies of a PhD. Also to Mark, Richard, and Simon– it’s hard to appreciate just how much I have learned over these years by beingsurrounded and collaborating with a group of people possessing such a breadth ofexperience and skill. This sentiment extends to everyone in the UBC experimentalcosmology lab – my lab-mates over the years: Mateus, Carolin, Deborah, Sidhant,Kwinten, Arash, Sofia, Josh, Pedro, Parham, and Looney with whom I shared theups and downs of grad school; Mandana and Don who keep the ship afloat and alwaysprovided support; and there are still more people I was lucky to work with, but I’llleave it at that. The CHIME team at large has been a wonderful group to be a partof during these formative years.I thank my family and friends at home in Montreal for their support and patienceover these many years studying at the other end of the country. My parents made itpossible for me to have this opportunity, and I am happy to see the pride they mayrightfully feel in this achievement. My brother Sam is my best friend and the tripswe shared are highlights of this period in my life. Hearing from Maurice regularlywas a precious reminder of the friendships that never faded.Finally, I am also lucky to have friends and (almost) family in Vancouver, andI’d like to acknowledge all the people who have made living here so wonderful. Idiscovered many things about the world and myself over these years, and it is thanksto those who brought me into the mountains and guided me on this parallel path oflearning. I am grateful for the adventures and daily life I shared with Carolin andLooney, who formed a little family with me, and to Ellie and Ronnie who gave me ahome here. Love strikes without warning, and the tumultuous months that producedthis thesis were tempered by a sweet embrace. Merci Cha.xviiChapter 1IntroductionThis thesis reports on research that I performed at the UBC experimental cosmologylab working on the Canadian Hydrogen Intensity Mapping Experiment (CHIME).During my time working on CHIME, I participated at every step on the way to makingthis experiment work, including assembling the instrument, building the software thatperforms the data acquisition, calibrating and understanding the data, and finallymaking a detection of the cosmological signal it aims to measure. By no means amI solely responsible for any of these efforts, my contributions being in collaborationwith a team of student, staff, and faculty colleagues at UBC and the other partnerinstitutions. CHIME has not yet completed its scientific mission, and challengesremain to be overcome, but already it has been forging a novel path forward to learnmore about our Universe, and particularly dark energy, as I will explain below.The nature of working on a large experiment over a number of years is such thatprobably the bulk of the time I accumulated working on CHIME is not reflected inthis thesis, as it consisted of completing tasks that were necessary for the success ofthe experiment – and built up my knowledge of it – but are not worth reporting indetail here. Instead, I have chosen a few of the most notable contributions I madeto the experiment as chapter topics, culminating in a detection of cosmological 21 cmemission in cross-correlation with the Lyman-α forest.In this introduction, I will provide some background in cosmology and radio in-terferometry that will be helpful to understand CHIME’s scientific objectives andthe instrumental design chosen to achieve them. I will also give an overview of theCHIME project and its scientific context.11.1 Cosmology1.1.1 The Geometry of SpacetimeIn general relativity (GR) spacetime is not a static, Minkowskian backdrop on whichphysics proceeds, but rather is allowed to stretch and curve dynamically in responseto the presence of matter and other fields, which in turn evolve as directed by thisnon-trivial geometry. John Wheeler put it succinctly: “Spacetime tells matter howto move; matter tells spacetime how to curve” (Wheeler & Ford, 1998). In detail,the geometry of spacetime in the Universe is very complicated: we experience thegravitational effect of the Earth, which along with the other planets and Solar Systemobjects is orbiting the Sun, itself one of billions of stars making up the Milky Way,which in turn is part of a group of galaxies separated by vast reaches of emptiness,and so on... But eventually, at large enough scales and on average, we expect everypart of the Universe to look the same, from whichever direction we look at it. Thisis known as the Cosmological Principle – that on large scales, the Universe mustbe homogeneous and isotropic. Some background concerning the geometry of thehomogeneous and isotropic Universe will be reviewed here, but for a more detailedtreatment, see e.g. Dodelson & Schmidt (2020); Baumann (2022).The FLRW MetricThe geometry of spacetime is described by a metric, g, which locally determines theinvariant interval between points in spacetime,ds2 = gµνdxµdxν , (1.1)where xµ is a coordinate in 4-dimensional spacetime indexed by µ and the Einsteinconvention is used to indicate the sum over repeated indices. The requirement thatthe Universe be homogeneous and isotropic on large scales demands that the metrichave these properties, which constrains its possible forms significantly. The mostgeneral metric that satisfies this requirement is known as the Friedmann-Lemaître-Robertson-Walker (FLRW) metricds2 = −c2dt2 + a2(t)[dr21− kr2/R2 + r2(dθ2 + sin2 θdφ2)], (1.2)2where the spatial part is expressed in spherical coordinates (r, θ, φ), a(t) is a freefunction of time called the scale factor, which changes the relative scale of the spatialpart of the metric, and k determines the sign of global spatial curvature, with acharacteristic scale R,k =+1 Spherical,0 Euclidean,−1 Hyperbolic.(1.3)There is an overall redundancy in the definition of r and a(t) that allows us to referencethe scale factor to its value today: a(t0) = 1.Recent constraints on the global curvature from measurements of the cosmic mi-crowave background (CMB) combined with measurements of large-scale structure atlater times suggest that the Universe is spatially flat to a 1σ accuracy of 0.2% –in other words, the radius of curvature R must be at least ∼ 20 times that of theobservable Universe (Planck Collaboration et al., 2020). Given that the geometryrelevant to our Universe is observed to be Euclidean to a high degree of accuracy, Iwill set k = 0 in what follows, which will simplify the mathematics.The time-dependence of the scale factor means that distances at different stagesin the expansion (or contraction) of space are measured differently, so it is useful todistinguish physical distances from the coordinates: xphys = a(t)x. The coordinateswithout the scale factor are referred to as comoving coordinates, and they correspondto the physical distance that would be measured today (when a(t0) = 1).RedshiftFree-falling test particles in this dynamical space follow geodesics, defined byd2xµdτ 2+ Γµαβdxαdτdxβdτ= 0 , (1.4)where the Γµαβ are the Christoffel symbols for the FLRW metric and τ is an affineparameter. The parameter τ can be chosen to be proper time, such that dxµdτ= P µ =(E/c,p)µ is identified with four-momentum. With this definition,ddτ=dx0dτddx0=Ec2ddt. (1.5)3If you were to compute the Cristoffel symbols for the µ = 0 component of Equa-tion (1.4), you would find (Dodelson & Schmidt, 2020)Ec2dEdt= −aa˙δijP iP j = − a˙ap2 , (1.6)where Latin indices i, j span only the spatial components, and p2 = a2δijP iP j definesthe physical momentum (which includes a factor of a2). For a particle of mass m,E2−p2c2 = m2c4, including for photons with m = 0. This implies that EdE = c2pdp,sodpdt= − a˙ap = −Hp , (1.7)where I have defined the Hubble parameter, H ≡ a˙/a. The solution to this equationisp ∝ 1a, (1.8)implying that momenta decay as space expands. For photons, p = h/λ, and thiscorresponds to the wavelength of light stretching with the expansion. Given thatwe mainly observe the distant Universe by receiving light that has traveled largedistances, this effect is observationally very important. Light that was emitted witha rest wavelength λe at a time when the scale factor was a(te) can be measured todayat a shifted wavelengthλ = λe1a(te)≡ λe (1 + z) , (1.9)defining the redshift z. To the best of our current knowledge, the scale factor is amonotonically increasing function of time, so it can be used as a proxy for time. Thisis thus also true of redshift, which, as an observable quantity, is the most commonlyused measure of cosmological time/distance into the past.Hubble’s LawFor distances that are relatively nearby, such that the redshift is small, the scalefactor can be expanded in a Taylor series about its value todaya(te) = 1 + (te − t0)H0 + ... (1.10)where H0 ≡ H(t0) = a˙(t0) is the Hubble constant, the rate of expansion today. Theredshift is thenz =a(te)− 1a(te)≈ (te − t0)H0 = dcH0 , (1.11)4where d is the distance traveled by light in that time (neglecting the expansion overthis short period). This equation is a statement of Hubble’s law, the observationthat the measured recession velocity of distant galaxies (interpreting the redshift asa Doppler shift) is proportional to their distance – historically the first evidence forthe expanding Universe. Peculiar velocities – i.e. the galaxies not being at restwith respect to the frame in which the expansion proceeds – also contribute to themeasured redshift.DistancesBecause the metric is changing in time, measuring distances to faraway objects inthe Universe can be ambiguous and there exist a few different definitions that canbe useful in different contexts. In all cases, I am considering distance along the lineof sight of an observer located here today. The first is comoving distance, which issimply the coordinate distance that a photon would travel to reach us from a givenredshift,χ =∫ t0t1cdta(t)=∫ 1acda′H(a′)a′2=∫ z0cdz′H(z′). (1.12)A common way of measuring distance is from the apparent angular diameter ofan object of known physical size l:∆θ =ld. (1.13)In an expanding space, this relation is not quite correct, because the physical size wasdifferent at the time of emission. It can be restored by considering distances on thecomoving grid,∆θ =l/aχ, (1.14)which leads to the definition of angular diameter distance:∆θ =lDA; DA ≡ χ1 + z. (1.15)DA is thus an observable that can be inferred from the measured angular size of astandard ruler of known physical extent at a given redshift.Finally, luminosity distance is the distance that can be inferred by measuringthe apparent magnitude of standard candles, i.e. sources of known luminosity. In theabsence of an expanding background, the apparent flux from a source with luminosity5L at a distance d isF =L4pid2. (1.16)There are a few corrections to this equation from the expansion of space. First, thearea collecting the flux is assumed to be measured at the present day, so the comovingdistance should be used in the denominator. Second, in a unit of time, the numberof photons streaming through a spherical shell at the observer will be less than wherethey were emitted because they must travel a further distance, by a factor of a.Finally, the redshift also reduces the energy of the photons by another factor of a.Luminosity distance is defined such thatF =L4piD2L; DL ≡ χ(1 + z) . (1.17)Luminosity distance can be measured by observing the apparent flux from sourceswith a calibrated luminosity, a method that has been notably applied to Type Iasupernovae.Both angular diameter and luminosity distances are observables that have beenused to measure comoving distance as a function of redshift. These measurementsprobe the expansion of space over cosmic time and constrain the components of thestandard model of cosmology.1.1.2 The Accelerating Universe and ΛCDMThe dynamics of spacetime are determined by its stress-energy content, via the Ein-stein equationGµν =8piGc4Tµν , (1.18)where Gµν captures the curvature of spacetime and is composed of derivatives of themetric, and Tµν is the stress-energy tensor, accounting for the contents of spacetime.Components of the Stress-energyIt can be shown (e.g. Baumann (2022)) that the requirements of homogeneity andisotropy constrain the form of the stress-energy tensor to that of a perfect fluid, whichis diagonal and characterised by its energy density ρ and pressure P . The covariantexpression of local energy and momentum conservation leads to the continuity equa-tionρ˙+ 3a˙a(ρ+ P ) = 0 . (1.19)6Together with an equation of state that relates density and pressure,P = wρ , (1.20)the continuity equation describes how a fluid will evolve in the expanding FLRWspace. There are three types of fluids that are relevant over the course of cosmichistory.• The first is non-relativistic matter, for which P ρ, i.e. w ≈ 0. In this case,and in the absence of other components, Equation (1.19) givesρ ∝ a−3 , (1.21)which describes the dilution of the energy density as the volume increases.• Next is radiation, or relativistic matter with m E. For a gas of relativisticparticles, the equation of state is P = 13ρ, and we getρ ∝ a−4 , (1.22)with the extra power of a accounting for the effect of redshift.• Finally is the case of Dark energy, which has w = −1, i.e. negative pressure:P = −ρ. This leads toρ ∝ a0 , (1.23)a constant energy density. We have not directly observed any physical fieldsthat behave this way, but there is convincing indirect evidence that dark energyplays an important role in the expansion of our Universe. More on this later.The different powers of the scale factor in the evolution of each component ofenergy density result in stages of cosmic expansion that are dominated by one com-ponent at a time. This is illustrated in Figure 1.1.Friedmann EquationsThe Einstein equations (1.18) simplify greatly for the case of perfect fluids evolvingon the FLRW metric. The remaining independent equations are called the Friedmann710 5 10 4 10 3 10 2 10 1 100a10 510 2101104107101010131016/c,0matterradiationFigure 1.1: Evolution of the density of matter, radiation, and Λ dark energyover a range of scale factor that spans periods of radiation domination,matter domination and the current era of dark energy domination. Thedensities are given in units of the critical density today (Equation (1.26)).equations,H2 =(a˙a)2=8piG3c2ρ− ka2(1.24)a¨a= −4piG3c2(ρ+ 3P ) . (1.25)Although I had previously fixed k = 0, I reintroduce it here to highlight how it isrelated to the total energy density ρ. Equation (1.24) defines a critical density forwhich k = 0:ρcrit =3H2c28piG. (1.26)The critical density today, ρcrit,0 = 3H20c2/8piG is a useful reference that is conven-tionally used to define the relative densities of the different components of the energymakeup of the Universe:ΩX =ρX(t0)ρcrit,0, (1.27)for component X. The statement that our Universe is observed to be flat is that themeasured densities relative to the critical density sum to 11− Ωr − Ωm − ΩΛ ≡ Ωk < 0.2% , (1.28)8where the subscripts indicate non-relativistic matter, radiation, and dark energy. Ωk,the energy density of curvature, is defined to be the residual of this sum from unity.Another thing to note about curvature is that in Equation (1.24), its contributioncan be seen to scale as a−2. This means that its small value today implies it musthave been even smaller in the past.The Friedmann equations, along with the continuity equation, the relative densi-ties of each energy component today, and the Hubble constant, fully determine thedynamics of the homogeneous and isotropic Universe (at least as long as there isno significant change in the constituents, e.g. relativistic particles becoming non-relativistic).Dark EnergyThe second Friedmann equation, (1.25), describes the acceleration of the expansionof space. It is negative as long as ρ > −3P . For a single component, this is equivalentto w > −1/3, which is true of matter and radiation, but not dark energy. This is whythe expansion was expected to be decelerating, until measurements of the luminositydistance to distant supernovae revealed that actually the expansion is currently ac-celerating (Riess et al., 1998; Perlmutter et al., 1999). Since then, observations of theanisotropies of the CMB (e.g. Planck Collaboration et al. (2020)) and other probeshave confirmed the accelerated expansion and bolstered the case that dark energy isan inescapable part of the makeup of our Universe.Any model of dark energy that has w < −1/3 could plausibly explain the accel-eration, but the simplest and currently most favoured by observations is the cosmo-logical constant. The Einstein equations (1.18) have freedom to allow an additionalterm −Λgµν without breaking the conservation of stress-energy (since ∇µgµν = 0) forany constant Λ,Gµν − Λgµν = 8piGc4Tµν . (1.29)Λ is known as the cosmological constant, and it contributes to the overall dynamicsof the Universe (it was used by Einstein to engineer a static, but unstable Universe).Although I have included it in the Einstein equation here, it is more common nowa-days to integrate it into the stress-energy tensor as an additional component, one withw = −1. Such a constant energy density does have a counterpart in quantum fieldtheory (QFT) – the vacuum energy – but its expected value is ∼ 120 orders of mag-nitude larger than what is observed in cosmology (Baumann, 2022). This enormousdiscrepancy is called the cosmological constant problem.9If the observed dark energy is not consistent with the cosmological constant ex-pected from QFT, then something else might be responsible for the accelerated ex-pansion. Two possibilities are often considered as an alternative to the cosmologicalconstant: it could be that our assumption that GR is a good description of gravityon these large scales is incorrect, with the solution being a modified theory of gravity;or the dark energy could be composed of a new field that has yet to be discovered.Although these are theoretical possibilities, there is currently no widely accepted pro-posal that solves these problems consistently. For a review of the problem and sometheoretical considerations, see Particle Data Group et al. (2020) Chapter 28. In theabsence of a specific prediction for deviations from a cosmological constant, the darkenergy is often parametrised in terms of the equation of state parameter w, since itdetermines the rate with which it becomes dominant at late times. Deviations fromw = −1, or a w that evolves with redshift are theory-agnostic models that can beconstrained by observation.Progress on this topic can be expected through continued and more precise obser-vations of the expansion. Although current measurements are entirely consistent witha cosmological constant (Planck Collaboration et al., 2020), the late times when darkenergy becomes important (Figure 1.1) remain relatively sparsely observed. Manyexperiments are currently underway to probe this period of the expansion with un-precedented precision, and will either confirm that w = −1 or discover deviationsaway from the cosmological constant model. CHIME is one of these experiments,and its scientific design will be discussed in Section 1.3.ΛCDMThe standard model of cosmology is known as ΛCDM – Λ refers to the cosmologicalconstant, and CDM stands for Cold Dark Matter. Dark matter is a non-relativistic(“cold”) component of the stress-energy content of the Universe that behaves justlike regular matter gravitationally, but has no other apparent interactions that wouldallow it to be detected. The nature of dark matter is another open question incosmology.Despite being composed mainly of things we do not understand, ΛCDM is anextremely good description of the observable Universe, fitting observations spanningbillions of years of evolution consistently. The CMB has provided the most constrain-ing measurements of the contents of the Universe, which is composed today mainlyof dark energy (ΩΛ = 0.69), CDM (Ωc = 0.26), and baryonic matter (Ωb = 0.05),with a trace amount of radiation (Planck Collaboration et al., 2020). Together with10galaxy redshift surveys that probe structure at late times, these data are remarkablyconsistent with the standard ΛCDM model (Alam et al., 2021).In recent years, however, some tensions between observations using different trac-ers have begun to emerge. The most publicised are theH0 and σ8 tensions, which referto two parameters in ΛCDM for which recent measurements derived from late-timeobservations appear to be discrepant with those inferred from the CMB. A review ofthe current tensions in observational cosmology is provided in Abdalla et al. (2022).Although it remains unclear if these tensions point to unidentified systematic errorsor a real failure of the model, they highlight the need for more precise observationsfrom a variety of methods in order to continue to systematically test ΛCDM.1.1.3 BAO Standard RulerAs we rewind the clock on cosmic expansion, the temperature and density increase.At sufficiently early times, the temperature was high enough that matter was fullyionised, composed of free electrons and protons that were tightly coupled to thesurrounding radiation. In this environment, the mean free path between scatteringevents for photons was very short, so the Universe was essentially opaque until itcooled sufficiently for hydrogen to form, an event known as recombination. Followingthis moment, photons could free-stream essentially unimpeded through the expandingspace. These photons make up the CMB, and can be observed today to form animage of the surface of last scattering – the furthest back in time we can probe withelectromagnetic radiation.Before recombination, the competing effects of gravity and pressure create a restor-ing force that can support acoustic waves in the photon-baryon plasma, with a soundspeedc2s =13 (1 +R), (1.30)where R = (3ρb)/(4ργ) is the baryon-to-photon ratio. These are baryon acoustic oscil-lations (BAO). The initial density perturbations seed travelling waves that propagateat the sound speed until recombination, when the radiation pressure suddenly van-ishes. The CMB photons that are released record these acoustic oscillations, whichleft a statistical imprint in the pattern of the anisotropies across the sky at a scalecorresponding to the distance traveled until recombination, the sound horizon:rs =∫ t∗0csdta(t), (1.31)11Figure 1.2: (Top Left) Full-sky image of the CMB measured by Planck. (TopRight) Average of all of the hot and cold spots in the CMB map on theleft. The ring that statistically emerges from the random fluctuationsin the map is a manifestation of the BAO. The plots in the top roware reproduced from Planck Collaboration et al. (2014). (Bottom) Theangular power spectrum of the CMB measured by Planck (reproducedfrom Planck Collaboration et al. (2020)). In harmonic space, the discretebump at the BAO scale appears as “wiggles”, with a spacing determinedby rs.where t∗ is the time of recombination. The quantity rs has been precisely measuredin the CMB. Figure 1.2 shows the BAO feature both in the map of CMB anisotropiesand in their angular power spectrum. The physics of the BAO is exquisitely wellunderstood and the detailed analysis of CMB measurements has provided the mostpowerful constraints on the ΛCDM model. A review of BAO physics is outside thescope of this thesis, but can be found in, e.g. Hu (1995).12In addition to leaving a trace in the anisotropies of the CMB, BAO are frozen intothe distribution of baryonic matter following recombination. This feature manifestsitself as a statistical correlation in the distribution of matter at separations of rs. Thescale rs ≈ 150Mpc is large enough that it is stable to gravitational collapse over theage of the Universe, and just tracks the expansion of space. This makes it an excellentstandard ruler – a feature with a known size that can be measured at different stagesof the expansion to follow the evolution of the scale factor a(t).The observational method that has been used most successfully to measure BAOin the large-scale structure of the Universe at late times has been the surveying ofgalaxy redshifts. In this approach, a large sample of distant galaxies are observedspectroscopically, and the spectra are used to precisely determine the redshift ofeach galaxy. The result is a point cloud of measured galaxy positions in a large3-dimensional volume, with coordinates of angular position and redshift. Since thegalaxies form where matter is present, and more frequently where the density is higher,they can be used as a biased tracer of the overall matter density in that volume ofthe Universe. The BAO scale can be measured statistically in such a sample, asillustrated in Figure 1.3. The Extended Baryon Oscillation Spectroscopic Survey(eBOSS) survey1, which was recently completed, has measured BAOs at redshifts of0–2 and found good agreement with the predictions of ΛCDM (Alam et al., 2021).1.1.4 Hydrogen Intensity MappingGalaxy redshift surveys have been very successful probes of the large-scale structureat late times, and next generation surveys like DESI2, or slightly different approaches– e.g. DES3, Euclid4, HETDEX5 – will expand the volume being mapped this way.Scaling up these surveys, however, is expensive and time-consuming. Over about adecade, eBOSS measured spectra of around two million galaxies, using an instrumentthat could isolate 1000 objects in its field of view. Its successor, DESI is expectedto measure ∼ 40 million galaxies over 5 years, with 5000 robotically positioned fibresfor isolating objects. Higher redshifts are also more difficult to observe because ofincreasing atmospheric absorption at larger wavelengths.A complementary approach to mapping large-scale structure is Hydrogen Inten-sity Mapping. The technique uses modest-angular-resolution observations of red-1https://www.sdss4.org/surveys/eboss/2https://www.desi.lbl.gov/3https://www.darkenergysurvey.org/4https://www.esa.int/Science_Exploration/Space_Science/Euclid5https://hetdex.org/13Figure 1.3: Illustration of the completed eBOSS survey. The map in the centreshows a projection of the measured positions of galaxies in angle andlookback time (redshift is also shown in light blue), with different coloursindicating different samples. The plots on the right show the correlationfunction evaluated for these samples spanning a range of redshifts, with abump at the BAO scale clearly visible. Figure reproduced from eBOSS.shifted 21 cm emission from the hyperfine transition of neutral hydrogen to trace thedistribution of hydrogen gas, and thus matter, in the Universe. Hydrogen intensitymapping allows the apparent angular and radial BAO scale to be measured throughcosmic history without the expensive and time-consuming step of resolving individualgalaxies.The Redshifted 21 cm LineThe spin coupling of the electron and proton in neutral hydrogen (HI) leads to ahyperfine splitting of the electron’s energy levels. The transition between the hyper-fine levels of the ground state of HI has an energy corresponding to a photon withfrequency 1420.406MHz, or wavelength of about 21 cm. The probability of this tran-sition occurring spontaneously in a given hydrogen atom is extremely low, but this14is balanced by the large cosmic abundance of HI in such a way that extragalactic21 cm emission (and/or absorption) is measurable in aggregate. The lack of compa-rably strong spectral lines at frequencies below 1420MHz and the optical thinnessof the hyperfine transition together imply that we can directly observe a redshift ofthe 21 cm line at radio wavelengths. This can then be related to a distance from theobserver. Thus, maps of the radio sky at different frequencies contain informationabout the distribution of HI at different cosmic times, and the spectral and angularfluctuations of these maps can provide us with a three-dimensional picture of thisdistribution (Peterson et al., 2009; Scott & Rees, 1990).Intensity mapping consists of measuring the integrated emission at a relativelycoarse angular resolution to produce a continuous map of the sky. For a radio telescopemeasuring redshifted 21 cm radiation, its typical angular resolution covers an area thatwould contain many galaxies, as well as gas in the intergalactic medium, and it detectsthe sum of the emission from all of these objects together. For most applicationsin cosmology, and in particular for measuring the BAO scale, resolving individualgalaxies is not necessary because the relevant scales are much larger (∼ 100 kpcversus 150 Mpc). Another benefit of this method is that because 21 cm emission isobserved in aggregate, even faint sources that would be not be detectable individuallycontribute.At late times, radiation from stars has largely ionised the inter-galactic mediumand the remaining HI mainly resides in dense clouds within galaxy halos (Villaescusa-Navarro et al., 2018). This period of reionisation is thought to have occured at red-shifts around 20 & z & 6 (precisely how this proceeded has not been experimentallyconfirmed), and before that the Universe was entirely neutral and devoid of luminousstars. Intensity mapping with the 21 cm line is one of few observational techniquesthat can access these periods in the history of the expanding Universe, and has the po-tential to map enormous volumes at high redshift, unlocking a wealth of informationabout the Cosmos.ForegroundsAt the radio wavelengths of redshifted 21 cm emission, the sky is dominated by fore-ground emission from our Galaxy and extragalactic point sources. The brightnessof 21 cm emission is expected to be ∼ 5 orders of magnitude fainter than the fore-grounds at z ∼ 1. The foregrounds, however, are characterised by smooth frequencydependence, typically originating from synchrotron emission that produces a power-law spectrum. On the other hand, the 21 cm signal amplitude tracks the density15fluctuations of matter, and its frequency spectrum corresponds to the redshift axisin a tomographic map of these fluctuations, spanning a vast distance along the line-of-sight. The spectrum of 21 cm emission is thus characterised by structure on smallfrequency scales. These very different spectral properties can be used to separate thesignal from foregrounds (e.g. Liu & Tegmark (2011)).Intrinsically, the signal and foregrounds are separable, but observing them with aphysical apparatus will necessarily introduce an instrumental response as well. Anyspectral structure that is imprinted on the foreground power by the instrument cansignificantly impede signal recovery. The large dynamic range means that even asmall amount of leakage out of the smooth spectral modes where foregrounds areconfined can be enough to overwhelm the signal. Precise instrumental calibrationand analysis methods that attempt to mitigate these effects will therefore be key tothe success of 21 cm cosmology.Existing Measurements of Cosmological 21 cm EmissionA first detection of cosmological large-scale structure using 21 cm intensity mappingwas reported by Pen et al. (2009) in a cross-correlation analysis of existing radiodata from the HIPASS survey and the galaxy redshift catalogue 6dFGS. Subsequentanalyses have used dedicated observations by the Green Bank, Parkes and MeerKATradio telescopes, in concert with galaxy catalogues from the DEEP2, WiggleZ, 2dF,and eBOSS surveys. Chang et al. (2010) reported a detection at z = 0.8 in cross-correlation with DEEP2 at 4σ using observations made with the Green Bank tele-scope, followed by Masui et al. (2013) with a more precise measurement in cross-correlation with WiggleZ, also at z = 0.8. More recently, Anderson et al. (2018) usedthe Parkes telescope to detect the cross-power spectrum of 21 cm emission with the2dF survey at 5.7σ and z < 0.1, and Tramonte & Ma (2020) performed a separateanalysis of these data in map-space. Li et al. (2021) made a detection with signal-to-noise of 3σ at z = 0.75 using observations by the Parkes telescope and combined withthe WiggleZ survey. Wolz et al. (2022) cross-correlated the Green Bank 21 cm obser-vations with three galaxy samples from the WiggleZ and eBOSS catalogues resultingin detections from these three populations at z = 0.8. Finally, a first detection of21 cm signal using data from the MeerKAT telescope was reported by Cunningtonet al. (2023) in cross-correlation with WiggleZ, with a significance of 7.7σ at z = 0.4.In principle, much more powerful measurements of large-scale structure are pos-sible with custom-built telescopes that are optimized for 21 cm observations. Large,compact interferometers (see next Section) are particularly well-suited to observing16the 21 cm sky because they can be tuned to measure the cosmological scales of in-terest and can straightforwardly be scaled up in order to achieve a high sensitivityand fast mapping speed at low cost. Currently operating interferometers that werepurpose-built for this application notably include CHIME6 (described in detail inSection 1.3), targeting BAOs at redshifts 0.8 < z < 2.5, and HERA7, which aims tomeasure the epoch of reionisation at higher redshifts. CHIME recently reported thefirst 21 cm intensity mapping measurements made with an interferometer, in cross-correlation with the eBOSS galaxy redshift catalogues (CHIME Collaboration et al.,2022b). HERA continues to push down upper limits on the power spectrum of 21 cmfluctuations in the epoch of reionisation, with the latest results placing constraintson certain models of reionisation (HERA Collaboration et al., 2023).In contrast to these targeted experiments, the SKA8 project will be a general-purpose observatory consisting of two extremely large arrays of radio receivers. Us-ing data from an SKA pathfinder instrument, MeerKAT9, Paul et al. (2023) haverecently reported a detection of the 21 cm power spectrum at redshifts around 0.4,with > 8σ significance. This is the first Hydrogen intensity mapping detection madewithout resorting to cross-correlation with an external dataset, an important step inthe development of this method for precision cosmology.1.2 Radio InterferometersA brief introduction to radio interferometry is given here to provide necessary back-ground to the discussion of the CHIME project that will follow. In the first part Iwill describe conceptually how radio interferometry works, and then give some moreformal definitions in the second part. An excellent and in-depth overview of interfer-ometry in the context of 21 cm cosmology can be found in Liu & Shaw (2020).1.2.1 Beyond the ApertureTraditional telescopes are typically characterised by the diameter of the reflector theyemploy to collect and focus light to a receiver. The angular resolution achievable fora given diameter d at wavelength λ is given by the familiar relationship∆θ =λd. (1.32)6https://chime-experiment.ca/en7https://reionization.org/8https://www.skao.int/en9https://www.sarao.ac.za/gallery/meerkat/17One way to derive this expression is to consider how light arriving from a distantsource is collected and diffracted by the aperture of a parabolic dish reflector withdiameter d, measured at its focus. This geometry is illustrated in Figure 1.4. For aplane wave with wavevector k incident on the plane of the aperture, the electric fieldin that plane isE(r) = E0eik·r , (1.33)assuming a magnitude E0 for the plane wave (and not considering polarisation).Invoking Huygen’s principle, every point on the aperture can be considered as asource, from which radiation will be reflected off the parabola. At the focus, thepath lengths from every point on the aperture are equal, and thus contribute only acommon overall phase, so the total electric field is proportional to an integral overthe aperture, with the direction dependent part given byEfocus ∝∫Ad2reik·r . (1.34)This can be recognised as the Fourier transform of the aperture illumination (whichwas assumed to be uniform in this case), resulting in an Airy pattern for the diffractiondue to a circular aperture. Generally, the Fourier transform of an aperture withdiameter d cannot be more compact than ∆k ∼ 1/d or ∆θ ∼ λ/d. This pattern isoften referred to as the beam shape or point spread function of the telescope.At radio wavelengths, λ is on the order of meters, so building reflectors withdiameters sufficient to achieve high angular resolution quickly becomes impractical.An alternative is to build separate telescopes of reasonable size and place them farapart to observe the same source simultaneously, then combine their signals in sucha way as to start to approximate the response of a much larger aperture instrumentin a technique called interferometry. In a sense, such a pair of telescopes is like adiscrete sample in the plane of a very large aperture.In the first case where we considered a small continuous aperture, the reflectorcollects the light from every point in the aperture and combines it at the focus, withthe geometry such that each path from the aperture to the focus has the same length.What is measured in that case at the focus is power, or the magnitude squared ofthe electric field. Note, however, that in Equation (1.33) the fields in the aperturecarry information about the direction of the incoming plane wave in their phase. Thisinformation is lost when they are summed at the focus, and the way you select whichdirection to observe with such a telescope is to physically point it. If we sample theaperture instead, and record the phase information, then we can coherently detect18Figure 1.4: A plane wave signal from a distant source incident with wavevectork on two telescopes: a) A parabolic reflector of diameter d with a singleantenna at its focus. Two rays are drawn to show the geometry, but allthe light collected by the aperture contributes to the signal measuredat the focus. b) A pair of antennas lying in the aperture plane, withtheir signals fed to a correlator. In this case, only the two rays drawncontribute to the measured signal. The blue lines represent surfaces ofconstant phase, and in red, a transmission line from each antenna to thecorrelator.radiation incident from any direction.The separations between receivers are called baselines, and if we think of them assamples in the aperture space, they are filling in the integral on the right hand sideof Equation (1.34). Taking the Fourier transform of the baselines thus emulates thecollection of light from the aperture at the focus and produces an angular responsepattern. A key difference is that we are free to choose any weights for the sum, so wecould for example arrange for signals from off-axis to add coherently – pointing thetelescope after the fact. This is a sketch of the technique of aperture synthesis, andthe resulting beam pattern is called the synthetic beam.To summarise the conceptual discussion in this section: interferometric measure-ments of radio emission from the sky can be thought of as discretely sampling anaperture that spans the largest separation between receivers, thus overcoming themechanical limitations of building extremely large aperture telescopes for radio wave-19lengths. To make it possible to coherently combine the signals from these discretesamples, their phase must be recorded, and the specific choice of weights used in thiscombination allows one to effectively point the array in any direction at the dataanalysis stage. The more baselines are available to fill in the synthetic aperture, themore precisely the synthetic beam response can be shaped.1.2.2 VisibilitiesIn this section, I will define mathematical expressions for the interferometric mea-surements that were just discussed conceptually. I will use the formalism developedin Liu & Shaw (2020).First, consider the signal measured by a single antenna exposed to the sky:si(t) =∫d2nˆdνAi(nˆ, ν) · (nˆ, ν)e2piiν(t−xi·nˆ/c) , (1.35)where (nˆ, ν) is the electric field density from direction nˆ at frequency ν, Ai(nˆ, ν)is the primary beam response of the antenna, and xi is the position of the antenna.The signal si measured for antenna i is the projection of the polarisation direction ofthe incident electric field density onto the beam response of the instrument to eachpolarisation axis, represented by the dot-product of those two vector quantities. Theintegral runs over angular position nˆ and frequency ν. The phase term depends ontime and the location where the fields are measured.The emission from the sky (nˆ, ν) has some important properties. It is producedby astronomical processes involving an enormous number of independent emitters(electrons, atoms, etc.), and what reaches us is the sum of their radiation. By thecentral limit theorem, their emission is thus well described as a Gaussian randomprocess. Additionally, these processes are incoherent, i.e. they are not correlatedin space because they originate from the motion of particles with random initialconditions (note that this is not true of the polarisation components, which can becorrelated due to, e.g., alignment of magnetic fields generating synchrotron radiation),and their correlations in time are stationary, i.e. statistically independent of time.These properties imply that the electric field density can be described statistically bya two-point function that is diagonal in space and frequency (Liu & Shaw, 2020)〈a(nˆ, ν)∗b(nˆ′, ν ′)〉 ∝ Tab(nˆ, ν)δ(nˆ− nˆ′)δ(ν − ν ′) , (1.36)where Tab is a coherency matrix that describes the correlations between polarisations20a and b. For a brief derivation of why stationarity implies spectral incoherence, seeAppendix B. Although polarisation is an important component of this formalism andI include it for completeness, I will not delve into it here because it is not essentialto this discussion and beyond the scope of the thesis. It will be sufficient to keep inmind that these measurements include polarisation, but we are mostly interested inthe total intensity, trT.Given the emission is incoherent spectrally, it is useful to consider the Fouriertransform of the signalsi(ν) =∫d2nˆAai (nˆ, ν)a(nˆ, ν)e−2piixi·nˆν/c . (1.37)Finally, recalling that the spatial phase is the key element of interferometry, it can berecorded by cross-correlating the signals from antennas at different locations xi andxj:Sij(ν) =〈si(ν)s∗j(ν)〉=∫d2nˆAai (nˆ, ν)Ab∗j (nˆ, ν)Tab(nˆ, ν)e−2piiuij ·nˆ , (1.38)where the angle brackets indicate an average (typically in time), and we have definedthe baseline between the two antennas in units of wavelength as uij = (xi − xj)/λ.These two operations – the frequency transform and the spatial correlation –combine signals from an array of antennas interferometrically. Performed in thisorder, they constitute an “FX”-type correlator, the instrumental implementation ofwhich is described in some detail in Liu & Shaw (2020) and in the specific case ofCHIME in Section 1.3.1. In a real instrument, the signals also include a noise term,which is added to the cross-correlations. The measured cross-correlations are calledvisibilities:Vij = Sij +Nij , (1.39)where Nij is the noise correlation matrix. Note that for separate signal chains thenoise is in principle uncorrelated, which means that the expectation value of Nij isdiagonal, i.e. the visibilities do not have a noise bias.1.3 CHIMEThe Canadian Hydrogen Intensity Mapping Experiment (CHIME) is a large radiointerferometer that was designed and built by a collaboration of researchers fromCanadian universities and research facilities (the core institutions being the University21of British Columbia, McGill University, the University of Toronto and the DominionRadio Astrophysical Observatory (DRAO)). The main motivation for CHIME is toobserve a period in the Universe’s history that can shed light on the nature of darkenergy, the existence of which has been theorised in order to explain the observedaccelerated expansion of space (Section 1.1.2). The frequency band of 400–800MHzcorresponds to 21 cm emission from redshifts 2.5 ≥ z ≥ 0.8, where the simplest modelfor dark energy predicts it becomes a non-negligible constituent of the total energydensity, driving the transition to acceleration. CHIME will measure the expansionrate over this period of time by using the BAO scale as a standard ruler. As shown inFigure 1.5, CHIME’s angular resolution of ∼ 40′ and frequency resolution of 390 kHzare well suited to measuring the BAO scale in 21 cm emission over this range ofredshifts. An overview of the CHIME project was published in CHIME Collaborationet al. (2022a), and some material from that work, including text and figures, was usedto construct the introduction to the experiment that follows.CHIME is located at the DRAO in Penticton, BC. The site is in the White LakeBasin, within the traditional and unceded territory of the Syilx/Okanagan people. Itis composed of four cylindrical reflectors oriented north-south, each instrumented with256 dual-polarisation antennas arranged uniformly along the focal lines separatedby about 0.3m. There are no moving parts and the sky is observed as it driftsthrough the telescope’s broad field-of-view that spans the south and north horizons,focused to within a few degrees of the meridian in the east-west direction. The signalsare fed to an FX-style digital correlator, where they are first Fourier-transformedinto 1024 frequency channels (across the 400–800MHz band) before being multipliedpair-wise to form a correlation matrix. The spatial correlations between antennascan be combined to reconstruct angular structures on the sky, effectively pointingthe telescope after the fact in any direction within the broad field of view of thecylindrical reflectors. This allows CHIME to rapidly map the sky at high sensitivityby observing the entire survey field every day.The CHIME telescope achieved first light in the fall of 2017, and has been ac-quiring data ever since. In the context of the nascent field of 21 cm cosmology thatwas introduced in Section 1.1.4, CHIME is a pathfinder experiment. The effortsdevoted to overcoming the challenges associated with the bright foregrounds will pro-vide valuable insight in order to help realise the enormous potential of this method foradvancing cosmology. Analysis efforts to date have been focused on instrumental cal-ibration and pipeline development, which has led to a recent detection of HI emissionin cross-correlation with Sloan Digital Sky Survey (SDSS) galaxy redshift catalogues22012345Angular Scale (degrees)100m Baseline20m BaselineBAO Scale1/2 BAO Scale1/3 BAO Scale0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0Redshift01020304050Frequency Scale (MHz)1420 946 710 568 473 405 355 315 284Frequency (MHz)Figure 1.5: The BAO scale (rs ≈ 150Mpc, comoving) compared to CHIME’s angular andfrequency resolution. (Top): The blue solid curve shows the angular scale associatedwith rs, while the other linestyles show the first few harmonics (corresponding tothe peaks of successive BAO “wiggles” in k⊥ Fourier space, located at multiples ofkBAO ≈ 2pi/rBAO). The shaded region shows the range of angular scales accessibleto CHIME as a function of frequency, for antenna baselines ranging from 0.3–100m.The grey straight lines show the angular resolution associated with feed separationsof 20m and 100m. (Bottom): the solid curve shows the frequency separation as-sociated with the line of sight BAO diameter for 21 cm radiation as a function ofredshift. The other linestyles indicate the frequency resolution required to resolvethe first two BAO harmonics in k‖ Fourier space. CHIME’s frequency resolutionis 390 kHz, which is oversampling the BAO scale, but useful to excise narrowbandinterference from terrestrial sources. For all curves, we take H0 = 70 km s−1Mpc−1,Ωm = 0.3, and ΩΛ = 0.7. In both panels the shaded region denotes the frequency andredshift coverage of CHIME. This figure and caption were reproduced from CHIMECollaboration et al. (2022a).23at ∼ 10σ significance (CHIME Collaboration et al., 2022b), as well as the detectionin cross-correlation with the Lyman-α forest that is reported in this thesis. In thesections that follow, the instrument is described in more detail, as well as operationsthat are performed on the data in real time during the acquisition stages, before it issaved to disk and archived.1.3.1 InstrumentThe CHIME telescope is pictured in Figure 1.6. The parabolic reflecting surface ofeach of the four cylinders is 20m aperture, 5m focal length and 100m long. 256dual-polarization antennas are placed along the central 78m of the focal line of eachcylinder, beneath an 88m× 0.65m groundplane. The focal lines are covered by andsuspended from the walkways visible along each cylinder axis. Signals are amplified ateach feed and brought by low-loss coaxial cables to receiver huts located in commer-cial radio frequency (RF) shielded rooms within customized RF-protective shippingcontainers, one located between the first and second cylinders, another between thethird and fourth. After band-defining amplification, analogue-to-digital conversion, atime-to-frequency transform and half of a corner-turn, signals are brought from thetwo receiver huts to additional RF rooms within the white shipping containers seenat the right in the photo, where the corner-turn is completed and a spatial correlationand other processing are performed. The data flow through the instrument is illus-trated in Figure 1.7, which is divided into four main components, three of which willbe briefly described here: the analogue chain; the F-engine; and the X-engine. The re-ceiver pipeline will be discussed in the following section. A more detailed descriptionof these systems can be found in CHIME Collaboration et al. (2022a).Optical Design and Primary BeamThe design of CHIME is focused on enabling the measurement of BAO across theredshift range where dark energy begins to impact the dynamics of the Universe.At the wavelengths of the corresponding redshifted 21 cm emission, sufficient angularresolution to resolve BAO features on the sky is easily achieved by a 100m baseline(see Figure 1.5).An east-west (EW) array of cylindrical, 100m-long reflectors each coupled to alinear feed array along its focus meets these needs. Such a system scans a north-south (NS) stripe of the sky interferometrically and observes most of the 3/4 ofthe celestial sphere visible from our site every day as the Earth turns. Given that24Figure 1.6: A photograph of CHIME looking north. The two people standingat the southeast corner (right, foreground) help to show the physicalsize. CHIME is located at the DRAO which is protected by law and theadjacent hills from terrestrial radio interference. In the background onecan see five dishes of the DRAO Synthesis Telescope and a solar radiomonitor. Figure and caption reproduced from CHIME Collaborationet al. (2022a).each feed in this system produces a response of ± 1 rad along the cylinder axis,choosing a reflector shape to be an f/D = 0.25 parabola allows the use of feeds withapproximately symmetric angular response patterns. At this f-ratio, the focus is levelwith the edges of the reflector, protecting the feed array from terrestrial radiation.The primary beam is the response to the sky (amplitude, phase, and polarisation)of a single feed, including its coupling to all other feeds along the focal line. Thefeeds are oriented such that the orthogonal polarisation axes are aligned east-west andnorth-south, referred to as the X and Y polarisations, respectively. In the absenceof multi-path effects, discussed below, the primary beam produces a nearly ellipticalillumination of the sky: ∼120 degrees long in the unfocused NS direction, alongthe cylinder axis, and a few degrees wide with frequency-dependent diffraction side-lobes in the EW direction, perpendicular to CHIME’s cylinder axis. This shape isillustrated in Figure 1.8.Multi-path and other coupling effects alter this simple description by as much25F-EngineADC + FT+ corner turnFocal line(Linear dual-polarized feed array+ LNAs )ReflectorX-EngineMultiply + AccumulateReceiverFlag + Calibrate + CompressCorrelatorFigure 1.7: Schematic diagram of data flow through CHIME. Signals focusedonto a linear array of broad-band dual-polarization antennas, are am-plified (each polarization separately) using room temperature receiverswith a noise performance below 30K that amplify and filter the signalsto 400–800MHz. The correlator is an FX design, where the F-engine dig-itizes and channelizes the signals from the 2048 analogue receivers andalso implements the majority of a corner-turn network that rearrangesthe channelized data for spatial correlation. The X-engine completesthe corner turn and performs the cross-multiplications and averaging tocompute the N2 spatial correlation matrix separately at each frequency.Finally, a receiver node collects the streams from every frequency andperforms additional processing in real time, including calibration andcompression steps. Figure and caption reproduced from CHIME Collab-oration et al. (2022a).26Figure 1.8: Model for the primary beam response of a single CHIME feed,evaluated at 476.6MHz. Multi-path effects due to reflections within thecylinder and coupling between feeds are not included in this model. TheMollweide projection spans the full 4pi steradian, centered on zenith atCHIME, with North at the top. The beam response is zero below thehorizon, and concentrated in a North-South stripe along the meridian.The East-West profile is characterised by the diffraction of the 20maperture in that direction.as 50% at some frequencies. The physical origin of the multi-path interference isradiation interacting with the focal-line assembly, which consists of the linear feedarray and a common ground plane. This can be understood most easily by consideringthe telescope in transmitting mode, which is just the time-reverse of the receivingmode. In this environment, a signal broadcast by a feed will reflect off the cylinderand a large fraction of that signal will go directly to the sky, but a small portionstrikes the focal plane assembly, where some is absorbed by a neighbouring feedand the rest is reflected and/or re-radiated by the assembly, eventually reaching thesky. The details of this latter interaction are complex and are still actively beingcharacterized. Nonetheless, the primary beam is the illumination on the sky onegets when these effects are accounted for. Since multi-path propagation is occurringwithin a 5m cavity (CHIME’s focal length), new interference fringes arise roughlyevery 30MHz in frequency.27Analogue ChainThere are 256 antennas on each of the four cylinders, sensitive to two orthogonalpolarisations, resulting in 2048 analogue signal chains in total. The antennas employa cloverleaf design that was optimized specifically for the 400–800MHz octave ofbandwidth (Deng, 2020; Deng & Campbell-Wilson, 2014). The signals are amplifieddirectly at the antenna feed on the focal line by a first stage low-noise amplifier whichadds about 20K of noise power across the band. 50m coaxial cables carry the signalsdown from the focal lines to the receiver huts (see photo, Figure 1.6). The 2048cables were cut to the same length to 0.1% accuracy to ensure the signals arrive atthe receiver huts with a uniform delay. A band defining filter and a second stageof amplification are applied at the receiver hut, prior to digitisation. The bandpassfilter is designed to have sharp edges with &30 dB suppression out-of-band to enablesampling in the second Nyquist zone, which would otherwise be confused with signalsfrom the 0–400MHz band.F-EngineThe F-engine is responsible for digitising the signals from the analogue chain andperforming a frequency transform to the 400–800MHz band. It is implemented onthe ICE platform (Bandura et al., 2016). Each ICEboard hosts a field programmablegate array (FPGA)-based processor and two analogue-to-digital converter (ADC)daughter boards, which receive a combined 16 analogue signals. The ADCs sampleat 800 million samples per second, and this together with the bandpass filter allowsthe 400–800MHz band to be measured in the second Nyquist zone. 2048 samples areaccumulated and transformed to 1024 frequency channels by way of a polyphase fil-terbank (PFB). The implementation of the PFB is detailed in CHIME Collaborationet al. (2022a), and computed by the FPGA. A PFB is used instead of a standardFourier transform because it results in frequency channels with more compact sup-port.A total of 128 ICEboards are needed to digitize the 2048 signals, and these aredivided between eight crates, which provide power, networking and cooling to theboards. There are four crates in each receiver hut, serving two cylinders. The nextstage of the correlator – the X-engine – performs the spatial correlation and is dis-tributed over 1024 graphical processing unit (GPU)s. Each GPU processes a singlefrequency, and needs access to all inputs, whereas each FPGA has all frequencies foronly 16 inputs. A significant shuffling of the data needs to happen during the transferto the X-engine, which has been implemented in four hierarchical stages of network28communication called the corner-turn: first involving intra-board communication,then intra-crate, followed by inter-crate, and a final shuffle between the GPU nodesof the X-engine. This reordering of the data operates at 6.6Tbit/s.X-EngineThe X-engine is composed of 256 computing nodes, each of which hosts four GPUsand one central processing unit (CPU). These are connected to the F-engine viaa fibre optic network, with one fibre per GPU. The bundles of fibres are strungthrough ducts that link the receiver huts to another pair of RF-shielded huts wherethe GPU nodes are housed. Each GPU computes the products of every pair of inputsto form the visibility matrix. This operation constitutes the bulk of the computing inthe CHIME correlator, and a highly optimised kernel was developed to perform thistask on this specific model of GPU. Visibility samples at the native 2.56 µs cadenceare accumulated by the X-engine into 31ms frames to be passed on to the receiverpipeline for processing. Before the accumulation in time is complete however, radiofrequency interference (RFI) contaminated samples are identified by computing thespectral kurtosis over 256 time samples (0.66ms) and all inputs to flag non-Gaussiansignals (Taylor et al., 2018). Above a given threshold, these sub-frames are excisedfrom the data and the weights are adjusted accordingly so that the excision can betaken into account in subsequent processing.1.3.2 Real-time Calibration and FlaggingThe data produced by the correlator is a time stream of raw visibilities at 31ms ca-dence for every frequency, referred to as the N2 visibilities, where N = 2048 is thenumber of inputs. More precisely, only the upper triangle of the visibility matrix iscomputed, which contains N(N + 1)/2 products. This is a large data rate, amount-ing to over 500GB/s, and even after averaging frames to a final cadence of 10 s (overwhich the amount of Earth rotation is small enough compared to the resolution of thetelescope), saving N2 visibilities to disk is not feasible. The data are compressed bytaking advantage of the significant redundancy in the baseline array: the rectangulararray of antennas contains many separations that have the same length and orien-tation (especially the short ones) and therefore nominally observe the same signalfrom the sky. These redundant baselines are averaged together to reduce the numberof products from ∼ 2 × 106 to around 16 × 103. Such a compression is impossibleto undo, so any non-redundancies or contaminated signals must be flagged out or29corrected beforehand.The calibration and flagging stages of the real-time data processing pipeline thatenable baseline compression are summarised in this section, as well as a stage thatestimates noise using the fast-cadence samples. A more detailed account of these stepscan also be found in CHIME Collaboration et al. (2022a). These stages of processingare distributed across the GPU node CPUs and a dedicated receiver node. Chapter 2will provide a more technical description of this pipeline and its implementation, todiscuss specific components for which I led the development, whereas here I will onlyfunctionally describe stages of the pipeline that are important to understand theCHIME data.Gain CalibrationThe analogue chain imparts an amplitude and phase to the signal, which in generalis different for every input. Referred to as the complex gain, it must be calibrated forevery one of the 2048 inputs before visibilities can be averaged together.The complex gain of each input is calibrated once per sidereal day by fitting amodel to the eigen-decomposition of the N2 visibility matrix during the transit of thebrightest radio source that is available at night. This scheme relies on the calibrationsource being the dominant component of the sky visible to CHIME during its transit,and the Sun is bright enough even off axis that this approximation is not valid duringdaytime. The primary calibration source is Cygnus A because it is the brightestradio point source in the sky between 400 and 800MHz (Perley & Butler, 2017). It isalso unresolved by the longest CHIME baselines and has a stable, well-characterizedspectral flux density. Cassiopeia A, Taurus A, and Virgo A are used as alternativecalibration sources when Cygnus A is transiting during the day. The complex gainsare scaled by the flux density of the calibrator source, such that application of thegains converts the visibility data to units of Jansky/beam. The flux density of thesesources was measured with the Karl G. Jansky Very Large Array (VLA) in 2014 and2016 at frequencies ranging from 220MHz to 48.1GHz. VLA legacy observations from1998 also exist at 73.8MHz for all sources but Casseopia A. These measurements areinterpolated to the CHIME band using the polynomial expressions provided in Perley& Butler (2017).Between calibration solutions, the gains will drift due to changes in the environ-mental conditions, mainly temperature. However, given that the entire instrument isexposed to the same environment, these changes are largely common-mode betweenall inputs. What matters for redundant baseline compression are relative errors in30the gains within the array. The common-mode drift in the gains is corrected using athermal model as part of later offline processing (CHIME Collaboration et al., 2022b).Bad Input FlaggingFor a variety of reasons – e.g. malfunctions, disconnection – some inputs need to beexcluded from the baseline collation. Some of these are permanently disabled, butothers can spontaneously exhibit out-of-the-ordinary response and eventually recover.The role of the flagging stage of the pipeline is to perform real-time identification ofcorrelator inputs that should be excluded from further analysis. This runs on a sepa-rate processing node and provides regular updates to the relevant stages of the receiverpipeline. It uses a variety of data products and housekeeping metrics to repeatedlyevaluate 10 different tests, with each test designed to identify malfunctioning or other-wise anomalous correlator inputs. These steps are described in CHIME Collaborationet al. (2022a).Fast-cadence Noise EstimateThe real-time pipeline estimates the variance of the visibility for each baseline, fre-quency channel, and ∼ 10 s integration by differencing the even and odd 31ms sub-integrations. Since the observed foregrounds do not change significantly on 31mstimescales, they cancel for this difference, leaving contributions from RFI and intrinsicradiometric noise. This “fast-cadence estimate” of the variance is propagated througheach stage of the real-time and offline pipeline. When I refer to “noise weights” inlater chapters that discuss CHIME data, I will be referring to this estimate.31Chapter 2Receiver Pipeline DevelopmentThe receiver pipeline is a component of the CHIME correlator that processes in realtime data streaming off the X-Engine (as explained in the introduction, and illustratedin Figure 1.7). It performs data reduction and compression steps in addition toflagging and calibration, before packaging the data into a final archive format andwriting it to disk. These stages are summarised in Figure 2.1 (reproduced fromCHIME Collaboration et al. (2022a)).The pipeline is built on the kotekan1 framework, a high-performance data pro-cessing pipeline developed specifically for CHIME but now made open-source andused by other radio experiments (Renard et al., 2021). It is written in the C++ lan-guage and has a modular design centred around buffers and processing stages that cansupport multiple parallel streams of data. Performance requirements are stringent inorder to keep up with CHIME’s large data rate. At its peak, the total streaming datarate (distributed across all GPU nodes) reaches 500GB/s.Multiple instances of kotekan are distributed across an array of computer nodesconnected by a high-speed network (Figure 2.1). The most numerous are the GPUnodes, of which there are 256, each housing four GPUs for a total of 1024. A GPUnode receives the data from all feeds at a single frequency and computes their cross-correlations, the N2 matrix of visibilities. These data streams are passed on to thenode’s CPU where they may be additionally processed and are integrated to a 5or 10 s cadence. The data from all the GPU nodes are sent over the network toa receiver node where frequencies are collected into a single stream. At this stage,various subsets of the data are formed, and flagging and gain calibration are applied.These are generated in real time based on information derived from the streamingdata and forwarded to a separate processing node. The data are finally written to1https://github.com/kotekan/kotekan32disk on the receiver node. The pipeline is described in much greater detail in CHIMECollaboration et al. (2022a).During the commissioning phase of CHIME, I was extensively involved in thedevelopment of various components of the pipeline, both in their design and imple-mentation. In this section I highlight a few parts for which I led the development.2.1 Final Data Compression StageThe data written to disk amounts to ∼2.7TB per day, before the final compressionis applied prior to archiving. A compression scheme was developed to take advantageof the slowly varying nature of the CHIME output in time. As the sky drifts throughthe telescope’s field of view, the signal is highly correlated between neighbouringtime samples, only changing in its lowest significant digits. Given this property ofthe data, significant compression improvement can be achieved by transposing thebitwise representation of a time stream so that the redundant most-significant bitsare adjacent. Then any compression algorithm can efficiently pack those data. Thisalgorithm has been implemented by Masui et al. (2015) in a code named bitshuffle.My task was to implement the processing stages to reorder the data, prepare it forcompression, apply the compression, and package it into a standard archive format.The CHIME data are multidimensional – with axes for time, frequency, and baseline.In order for the bitshuffle compression to be effective, the data must be orderedsuch that time is the fastest-varying axis, i.e. adjacent time samples are packed insequence. However, by necessity data are acquired with time as the slowest-varyingindex, so the first step is to transpose the data. For this reason, the compressionpipeline cannot operate in real time with the rest of the receiver pipeline, and insteadwaits until a certain amount of data (a fixed-size file) is available to be transposed.Although the most-significant bits change very little, at the other extreme theleast-significant bit (LSB)s are dominated by noise, which does not compress. Tomitigate this and further improve the compression ratio, we truncate the LSBs, givenan estimate of the noise level in the data, and a specified limit on the truncation noisethus introduced. A routine was developed for this purpose that truncates floatingpoint numbers in their binary representation.2.1.1 Floating Point TruncationIn accordance with the Institute of Electrical and Electronics Engineers standard(IEEE, 2019), 32-bit floating point numbers are represented in three parts: the sign33x 1024Merge singlefrequency streamsPrimaryData FlowCalibration/FlaggingSecondaryData FlowNetworktransferN2 @ 30ms Integration, noiseestimation, gatingGated data @ 5sNormal @ 5s Time integration to10sEigenmodecalculationDisk outputN2 all freqSubset 26mproductsApply feed flags Apply gainsBaseline collationGain estimateMerge singlefrequency streamsMerge singlefrequency streamsDisk outputSubset 26mproductsMerge singlefrequency streamsfreq fromother GPUsfreq fromother GPUsfreq fromother GPUsDisk outputDisk outputFlagging brokerExternal GainInformationADC samples, amplifierstatus and other externalinfoCalibration brokerReal time mapmakerGPU node: singlefrequencyProcessing nodeShared memoryFreq subsetSelect eigenmodesand timing dataReceiver node: allfrequenciesFigure 2.1: Schematic of the CHIME receiver pipeline (reproduced fromCHIME Collaboration et al. (2022a)). The inner rectangles representdata-processing tasks and the ovals intermediate data products. Theouter rectangles indicate on which computer node the tasks are operat-ing. The input of the pipeline is the N2 correlation matrix computed onthe GPUs at a 31 ms cadence. The pipeline branches into multiple datastreams that undergo various processing. In the end four data productsare output to disk, and a separate data stream is forwarded to the real-time map maker. The processing tasks that are described in this sectionare highlighted with a red dashed outline.34bit (1 bit), the exponent (8 bits), and the mantissa (23 bits). The mantissa encodesthe fraction in base-2 and is scaled by the exponent to obtain the numerical value:(−1)s 2e−127(1 +23∑i=1b23−i2−i), (2.1)where s ∈ {0, 1} is the sign bit, e ∈ {1, 2, ..., 254} is the exponent, and b23−i ∈ {0, 1}are the successive bits of the mantissa. The extreme values of the exponent, {0, 255},are reserved for special values: zero and infinity/“not a number”, respectively.In order to truncate the bitwise representation of a 32-bit float to a given decimalprecision, we need to find a new float with the maximum number of leading zerospossible without deviating from the original value by more then the required precision.Care must be taken to handle cases where the truncation will require a change to theexponent (e.g. 3.8→ 4), and to avoid introducing a statistical bias to the data.The same problem applied to integers is more straightforward. In that case, onecan identify the largest power of two that is smaller than the allowed error and maskthe bitwise representation to that level. Therefore, the first step of the procedure fortruncating floats that I developed is to cast the mantissa and the precision to integers.1. Interpret the sequence of bits of the mantissa as an integer, corresponding to∑23i=0 b23−i223−i, where b23 = 1 is implicit in the float definition.2. The previous step is equivalent to multiplying the mantissa by 223, and we haveremoved the exponent, i.e. divided by 2e−127. To cast the precision to an integeron the same scale, multiply it by 223+127−e.3. Truncate the integer representation of the mantissa to the scaled integer preci-sion.4. Cast the truncated value back to a float by performing the inverse of the op-erations described above. Cases such as truncation to zero must be accountedfor.If the truncation procedure is unbiased, the truncation errors should be strictlysmaller than the required precision (or more precisely, the nearest fractional power of2) and uniformly distributed within. This was verified for the float truncation algo-rithm by truncating large arrays of random numbers and analysing the distributionof the errors. An example is shown in Figure 2.2.350.004 0.003 0.002 0.001 0.000 0.001 0.002 0.003 0.004truncation error02505007501000samples per binprecisionFigure 2.2: Distribution of truncation errors for a random sample of N = 10532-bit floats uniformly distributed between 0 and 1, plotted as a 100-binhistogram. The truncation precision was set to 2−8, as indicated by thedashed vertical lines. The variance between bins is consistent with auniform distribution, with a probability of measuring a higher varianceof 0.7.Assuming the sample of truncation errors is drawn from a uniform distribution,the variance of the number of samples that fall into each histogram bin should begiven by the binomial distribution:σ2 = Np(1− p), (2.2)where N is the number of samples and p, the probability of landing in that bin, is justthe inverse of the number of bins. The sample variance measured from the histogramfollows a χ2 distribution, which we can use to evaluate the probability of observinga variance exceeding that value. For the data shown in Figure 2.2, with N = 105samples in 100 bins, I find a probability of measuring a higher sample variance of 0.7,which suggests the truncation errors do not deviate from a uniform distribution.Truncating the data to a given precision e introduces additional noise, which needsto be characterised. For a signalv = v0 + n+ nt , (2.3)with true value v0 corrupted by measurement noise n, we can represent the errorintroduced from truncation as an additional noise term nt. As discussed, nt is uni-formly distributed within ±e with zero mean, which has variance e2/3. Since the36measurement and truncation noise are uncorrelated, their variances simply add:〈(v − v0)2〉=〈n2〉+e23. (2.4)Given an estimate of the measurement noise, we use this relationship to pick thetruncation precision such that the fractional increase in the noise variance is small(typically 10−3).2.1.2 Software ImplementationBecause of the memory and computational requirements of reprocessing archive-files’worth of data at a time, the compression and archiving stages are run on a dedicatedcomputer node. This added to the complexity of the development, since it required anetwork layer within the archiving pipeline to connect the receiver node and archivingnode. The archiving pipeline was implemented as a series of kotekan stages thatoperate on a given acquisition file:1. read data from disk (receiver node);2. network transfer (receiver node – archiving node);3. transpose the data so time is the fastest-varying index (archiving node);4. truncate the data (archiving node);5. package into archive file format, bitshuffle compress, and write to disk (archiv-ing node).The kotekan pipeline described above operates on a single file at a time. Ahigher-level program was also implemented, in Python, to operate as a daemon thatmonitors the receiver node for new files, spawns kotekan processes on both nodesto process individual files, moves them to their final destination, and deletes theprocessed acquisition files.2.1.3 Compression PerformanceApplied to the main dataset of the CHIME archive files, this scheme typically achievesa compression ratio of ∼ 3 − 4. For comparison, an arbitrary CHIME data file wasselected and recompressed using the LZF2 algorithm. Using the nominal bitshuffle2http://oldhome.schmorp.de/marc/liblzf.html37algorithm the compression ratio for this file is 3.4, compared to 2.2 using LZF. Thissignificant improvement in performance has a considerable impact on the storagerequirements for the > 2 PB of (compressed) CHIME data accumulated so far.2.2 Pulsar GatingObservations of bright point sources have been used for the calibration of the tele-scope’s time-dependent gain and its directional beam response. A limitation of thismethod (for beam calibration in particular) is confusion noise, which originates fromthe presence of a background of numerous dim sources in the vicinity of the brightsource that are also captured by the telescope’s broad beam response. Because it isa systematic effect, confusion noise cannot be mitigated by repeated observations.Measurements of pulsars can be made free of confusion noise, since they providea signal that is chopped in time at a fast cadence (. 1 s). Gating is the process ofsorting the fast-cadence time stream into pulsar on- and off- bins, and taking theirdifference to subtract the background sky. In principle this can provide a very cleanmeasurement of the pulsar as seen through the telescope beam. In practice there arechallenges associated with such a measurement, including the variation in time andfrequency of the pulsar signal, which will be discussed in Section 3.Since the pulsars are switching at a rate much faster than the CHIME data sam-ples are written to disk, gating must be implemented in the real-time pipeline. Iparticipated in adding this feature to the stage that accumulates the intermediatehigh-cadence (31 ms) frames into the final sample length. My main contribution wasto implement the time-of-arrival predictions for the individual pulsar pulses, whichneed to be known precisely in order to determine the on- and off- gates. The inputto the gating code is a set of timing parameters that were derived for each observedpulsar by the CHIME/Pulsar team (CHIME/Pulsar Collaboration et al., 2020). Ialso developed the tools to start, stop and schedule pulsar observations.2.2.1 Gating SchemeThe gating is inserted at a stage in the receiver pipeline called accumulation, in whichsamples at a 31ms cadence are further averaged in time to a final frame length of 5 s.Given N samples per frame (here N ≈ 166), this operation isv =N∑kvk =N∑k(bk + pk + nk) , (2.5)38where we have decomposed the individual samples into background, pulsar, and noiseparts respectively.Within this frame there will be many pulses, each of which will vary in intensitydue to scintillation, but on average typically have a regular profile in time. If weknow this profile, we could weight the samples when binning them into gates inorder to optimise signal-to-noise. We do not necessarily have these profiles for allthe pulsars observed, and in order to avoid additional complexity, we instead use auniform weighting, wherein the pulsar is either considered “on” or “off”. This meansthat for a subset of M < N samples, the pulsar will be “on”, and these samples willbe accumulated separately,von =∑j∈onvj. (2.6)It remains to combine v and von to subtract the background. Assuming that thebackground is not changing significantly over the length of a frame, the average valuein each gate is the same, andvon − MNv =∑j∈on(pj + nj)− MNN∑k(pk + nk) (2.7)=∑j∈on(pj − MNpj)+ (noise), (2.8)where we have used the assumption that the pulsar signal is absent outside the gate.Dividing by the overall factor, we get the final expression for the gated visibilityvg =1N −M(von − MNv)=1N∑j∈onpj + (noise). (2.9)In reality, upstream pipeline stages that operate at even faster cadence excise short-duration RFI from within the 31ms samples, which means that each one should notbe given the same weight. The last equation is easily adapted to account for thisby weighting the M and N by the complement of the fraction of excised data in thegated and ungated accumulations.As to the noise term that was not expanded in the previous equations, its varianceis related to the variance in the ungated visibilities (σ2) byσ2g =MN(1− MN)σ2. (2.10)39This result can be obtained by counting the number of independent samples thatcontribute to a gated visibility, weighting them by (1−M/N) when they are sharedbetween gates andM/N otherwise, and combining them assuming the noise averagesdown according to the central limit theorem. In the accumulation code, the varianceis estimated by accumulating squared differences of even and odd samples at the fast30ms cadence. This procedure works under the assumption that the signal meanis not changing appreciably from one sample to the next, such that the differencescan be attributed to noise. For a pulsar however, the signal level can be fluctuatingsignificantly on this timescale, which may contribute additional power to the noisevariance estimate.2.2.2 Timing ModelThe scheme for gating described in the previous section requires a prediction for thetimes of arrival of individual pulses in order to sort them into the correct gate. Thank-fully, pulsars are extremely regular and the research community who studies them hasdeveloped precise models of their timing. These models depend on parameters thatdescribe both features of the pulsar (mainly its location, period, and derivatives intime) and the location of the observatory where the pulses are to be measured (theDRAO in this case). The current standard in the field is the TEMPO2 software pack-age, which implements these calculations accounting for “the effects of a binary orbitalmotion, the secular motion of the pulsar or binary system, interstellar, Solar systemand ionospheric dispersion, observatory motion (including Earth rotation, precession,nutation, polar motion and orbital motion), tropospheric propagation delay, and grav-itational time-dilation due to binary companions and Solar system bodies” (Edwardset al., 2006). The uncertainty quoted on the timing solutions is ∼1 ns (though thereare probably additional uncertainties associated with the parameters derived from fitsto pulsar observations performed by the CHIME/Pulsar team to calibrate the modelsto the location of the DRAO site). The timing requirements for gating at 30ms areeasily met by the accuracy of the models.The pulsar and observatory parameters input to the model were derived by theCHIME/Pulsar team from their observations using the pulsar tracking beams of theCHIME telescope (CHIME/Pulsar Collaboration et al., 2020). When a gated pulsarobservation is scheduled for CHIME, a timing model that is valid at the time of theobservation is derived using TEMPO2 and the pulsar parameters. This model is passedto the accumulation stage of the pipeline, as described below.402.2.3 Software ImplementationThere are multiple components to the software implementation of the gating:• modifying the accumulation code to support gating as described in Section 2.2.1and handle splitting the data streams into gated and ungated parts;• creating a tool that can generate gating predictions from TEMPO2 output;• adding the capability to toggle this mode on and off and schedule observationsof specific pulsars.The gate timing calculation was implemented as a tool within kotekan. It isinstantiated given a timing model, and provides a function that returns whether thepulsar is on or off at the requested time. This is determined by computing the nearesttime of arrival, shifting it in time to correct for dispersion in the interstellar mediumat that frequency, and padding it by the width of the pulse in time. The dispersionmeasure and width (the width at 10% power is used) are parameters that are passedalong with the timing model.Because observations of pulsars are not continuous and each individual pulsarrequires a different timing, the gating configuration must be allowed to change on thefly, which is supported by kotekan through a REST interface3. A suite of command-line tools were developed in Python to facilitate and schedule these updates. The toolsgenerate a timing model for the requested pulsar and observing time using TEMPO2and communicate the new parameters to the kotekan process at the appropriatetimes.2.3 Real-time MapsCHIME measures visibilities (Section 1.2.2), which are related to the brightness ofthe sky at a given moment by an integral of the sky brightness modulated by theantenna response, for a baseline separation bij:vij ∝∫d2nˆT (nˆ)Ai(φ, nˆ)A∗j(φ, nˆ)e2piinˆ·bij/λ, (2.11)where T (nˆ) is the brightness temperature (i.e. the specific intensity) of the sky in thedirection nˆ and Ai(φ, nˆ) is the beam response of feed i in the direction nˆ at sidereal3The interface designed for the World Wide Web. See https://www.ics.uci.edu/~fielding/pubs/dissertation/rest_arch_style.htm.41time φ (i.e. the rotation angle of the Earth relative to the celestial sphere). The phasefactor in the exponential accounts for the difference in arrival times of the signal ateach feed given by the projection nˆ ·bij/λ. I consider only the total intensity here forsimplicity, but a fuller treatment should also include the response of the instrumentto polarised components of the sky emission.The raw visibilities are difficult to interpret at a glance and they are very numerous(about 16000 baselines are recorded), which makes it difficult to use them directlyto monitor data quality in real time. Much more useful for this purpose is a rollingmap of the sky as it is being observed by the telescope. A quick way to generate amap from CHIME visibilities is to make a “ringmap”, a reconstruction of the sky asit drifts through the CHIME beam in time, described in more detail below.The shape of the CHIME beam is narrow in the EW direction and elongated NSdue to the cylindrical shape of the reflectors. At each time sample, we thus recordvisibilities where the integral in Equation (2.11) is dominated by a narrow regionaround the meridian. For illustration, I will collapse the integral entirely onto themeridian and suppress the time dependence (taking the sky direction nˆ to be onmeridian at a given sidereal time):vij ∝∫d (sin θ)T (θ)Ai(θ)A∗j(θ)e2pii(sin θ) yij/λ, (2.12)where θ is the angle from zenith along the meridian and yij is the projection of thebaseline along the NS direction. This makes it apparent that the visibility is a Fouriertransform of the sky, with the baseline in units of wavelength y/λ conjugate to the skycoordinate sin θ. We can thus obtain an estimate of the sky brightness (modulatedby the beam) at a pixel mk by taking the inverse transform:mk =∑ijvije−2pii sin θkyij/λ=∫d (sin θ)T (θ)∑ijAi(θ)A∗j(θ)e2pii(sin θ−sin θk)yij/λ.(2.13)The summation in the last line is called the synthesized beam. Its form is determinedby the discrete sampling of the aperture by the visibilities and approaches a δ-functionin the limit of a large number of baseline separations. In practice the EW componentsof the baselines are also averaged together, but since they are few, the EW synthesizedbeam is only marginally more localized than the primary cylinder beam. This aspectis neglected in the simplified picture just presented, as the integral was collapsed onto42the meridian.The operation in Equation (2.13) is performed to generate every pixel of a ringmap.It is a discrete fourier transform (DFT), where for a given choice of pixels, the weightsdepend on the wavelength λ. For CHIME, there are 255 NS baseline separations, so agrid of 511 pixels captures all of the information available. The live-ringmap pipelinestage generates these DFT weights for every frequency and performs the matrix mul-tiplication to form an array of pixels at every time sample. A 24h buffer of ringmapsstored in memory can be queried via a REST interface. The ringmaps are displayedin a web-viewer that I initially developed. Only 64 frequencies are directed to thelive-ringmap stage due to limited computational and human resources available togenerate and look at the maps.The live-ringmap pipeline runs as a separate instance of kotekan on the processingnode (see Figure 2.1) in order not to interfere with the data acquisition process onthe receiver node. The 64 frequencies chosen for the maps are communicated overthe network to the processing node, where they are processed into ringmaps, whichcan then be served to clients over the network.An example ringmap is shown in Figure 2.3. The format nicely illustrates theway a cylindrical transit telescope like CHIME observes the sky. Specific features aredescribed in the caption.2.4 ConclusionThe receiver pipeline is a crucial component of the CHIME telescope, which is re-sponsible for the last stages of the data acquisition, performing real-time processingbefore the data are archived. I contributed to its development extensively with soft-ware design and implementation, and in this chapter I have described in detail threepipeline stages for which I led the development: the final compression and packingof the data files; a gated acquisition mode that enables pulsar observations for beamcalibration; and the generation of maps of the sky in real time.The impact of these contributions has been significant for the CHIME project. Todate over 2PB of CHIME data have been acquired, all of it compressed and packedinto archive files by the code that I wrote. The live sky maps are used on a dailybasis to monitor telescope performance and data quality. The gated acquisition modehas been used for regular holographic observations of pulsars (Section 2.2), but alsorepurposed for other calibration tasks – notably a recent experiment that used a droneflying over the telescope broadcasting a known signal to help measure its beam, which43Figure 2.3: Screen-grab of a ringmap as displayed by the live-viewer. The horizontalaxis is time and the vertical axis projected angle along the meridian. Forsources on the meridian at the time, these map to right ascension (RA) anddeclination (Dec), respectively. When sources drift away from the meridianthey remain in the side-lobes of the beam, which can be seen clearly for thebrightest sources as an arc around their transit. In this projection, solid angleis mapped to an area that scales like cos(θ + lat.)−1, which reflects the factthat higher declinations take longer to drift through CHIME’s field of view.Zero-length baselines are excluded from the sums because they are subject toadditive auto-correlation noise, and this missing information results in everyvertical column having zero mean. This is why bright sources produce negativevertical bars. The galaxy is clearly visible, but so is a copy shifted vertically.This is an alias due to the sub-Nyquist sampling of the baseline separations atthis wavelength (the sampling scales like y/λ) which manifests itself as copiesat high and low sin θ (see Appendix C). Finally, the brightest arc that spansover half of the day is the Sun (N.B. daytime data are not presently used forcosmological analyses).44synchronised the data acquisition to the injected signal with the gated mode.45Chapter 3Beam CalibrationThe large dynamic range between foregrounds and the 21 cm signal sets stringentrequirements on the precision of instrumental calibration (Section 1.1.4). For CHIME,calibrating the primary beam poses a particular challenge as a result of the cylindricaltransit telescope design. Multi-path effects and coupling between feeds are significantwithin the cylinders, which leads to a complicated frequency response (a generalsummary of the CHIME primary beam was given in Section 1.3.1). Because thetelescope has no moving parts and a very large field of view, it is not possible to scanit across celestial calibration sources to measure the beam the way steerable dishescan be.In this chapter I will describe work that I did towards characterising the CHIMEprimary beam. Note that this does not nearly constitute a complete summary of beamcalibration efforts for CHIME, which are being carried out by a number of people,using a variety of approaches. Some of these have been published: a model derivedby fitting visibilities to catalogues of radio point sources is described in CHIMECollaboration et al. (2022b) and also summarised in Section 4.5.4, and measurementsof the beam using the Sun as a calibrator were presented in Amiri et al. (2022).3.1 Holography Data ProcessingHolography is an established technique for making accurate measurements of theamplitude and phase of antenna beams at radio frequencies (e.g. Bennett et al. (1976);Scott & Ryle (1977); Baars (2007)). CHIME uses this technique by tracking a celestialsource with a nearby moving telescope while the source transits through the stationaryCHIME beam. The correlation between the signals from each stationary feed andthe tracking reference telescope picks out only what is common between the (large)46Figure 3.1: Photo of the DRAO site, with CHIME visible on the left and the26m John A. Galt telescope on the right (in stow next to the mainbuilding). The Galt telescope is separated from CHIME by 230m, nearlydue East.CHIME primary beam and the more compact tracking beam, tracing the responseof CHIME along the path of the source. For CHIME holography, the John A. Galt26m telescope, located 230m east of CHIME, is used as the tracking system. TheGalt telescope and its position relative to CHIME are shown in Figure 3.1. For theseobservations, a 400–800MHz dual-polarization modified CHIME receiver is mountedon the Galt telescope (Berger et al., 2016). The resulting cross-correlations yieldCHIME’s co-polar and cross-polar far-field beam response (amplitude and phase) perfeed, per frequency, along a track in hour angle at the declination of each observedsource. A holographic observation of Cyg A is plotted in Figure 3.2, where the detailedbeam side-lobe information and phase measurements are highlighted.Holography observations with the Galt telescope are the most accurate beammeasurements available to CHIME, sensitive to amplitude and phase for every feed.Their limitation is sparseness within the large CHIME primary beam, since they onlyprovide a measurement along the east-west trajectory that is followed by individualpoint sources as they transit over the telescope. These tracks are illustrated for a fewof the brighter sources in Figure 3.3.To date, about 1800 observations of a total of 23 different sources have beenacquired. In its raw form, these data are not easily interpretable as a beam mea-surement. To make these measurements more accessible to researchers (including4715 10 5 0 5 10 15103102101100NormalizedamplitudeCo-pol amplitude15 10 5 0 5 10 15103102101100NormalizedamplitudeCross-pol amplitude4 2 0 2 40.00.20.40.60.81.0NormalizedamplitudeCo-pol amplitude (zoomed)3 2 1 0 1 2 30.0000.0250.0500.0750.1000.1250.150NormalizedamplitudeCross-pol amplitude (zoomed)4 2 0 2 42⇡ 32⇡⇡⇡/20⇡/2⇡32⇡2⇡Phase(radians)Co-pol phase (referenced to transit)3 2 1 0 1 2 32⇡ 32⇡⇡⇡/20⇡/2⇡32⇡2⇡Phase(radians)Cross-pol phase (referenced to transit)Degrees on SkyCylinder ACylinder BCylinder CCylinder DFigure 3.2: The CHIME X-polarized beam response at 717MHz from the holographicmeasurements of a Cyg A transit, taken on 2018 Sep 28. Each panel shows themedian response taken over all feeds within a cylinder, normalized such that the co-polar response is 1+0i at transit. (Upper left): The median co-polar amplitude of allfeeds per cylinder, normalized by the Gaussian-fit peak height, over the full extentof the observation, converted to degrees on sky, HA· cos (δCygA). (Upper right): Themedian cross-polar amplitude from the product of a CHIME feed with the oppositepolarization on the Galt receiver. The data have been scaled by the same factor asapplied to the co-polar response, so the curves give an indication of the level of cross-polarization in the beam. (Middle left): Same as the upper left panel, but zoomedto a smaller hour angle range and plotted on a linear scale. (Middle right) Same asthe upper right panel, but zoomed to a smaller hour angle range and plotted on alinear scale. (Lower left): The median co-polar phase as a function of scaled hourangle, taken over all feeds in a cylinder. (Lower right): Same as the lower left forthe cross-polar phase. The phase difference between cylinders, after accounting forphase wrap, is only large near the first zero-crossing of the field. The gray bands inthe amplitude plots indicate the standard deviation over all the Cyg A holographytracks of Cylinder A’s median feed response (Cylinder A is representative). Figureand caption reproduced from CHIME Collaboration et al. (2022a).48Figure 3.3: (Left) Trajectories through the CHIME beam of six of the brightestholography sources. The coordinates are direction cosines, with North atthe top. The horizon is drawn in black. The colour of the tracks repre-sents the amplitude of the beam measurement at that location. (Right)Corresponding beam amplitude measurements (phase is observed butnot shown), for a single feed, at frequency 717MHz in the Y polarisa-tion. This figure was generated by Alex Reda for a CHIME holographypaper in preparation.myself), I have implemented a data processing pipeline to prepare them for anal-ysis. The pipeline includes steps that correct for the interferometric phase due tothe Galt-CHIME baseline, correct for decorrelation effects due to the large delay ofthis baseline, interpolate the measurements onto a common grid, and package theobservations into a format that is well-suited to further analysis and processing inthe analysis pipeline framework. A brief description of the CHIME analysis pipeline’ssoftware framework will be provided in a later chapter (Section 4.5), while here I willfocus on the function of each stage of the holography processing.493.1.1 Fringestopping and RegriddingTo facilitate the comparison of beam measurements from multiple observations andsources, the time-ordered holography data are resampled onto a common grid in hourangle (HA) relative to the transit time of the source on the telescope’s meridian.This regridding is implemented as an inverse Lanczos interpolation, explained inAppendix A, and takes the data from a 5 s cadence to approximately 20 s. The chosengrid spans ±30◦ of HA, which is broad enough to cover most holography observations.Prior to downsampling in HA, a phase correction is applied to undo the fringingin the data produced by the changing delay between the CHIME and Galt feeds asthe source moves overhead (“fringestopping”). Without this correction, averagingadjacent time frames would wash out the signal. The phase correction is2pinˆ · dkλ, (3.1)the projection of the direction to the source nˆ onto the baseline between the Galtand the kth CHIME feed in units of wavelength λ. For a source at 0◦ declination onthe meridian, this corresponds to >4.5 rad of phase over a period of 20 s at 600MHz,enough to average out the signal. On the other hand, the full width at half maximumof the beam at that frequency is >3◦, an amount of rotation it takes the Earth ∼ 720 sto make, so the 20 s sampling is sufficient to capture the variation of the signal intime.In order to apply this correction, the position and geometry of the Galt telescopedo not need to be known exactly. Any error in the baseline distances dk will introducea phase error in the beammeasurement, but as long as it is small it will not degrade thesignal and it can be identified and corrected in later stages of processing if necessary.3.1.2 Decorrelation CorrectionCHIME contains 60m of coaxial cable, while the Galt telescope signal traverses 405mof cable. The result is that the signal from the Galt telescope is delayed relative toCHIME’s signal by 1.32 µs. Together with the geometric delay due to the projectionof the source position onto the CHIME-Galt baseline, this adds up to a delay that is asignificant fraction of the integration length in the frequency transform – 2.56 µs. Thismeans that the signals are not entirely overlapping for the CHIME and Galt inputsand this creates a decorrelation effect when the product of the two is taken, themagnitude of which depends on the delay. Since the geometric delay is changing overthe course of an observation, this modulates the beam measurement as a function of500.0 0.5 1.0 1.5 2.0 2.5delay ( s)0.00.20.40.60.81.0decorrelationHolographyCHIMEFigure 3.4: Magnitude of decorrelation in the PFB as a function of delay. Thecoloured regions show approximately the range of delays that are relevantto the baselines internal to CHIME, and to those that involve holographywith the Galt telescope. Although decorrelation is mostly negligible forCHIME baselines, it is significant and importantly has a slope over therange of delays in holographic baselines.HA, an effect that is observable in the holography data. We calculate the magnitudeof this decorrelation given the details of the frequency transform and the measureddelay between the receivers to correct it at this stage of the processing pipeline. Themagnitude of the suppression due to this effect is shown in Figure 3.4.3.1.3 Data PackingHolographic visibilities are products of a CHIME input and one of the two Galt po-larisations. The raw holography observations are packed such that all visibilities arestored along a single axis, which makes it difficult to identify the relevant measure-ments in an analysis that is interested in specific feeds or polarisations. The laststep of the holography processing is to expand and reorder the visibility axis intotwo dimensions: a co/cross-polar axis; and an ordered list of inputs. This makes itmuch easier to access beam measurements for individual feeds. I also created datastructures specific to the holography observations that provide convenience methodsfor, e.g., the location of the observed source on the sky.513.2 Analysis of Unexpected Slopes in HolographyAmplitudesIn this section I discuss an unexpected feature in the holography observations. Iwill describe its defining characteristics and explore possible explanations, but con-clude that none can consistently explain it. This account may be useful for futureinvestigations.It has been noted that the amplitude of holographic visibilities at the time whenthe source is on meridian scale linearly with north-south baseline length, with a slopethat changes sign for sources on either side of zenith. This effect is not explained bydecorrelation due to the significant delay between the Galt and CHIME signals, norby the effect of finite source extents being resolved by the long baselines. Figure 3.5shows the effect in unprocessed and uncalibrated holography data.3.2.1 Characterisation of the SlopesHolographic observations of the three bright sources for which calibration beam ratiosare available – Cyg A, Cas A, and Tau A – are considered in this analysis: observationsnumbered respectively 3541, 3816, and 3768.CalibrationThe following calibration steps are applied to each observation.• As described in the previous section, a decorrelation correction is applied aspart of the holography processing pipeline. The correction depends on baselinelength and source declination but not on frequency.• Gains derived in the real-time pipeline for the CHIME feeds are applied to thedata. These correct for feed-to-feed gain variation and reference the amplitudeas a function of frequency to the beam response in the direction of Cyg A onmeridian. The gains are defined per-frequency and per-feed.• A beam ratio that was derived as part of the real-time pipeline calibrationprocedure is applied to correct for the different response of the CHIME primarybeam in the direction of the different sources. It is defined per-frequency andper-feed, for each source. By definition the beam ratio for Cyg A is identically1.52Figure 3.5: Amplitude of a holographic observation of Tau A on meridian,taken on 2018/10/24, at 714 MHz. The horizontal axis is feed positionfrom south to north within each cylinder, ordered from A to D. Onlyinputs for the X polarisation are plotted. Note the decreasing trend asa function of feed within a cylinder.• The catalogue flux (from Perley & Butler (2017)) for each source is dividedout to remove the intrinsic amplitude variation with frequency. This correctiondepends only on frequency and source.Assuming all of these corrections are accurate, they should remove all knowninstrumental effects due only to the CHIME signal chain. Any remaining modulationin the visibilities as a function of feed and frequency can therefore be attributed toeffects particular to the holography. Figure 3.6 shows some on-meridian data post-calibration.Slope FitsFor each source and for each cylinder independently, a linear model was fit to theamplitudes as a function of NS feed position. The resulting slopes are shown in53Figure 3.6: Amplitude on meridian of calibrated data for a holographic obser-vation of Tau A at 714 MHz, as a function of X-pol input. The colouredlines are fit to the data within each cylinder independently. Within eachcylinder, feed index increases from south to north. Note that points ly-ing on the origin indicate missing data and were excluded from the fit.Figures 3.7–3.9. The trend with frequency is fairly smooth, with the amplitude ofthe slope growing as frequency increases. All four cylinders are largely consistentfor a given source. Between sources, it should be noted that the sign of the slopeschanges for sources north and south of zenith (compare Cyg A and Tau A to Cas A).Otherwise, with only three sources to compare it is hard to identify a trend.SVD AnalysisA singular value decomposition (SVD) was evaluated on the feed-frequency array of(complex) amplitudes for each source. The first two modes are shown in Figures 3.10and 3.11 for Tau A. The subsequent modes do not appear meaningful, but rathercapture RFI- or outlier-feed-like features. This is demonstrated in Figure 3.13 wherethe residuals of the data and the two-mode reconstruction are plotted. The SVD54Figure 3.7: Slopes at each frequency for linear fits of amplitude on meridianas a function of NS feed position within a cylinder. Each cylinder andfrequency were fit independently. Note that points of zero slope corre-spond to frequencies that were masked due to contamination from RFI.modes for the other sources appear qualitatively similar, and are not plotted for thesake of brevity. The expansion of the first two modes individually into the feed-frequency space is shown in Figure 3.12Note that a residual delay of ∼ −90 ns was discovered and corrected for whilederiving the SVD modes. It has been removed prior to making any of the plots inthis section. This is likely an error in the delay to the Galt feed that is used in theprocessing pipeline.A possible interpretation for the two modes is that one is needed to describe theslope within cylinders, and a second one to account for the steps in the total baselinelengths that occur between cylinders, which could lead to finite source extent effects(as will be discussed later). This interpretation is supported by the fact that one modeis sufficient to capture most of the meaningful structure if the SVD is performed on55Figure 3.8: Slopes at each frequency for linear fits of amplitude on meridianas a function of NS feed position within a cylinder. Each cylinder andfrequency were fit independently. Note that points of zero slope corre-spond to frequencies that were masked due to contamination from RFI.each cylinder independently, i.e. they are degenerate within cylinders. The result ofsuch an SVD is presented in Figure 3.14 for Tau A.3.2.2 Expected EffectsIn the next few paragraphs, I explore two effects that are expected to be present inthe holography and could have produced a feature similar to the one under consider-ation here. I then show how these explanations are not consistent with the observedcharacteristics.DecorrelationThe decorrelation due to the large delays between the CHIME and Galt signals thatwas discussed earlier has some resemblance to this observed effect: as the baseline56Figure 3.9: Slopes at each frequency for linear fits of amplitude on meridianas a function of NS feed position within a cylinder. Each cylinder andfrequency were fit independently. Note that points of zero slope corre-spond to frequencies that were masked due to contamination from RFI.length increases, the delay and the degree of decorrelation change. However, the slopethat is expected due to decorrelation has the opposite sign to what is observed in thedata. The suppression from decorrelation calculated for a few sources on meridianis plotted in Figure 3.15. One might wonder if the decorrelation correction appliedin the pipeline is wrong and actually introducing these slopes, but this is not thecase: the slopes are present in the uncalibrated data, as shown in Figure 3.5, and thecorrection only amplifies them slightly (note that the decorrelation was first identifiedby its modulating effect over the course of a transit, which is successfully accountedfor by this correction).57Figure 3.10: First SVD mode for Tau A, with the frequency-dependent modein the top panel and the feed mode beneath. Within each cylinder, feedindex increases from south to north. The feed dependence is charac-terised by slopes within cylinders and steps between them, while thefrequency dependence is fairly smooth.Extended Source EffectsThe large baselines between the Galt telescope and CHIME have the potential toresolve out the angular extent of bright sources that are usually assumed to be point-like. There is evidence for this effect in the steps observed in the amplitudes forholography on meridian (Figure 3.6), which occur between cylinders where the EWbaselines change in increments of ∼ 20m. This could potentially also explain theobserved slopes within cylinders if the baselines were monotonically increasing inlength as feed positions move north within a cylinder. This is not the case, as theGalt telescope is located approximately in line with CHIME along the NS direction,so the baseline length turns over in the middle of the cylinder as shown in Figure 3.16.This effect could still be important for holography in the sidelobes (and is cur-58Figure 3.11: Second SVD mode for Tau A, with the frequency-dependent modein the top panel and the feed mode beneath. Note that the magnitudeis not shown in the bottom panel. Within each cylinder, feed indexincreases from south to north. The singular value relative to the firstmode is given in the title. Subsequent modes do not appear meaningfulfor this analysis.rently unaccounted for) because the projected length of the baselines changes as thesource moves across the sky, which will modulate the amplitude of the suppression.Figure 3.6 demonstrates that changes of ∼ 20m in the baseline length modify theamplitude by a large fraction for Tau A, and a commensurate change in the projectedbaseline length will occur over the course of the transit. Correcting for this effectwould require a model for the angular extent of each source.In order to obtain an order of magnitude estimate for this effect, the simplestmodel for the source extent we can consider is a Gaussian centred on the sourceposition with a width σ:T (nˆ) ∝ e|nˆ−nˆ0|2/2σ2 . (3.2)59Figure 3.12: Magnitude of the first and second modes of the Tau A SVD ex-panded into feed-frequency space. Within each cylinder, feed indexincreases from south to north.Assuming neither the Galt nor CHIME beam varies significantly over the extent ofthe source, a holographic visibility for a baseline u can be written asv ∝∫d2nˆT (nˆ)e2piiu·nˆ. (3.3)Since the extent of the source is very small, we can use the small angle approximationto evaluate the integral in Cartesian coordinates, centred on the position of the source,nˆ = nˆ0 + x:v ∝∫d2xe|x|2/2σ2e2piiu′·x ∝ e−2pi2|u′|2σ2 , (3.4)where we have used the result for the Fourier transform of a 2-d Gaussian, u′ is thebaseline projected onto the tangent plane at the location of the source (on meridian):u = (u, v)→ u′ = (u, v cos θ), θ being the angle of the source away from zenith. Theresult of this calculation for a few different source extents are plotted in Figure 3.16.The Figure shows that even a moderate source extent can introduce an O(1) ef-fect to the measurements, which deserves further investigation and may need to beaccounted for in future analyses aiming to interpret the holography as beam mea-surements for certain sources.60Figure 3.13: Residuals from the Tau A amplitude data to the two-mode re-construction. The horizontal bands correspond to known RFI contam-inated regions of the spectrum, and other bright features also resemblesporadic RFI. Within each cylinder, feed index increases from south tonorth. No features coherent within cylinders are apparent.3.2.3 Other ChecksThe following potential contributors to this effect were investigated, but none werefound to be significant.• CHIME auto-correlations: No slopes were found in the auto-correlation ampli-tudes on meridian for a recent transit of Tau A.• Digital gains: No such feature was found in the digital gains that are appliedto scale the data prior to multiplication and accumulation in the correlator.• CHIME gains: No such feature was identified in the CHIME real-time calibra-61Figure 3.14: Single SVD mode for an independent decomposition of each cylin-der for Tau A. The solid/dashed lines are real/imaginary parts. Feedindex increases from south to north.tion gains derived daily from bright point source transits.• Geometric delay and decorrelation: The total delay as a function of feed posi-tions was measured in a delay transform of an unprocessed (i.e. not fringestopped,no decorrelation correction) holography observation of Tau A and the resultmatches what is assumed in the decorrelation calculations.• Asymmetry in HA was checked for by folding the transits around HA=0 andno convincing trend was found. The decorrelation effect is asymmetric in HAand it appears to be corrected for adequately.• A rotation of the cylinders away from true North would introduce a linear shiftin the centroid of the beam as a function of NS position. This was ruled out asan explanation by fitting a gaussian shape to the main lobe at every feed. Thecentroid was not found to have a linear trend with feed.62Figure 3.15: Top: Expected delay for different sources on meridian as a func-tion of feed position. Within each cylinder, feed index increases fromsouth to north. The total delay is composed of a fixed delay to theGalt feed and a geometric delay determined by the projection of thesource position onto the NS baseline. Bottom: Corresponding resultof decorrelation calculation for the suppression, referenced to the firstfeed. Note that this effect produces slopes in the amplitude, but thesehave the opposite sign to those observed in the data.• Misalignment of the CHIME array plane relative to astronomical zenith pro-ducing a feed-dependent vertical offset was ruled out by looking for this featurein the delay transform of CHIME gains.Although this analysis has failed to identify the source of the observed featurein the holography measurements, it has ruled out a number of candidates. Furtherinvestigation will be required to resolve this question, and the results described hereshould be used as a starting point for those efforts.63Figure 3.16: Top: Suppression factor due to resolving an extended source as afunction of feed position, calculated based on a Gaussian model for thesource extent at a few different widths (σ in the legend), and referencedto the last feed. Within each cylinder, feed (or baseline) index increasesfrom South to North. Bottom: The magnitude of the projected base-line length as a function of feed position. This is the quadrature sumof the EW and NS components, projected onto the tangent plane atthe location of Tau A. Within a cylinder, the length turns over becausethe Galt telescope is located inside the range of CHIME NS positions.This effect also produces a modulation of the measured amplitude, butit does not correspond to the linear trend that characterises the obser-vations.643.3 Pulsar HolographySince pulsar gating was implemented in the real-time pipeline (Section 2.2), holog-raphy observations between the Galt and CHIME telescopes have been proceedingregularly. To date, over 300 observations have been recorded, for over 20 pulsars.Deriving beam measurements from these observations has so far proved challeng-ing, partly due to the large variability of the pulsar signal. Most pulsars exhibitsome level of scintillation, variability of their observed brightness in both time andfrequency due to scattering in a turbulent interstellar medium surrounding the source(Narayan, 1992). Pulsar scintillation can vary on a timescale of minutes and frequencyseparations of a few MHz. See Figure 3.17 for an example. To measure the variationdue to the CHIME beam, the intrinsic variability of the source must be known in-dependently. Thankfully, this information is recorded in the auto-correlations of theGalt telescope as it tracks the source in holography.However, scintillation also greatly suppresses the already faint pulsar signal atcertain times and frequencies in any given observation, so multiple observations ofthe same source must be combined in order to fill in the track through the beam.This introduces another potentially confounding source of variability that must becorrected for: instrumental gains. The amplitude and phase response of the telescopeslowly drifts over the course of a day, so observations from different days cannot becompared before a calibration is performed. For CHIME, this is done daily (Sec-tion 1.3.2) and the gains are available for use in holography, but gains for the Galttelescope needed to be derived as part of the analysis.One reason to pursue gated pulsar holography that has already been mentionedis its potential to mitigate confusion noise, by subtracting the background measuredbetween pulses, but another is the possibility of calibrating the polarisation responseof the telescope. Strongly polarised celestial point sources are rare – of the bright radiopoint sources, Tau A is the most polarised with ∼ 7% linear polarisation at 22GHz(Weiland et al., 2011) – but pulsars often have even higher amounts of polarisation.For example, the pulsar B1929+10 is nearly 100% polarised at 610MHz and manyothers have > 50% polarisation (Stairs et al., 1999).In this section I will detail an analysis method that I developed to measure theCHIME beam using gated pulsar holography, including calibrating the CHIME andGalt polarisation response. Due to lack of time and the prioritisation of other projects,I have not pursued this analysis on actual data very far, apart from a preliminaryattempt that shows some promise.65Figure 3.17: Dynamic spectrum from a gated observation of the pulsarB0329+54 using the 26m Galt telescope on August 8th 2019. Onlythe Galt autocorrelation data are shown here (i.e. no CHIME inputsare involved). The gated pulsar signal exhibits a scintillation patternthat varies quickly in time and frequency. Horizontal bands that areempty or corrupted correspond to masked and unmasked RFI.3.3.1 Disentangling the CHIME and Galt GainsA measured holography visibility is the product of a CHIME input and one of theGalt polarisation inputs. Both of these include a complex gain term that accounts forthe amplitude and phase imparted on the signal by the respective analog electronics.Schematically,vij = (gCHIMEi si)(gGaltj sj)∗ , (3.5)where the g terms are the gain and s signal, and the indices represent the specificCHIME and Galt inputs being correlated. The CHIME gain is calibrated daily, so inwhat follows we will assume that it has been perfectly corrected, and only the Galtgain is unknown.For a given CHIME feed, we will write the gated visibility during a pulsar holog-66raphy observation (after gain calibration) asvj(φk) = g∗j b(φk)pj(φk) + njk , (3.6)and the corresponding Galt autocorrelation asaj(φk) = gjg∗jpj(φk) + najk , (3.7)where pj(φk) is the pulsar signal in the jth observation at hour angle φk, and gj is theGalt gain. The goal is to recover the beam profile b(φk), but this is made more difficultby its degeneracy with the gain term in any single observation. In these expressions,we have neglected the extent of the Galt beam, because it is much smaller than theCHIME primary beam. We also have not included polarisation – implicitly only co-polar visibilities are considered here – but we will extend the procedure developed inthis section with polarisation later on by building on each successive section.A possible way to break this degeneracy is to use the fact that the gain andbeam factors vary along different directions: the beam profile is unchanged from oneobservation to the next, and the gain factor is (assumed to be) constant throughout agiven observation. We can capture the variation in the gain by a nuisance parameterwhen fitting for the beam. Consider a model for the visibilityvˆjk = αjβkajk , (3.8)where the φ dependence has been expressed as a second index, k, and the parameterscorrespond toαj = g−1j , βk = b(φk) . (3.9)ForN observations composed ofM samples each, we haveM+N parameters andMNsamples, so we should have enough information to constrain them all independently.The model in Equation (3.8) is non-linear in the parameters, and can be fit using astandard algorithm like Levenberg-Marcquardt.Least Squares SolutionWe can optimise the parameters in Equation (3.8) by minimisingχ2 =∑jk|vjk − αjβkajk|2σ2jk. (3.10)67However, because both vjk and ajk are corrupted by noise, the noise weights will inprinciple depend on the true value of the parameters:σ2jk = σ2v;jk +∣∣∣αˆjβˆk∣∣∣2 σ2a;jk , (3.11)where σ2v;jk =〈njkn∗jk〉and σ2a;jk =〈najkna∗jk〉are the noise variances of the visibilitiesand Galt auto-correlations, respectively. Accounting for this can be done by updatingthe weights at every iteration of the fitting routine.In the previous expressions, the beam profile is parameterised at every point inhour angle, but we do not expect that many degrees of freedom in the actual beamshape. We can reduce the parameter space and improve the convergence of thefit while also decreasing the risk of overfitting by choosing an appropriate basis todecompose the fitted profile. Fourier modes are an obvious choice. We expand thebeam profile asβm = b(φm) =∑lγl e2piilm/M , l ∈ Z, −L/2 < l < L/2 , (3.12)where the γl are coefficients to be optimised, and there are L of them, with L < M .A suitable choice for L needs to be made based on our expectation for the shape ofthe beam.Another thing to note is that the model in Equation (3.8) is completely insensitiveto an overall multiplicative factor shared by all parameters. In other words, thetransformationαj → cαj, βk → c−1βk (3.13)leaves αjβk unchanged. This introduces a degeneracy direction in the fit, which canbe resolved by fixing one of the gain parameters to 1.The result of this procedure applied to measurements of the brightest pulsar ob-served in CHIME holography are shown in Figure 3.18. The derived beam shape is ingood agreement with independent measurements from regular holography of nearbysources. This analysis has not yet been attempted on other pulsars, for which thesignal-to-noise is lower. I also found that RFI in some of the observations would derailthe fits, so I had to remove those. More targeted flagging of RFI would be necessaryto make full use of this dataset.68Figure 3.18: Amplitude and phase of the beam profile fit to 25 gated hologra-phy observations of the pulsar B0329+54 using the method describedin this section, for a single CHIME input (which is representative of theothers). Also plotted are independent measurements of the beam forthat input from regular holography of sources with nearby declinations(indicated in the legend). Note that the smoothness of the pulsar mea-surement is imposed by the restriction of the fitted profile to 16 Fouriermodes.3.3.2 Separate Gain FitThe procedure outlined in the previous section requires a non-linear search over aparameter space with M + N dimensions per input, which may not scale very wellcomputationally if we want to fit for all inputs simultaneously (in order to obtain aconsistent solution for the gains). Since the Galt gains are common across visibilitiesfor all CHIME inputs, we could instead consider fitting for the gains together withan average beam profile, followed by a per-input fit. Consider the average visibilityv¯jk =∑iwiv(i)jk , (3.14)69where the i index runs over input and the weights wi are inverse variance, from thefast-cadence estimate in the real-time pipeline. Looking back at Equations (3.6) and(3.7), we see that ajk is common to all inputs, and we recover the model formvˆjk = αjβ′kajk, (3.15)where nowβ′k ∼∑iwibi(φk) (3.16)is estimating the average beam, and we can use the procedure in the previous sectionto fit for it and the gain parameters αj ∝ g−1j .Given the gain solution, Equation (3.8) becomes linear in the parameters βk. Wecan construct an estimator for the parameters by taking the correlation between vjkand the gain-corrected ajk:βk =∑j wjvjkα∗jajk∑j wj∣∣α∗jajk∣∣2 . (3.17)If we take the weights to be inverse variance of the auto-correlations, σ−2jk , and expandvjk and ajk, this can be rewritten asβk = b(φk)∑j |gj|2 p2jk/σ2jk∑j(1 + |gj|2 p2jk/σ2jk) , (3.18)which has a form reminiscent of the Wiener filter (Wiener, 1949).3.3.3 Including PolarisationTo add polarisation to the model in Equation (3.8), we start with an expression forthe two polarisation inputs of a given feed in terms of the electric field on the skyεa(nˆ):Fq =∫d2nˆJaq (nˆ, φ)εa(nˆ)e2piinˆ·u , (3.19)where Jqa(nˆ, φ) is the Jones matrix for the CHIME feed at hour angle φ. In theJones calculus treatment of polarised radiation, a Jones matrix is an operator thatrepresents the effect of an optical element on a polarisation vector (van Straten, 2004).Thus Jqa projects the sky polarisation vector εa onto the measured polarisation basisFq. A visibility is thenVpq =〈FpF∗q〉=∫d2nˆJap (nˆ, φ)〈εa(nˆ)ε∗b(nˆ)〉J b∗q (nˆ, φ)e2piinˆ·∆u . (3.20)70For holography, the second Jones matrix is the Galt telescope’s, Gap, and the integralcollapses approximately onto the location of the tracked source due to the compactGalt beam extent. In matrix notation, a fringestopped holographic visibility can beexpressed asV = J〈εε†〉G† = JpG† , (3.21)where p = 〈εε†〉 is the correlation matrix of the electric field on the sky. The Galtautocorrelations areA = GpG†. (3.22)With these definitions, our model becomesVˆjk = JkG−1jk Ajk , (3.23)where the capital letters indicate 2×2 polarisation matrices, and the indices representsamples in hour angle and holographic observation, as before. If we take our polarisa-tion basis to be celestial unit vectors, then, neglecting leakage between polarisations,the equatorially mounted Galt feed tracks the rotation of the Earth and G−1jk ≈ G−1jdoes not depend on hour angle and is approximately diagonal. We haveVˆjk = JkG−1j Ajk , (3.24)which is very similar to Equation (3.8), except that the terms are matrices. We shouldbe able to apply the same non-linear least squares procedure to fit for the elementsof the matrix Jk, parametrised by a compressed basis of Fourier modes, and the gainfactors in the diagonal elements of G−1j .To adapt the estimator in Equation (3.17) to the matrix problem, we assume thatwe have solved for the gains using the non-linear method on the average over inputsto obtain αj = G−1j , now a matrix. An estimate of the Jones matrix is thenJˆk =∑jVjkwj (αjAjk)†[∑j(αjAjk)wj (αjAjk)†]−1, (3.25)where the wj are once again inverse-noise-variance weights.3.3.4 Calibrating the Galt Jones MatrixIn the previous section, we assumed that the Galt Jones matrix was diagonal, whichwe might expect to be a good approximation but have not verified. Given that many71pulsars produce polarised emission, we can also use these observations to determinethe polarised response of the Galt feed.Relaxing the assumption of diagonality, we write the Jones matrix asGj = gjG0, (3.26)where gj is a diagonal gain matrix that varies from day to day, and G0 is assumed tobe fixed and encodes the polarisation response of the Galt feed. With this change, thefit procedure described in the previous section is the same, except that we identifythe parameters αj = g−1j , and also G0 must be included in the expected value of theestimator (where previously we assumed G0 = I):Jˆk → JkG−10 . (3.27)Measuring G0 is therefore necessary in order to recover the CHIME Jones matrix onits own.The model for the holographic visibilities and auto-correlations wasVjk = JkpjkG†0g†j ; Ajk = gjG0pjkG†0g†j , (3.28)where pjk is the pulsar signal matrix. In order to determine G0, we will focus onthe autos, where the signal to noise is higher over the entirety of the observation,since there is no modulation by the CHIME beam. After fitting for the gain param-eters αj, we can apply the correction and sum over samples to obtain an averageautocorrelation:A¯ =∑jk wjkαjAjkα†j∑jk wjk= G0p¯G†0, (3.29)where we have defined p¯ to be the weighted average of the pulsar signal.Now it remains to isolate G0. A similar problem was worked out in a PhD thesisby Gregory Paciga (Paciga, 2013), and I used this as a starting point to develop thefollowing procedure. We diagonalise the signal matrixp¯ = PΛP † (3.30)and the measured density matrixA¯ = UΛ′U † = G0PΛP †G†0, (3.31)72for which the eigenvalues are the same up to a multiplicative factor Λ′ = λΛ. We canidentifyGˆ =√λUP † (3.32)as a solution that satisfies Equation (3.29). Such a solution, however, is not unique.This can seen by considering the expansion of the signal matrix in the basis of Paulimatrices, which span the space of 2× 2 matrices – a useful formalism for expressingpolarisation detailed in van Straten (2004):p¯ = S0 + S · σ = T +Qσ1 + Uσ2 + V σ3. (3.33)In this basis, the product of two matrices can be expanded as (van Straten, 2004)AB = (a σ0 + a · σ)(b σ0 + b · σ) = (ab+ a · b)σ0 + (ab+ ba+ ia× b) · σ. (3.34)The matrices A and B therefore commute when a× b = 0, i.e. a and b are colinear.Using this identity, we can insert a unitary matrix R = r0+r ·σ into Equation (3.29)asA¯ = G0Rp¯R†G†0, (3.35)without changing the value of A¯ if R is such that r × S = 0, i.e. commutes with p¯.This demonstrates that the solution Gˆ is not unique.In order to fix this degeneracy, we need a second polarised source with a polari-sation vector that has an orthogonal component to the first: S × S2 6= 0. With thisconstraint, a solution Gˆ that solves for both sources can no longer be modified by asingle matrix R without disagreeing with one of the measurements.Another potential ambiguity lies in the normalisation of the Jones matrix. Anoverall phase factor is unobservable as it gets cancelled out in Equation (3.29), butthe overall amplitude needs to be set. In this context, we are not interested in thisamplitude (it can be absorbed into the gain), so we remove these two degrees offreedom by normalising by the determinant:|G| = 1. (3.36)Given two sources with known polarisations p¯1 and p¯2, and their measured auto-correlations A¯1 and A¯2, the following procedure can be used to estimate the Jonesmatrix of the Galt telescope. We diagonalise each matrix and combine the eigenvec-73tors to form two incomplete Jones matrix solutionsG1 = U1P†1∣∣∣U1P †1 ∣∣∣−1 ,G2 = U2P†2∣∣∣U2P †2 ∣∣∣−1 . (3.37)Each of these is ambiguous up to a unitary matrix generated by the polarisationvectorRi(φi) = cosφiσ0 + i sinφiSˆi · σ, (3.38)where Sˆi is the normalised vector of Stokes polarisation components of p¯i and φi is afree parameter. The parameters {φ1, φ2} must be optimised in order to satisfyG1R1(φ1) = G2R2(φ2) = Gˆ, (3.39)and obtain a consistent solution.The procedure developed in this section has not yet been put to the test on gatedholographic observations of pulsars using CHIME and the Galt telescope. Identifyingthe most promising measurements, vetting them for RFI or other contamination, andimplementing the analysis in the CHIME pipeline framework is ongoing and expectedto be continued as future work.3.4 ConclusionHolography is a valuable tool for calibrating the CHIME beam as it provides a highsignal-to-noise measurement of the beam in both amplitude and phase, albeit onlyalong the one-dimensional tracks of transiting bright celestial sources. In this chapterI described the pipeline that I developed to process these observations into beammeasurements that can more easily be used for analysis. The next stages of thispipeline are actively being developed by other members of the collaboration, includingcombining multiple observations of the same sources to improve the signal-to-noise ofthe beam measurement. This is important for the fainter sources that produce noisyindividual measurements, and will greatly extend the reach of holography across thelarge CHIME primary beam. A paper describing holography with CHIME, includingthe full pipeline, is currently in preparation.I also outlined an analysis framework that I developed to produce measurementsof the beam from gated pulsar holography. Pulsars have the potential to provideunique constraints on the CHIME beam because they are often strongly polarised,74and few avenues exist to calibrate the polarised response of the telescope. AlthoughI have not carried out the analysis of the gated holography using this method be-yond a preliminary test, I have shown that in principle observations of at least twopolarised sources in holography can enable the measurement of the Galt telescope’sJones matrix, and in turn the CHIME polarised beam. Understanding the propertiesof the polarised emission from the reference pulsars will be essential for this methodto be successful, and external measurements from pulsar surveys will need to be iden-tified and integrated into the analysis. Another obstacle to extending this analysisto more data is the presence of RFI that must be identified and flagged. A task thatdoes this has been implemented for the ungated holography and could be extendedto the gated holography in the future. The signal-to-noise of these measurements isrelatively low, so combining multiple observations will be necessary. Although theseand likely other unforeseen challenges will need to be overcome, producing polarisedbeam measurements from these observations presents a unique opportunity for thecalibration of CHIME that would be worthwhile to pursue in the future.Finally, I reported on an unexplained feature in the holography measurements.Although I have not found an explanation for the phenomenon, I have documented itsdefining features and ruled out some potential effects that have similar characteristics.Future investigation is warranted, which may be guided by the analysis presented inthis chapter.75Chapter 4Cross-correlation of CHIME andLyman-α forest MeasurementsThe challenge of separating the 21 cm signal from the extremely bright foregrounds,and the associated calibration requirements, can be avoided to some extent by com-bining CHIME data with observations from other surveys in cross-correlation. Datafrom a different instrument may be affected by different systematics, and in particularif it is observing another tracer of large-scale structure, will be affected by differentforegrounds. Where the fields of the two surveys overlap, only the cosmological sig-nal will be common between them, whereas the foregrounds and other sources ofsystematics will be uncorrelated.For CHIME, such an analysis has already been performed, by cross-correlatingagainst the Extended Baryon Oscillation Spectroscopic Survey (eBOSS) catalogueof distant quasars and galaxies (CHIME Collaboration et al., 2022b). The analysiswas performed by stacking CHIME observations on the known locations of eBOSSobjects, whose environments are also rich in HI. The resulting boost in signal-to-noise was sufficient to overcome the residual foreground contamination present in theCHIME data, as currently filtered, and detect 21 cm signal with > 10σ significancein the lower part of the redshift range accessible to CHIME (0.8 < z < 1.2). Thisconstituted the first detection of 21 cm signal from CHIME and even more significantlythe first detection of cosmological 21 cm signal from any inteferometer.Since that result was published, I have carried out a cross-correlation of CHIMEdata with another tracer of large-scale structure: observations of the Lyman-α forestin quasar spectra from eBOSS. These overlap the CHIME data at higher redshifts,enabling a detection at an average redshift of z¯ = 2.3 (1.8 < z < 2.5), the mostdistant 21 cm measurement reported in the literature to date. In this section I will76describe the work that led to this result.4.1 Motivation/BackgroundThe Lyman-α forest refers to a dense collection of absorption lines measured in thespectra of distant quasars. These are generated by clouds of neutral hydrogen alongthe line of sight to the quasar that have a large opacity to Lyman-α photons. Theyabsorb the continuum emission of the quasar at the Lyman-α wavelength, whichis shifted with respect to the quasar’s restframe as the radiation travels through theexpanding Universe. The measured wavelength of the lines thus provide the redshift ofthe absorbing HI, and their amplitude tracks the optical depth, which is proportionalto the density of HI (Croft et al., 1998).The restframe wavelength of Lyman-α emission is 121.567 nm, which can be ob-served using optical telescopes from the ground at redshifts & 2, corresponding to red-shifted wavelengths & 360 nm (where atmospheric UV absorption becomes limiting;du Mas des Bourboux et al. (2020)). Since it probes a skewer along the line-of-sightto a quasar, the Lyman-α forest measures HI in the inter-galactic medium (IGM).The largest Lyman-α forest survey to date was published by eBOSS in their finaldata release DR16 (du Mas des Bourboux et al., 2020), and has significant overlapwith the redshift range and sky coverage available to CHIME, allowing for a cross-correlation analysis to be carried out. Lyman-α forest measurements are characterizedby absorption of radiation whereas the 21 cm signal is in emission, so we would expectthe cross-correlation with CHIME maps to be negative at small separations, which isa unique signal to search for. This is a feature that has been highlighted by Carucciet al. (2017), who also emphasize the usefulness of such a cross-correlation to mitigateforeground contamination and break degeneracies that arise between the cosmologicalobservables and biases specific to these tracers.4.2 Data4.2.1 eBOSSThe Lyman-α forest measurements used in this analysis are those from the eBOSSsixteenth data release (SDSS DR16; du Mas des Bourboux et al. (2020)). The dataproducts consist of fractional flux transmission along the line of sight to backlight772.00 2.25 2.50 2.75 3.00 3.25 3.50 3.75 4.00z0123count1e6eBOSS pixelsCHIMEFigure 4.1: Histogram of redshifts of eBOSS Lyman-α forest samples and theiroverlap with the CHIME data.quasars, defined asδq(λ) =fq(λ)F¯ (λ)Cq(λ)− 1, (4.1)where fq is the measured flux, Cq is the unabsorbed quasar continuum and F¯ is themean transmission. The product F¯ (λ)Cq(λ) was modelled and fit to the data in orderto produce the catalogue of δq spectra as explained in du Mas des Bourboux et al.(2020). Apart from interpolating onto the CHIME frequency band (Section 4.3.1),no additional processing is performed on the catalogue provided by eBOSS. A fewrandomly selected spectra are shown in the right column of Figure 4.10.There are 210 005 background quasars in the eBOSS Lyman-α catalogue, at red-shifts 2.1 < z < 4, with spectra probing structures down to z > 1.96. The CHIMEdata go down to 400 MHz, corresponding to 21 cm radiation redshifted at z = 2.5.The resulting overlap includes a total of ∼ 3 × 107 spectral samples in the eBOSSmeasurements, with a redshift distribution shown in Figure 4.1. We include the fulleBOSS spectra in this total and the analysis, including the regions blueward of theLyman-β wavelength (∼ 102 nm) at the quasar redshift (referred to as the Lyβ regionin the eBOSS papers). This part of the spectrum is treated separately in their analy-sis because it can include absorption from both transition levels, introducing a shiftedcopy of the measured large-scale structure, but since we are cross-correlating with anexternal dataset, there is no concern of contamination from Lyβ correlations. Thequasars are distributed throughout the SDSS North and South Galactic Cap fields(NGC and SGC), which are fully contained within the sky area observed by CHIME,7890 0180Horizon at CHIME0 570number of quasarsFigure 4.2: Number of background quasars in the eBOSS Lyman-α sample,binned according to their location on the sky in celestial coordinates.The horizontal dashed line shows the southern horizon local to CHIME.The Galactic plane is drawn as a solid line.as illustrated in Figure 4.2.4.2.2 CHIMEThe CHIME data that wer used in this work are composed of nighttime observationsfrom 88 days spanning the calendar year 2019. This amounts to ∼ 1000 hours oftotal observation time, distributed approximately uniformly over the sidereal day.The full frequency range observed with CHIME is 400–800MHz, targeting 21 cmemission from redshifts 0.8 < z < 2.5, but in this analysis I use only the lowestquarter of the band, i.e. 400–500MHz MHz or z > 1.8, which has overlap with theLyman-α forest measurements. This is a similar set of observations to what was usedin CHIME Collaboration et al. (2022b) to detect 21 cm emission in cross-correlationwith SDSS galaxy catalogues. However, that analysis considered data from the upperpart of the frequency band, 585–800MHz, whereas this work is the first report on thehigh-redshift end at 400–500MHz.The lower part of the band is generally more affected by intermittent RFI, andthe overall foreground power is greater due to the negative spectral index of the syn-chrotron radiation that is the dominant component of the sky at these frequencies.790306090120150180Time [sec]145140135130125120Power [dBm]400 450 500 550 600 650 700 750 800Frequency [MHz]140130120110Power [dBm]Figure 4.3: RFI measured by CHIME antennas at the DRAO on May 4 2019.(Top) A spectrum against time showing power in the X polarisationmeasured using a pair of CHIME feeds on adjacent cylinders. The fre-quency resolution is 195 kHz, obtained by direct transform of the rawdata rather than through the CHIME acquisition pipeline. The timeresolution is 100ms. (Bottom) The median value over the time rangeshown in the top panel. This Figure was reproduced from CHIME Col-laboration et al. (2022a).This is apparent in the spectrum displayed in Figure 4.3. Another difference betweenthe lower and upper bands is the presence of spatial aliasing in the data (see Ap-pendix C). Although this simplifies the analysis in principle in the lower band, alias-ing in the upper band was effectively removed by deconvolving the beam response ofthe instrument during the map-making process described in CHIME Collaborationet al. (2022b), so it should not be expected to be a significant difference between theseanalyses.The acquisition process, including real-time processing, was described earlier inSection 1.3.2. The subsequent offline processing pipeline that was developed to turn80these observations into measurements of the 21 cm emission along the lines-of-sightto the background eBOSS quasars will be discussed in Section 4.5 of this Chapter.4.3 Cross-correlation MethodThe cross-correlation function of the 21 cm brightness temperature T21 (i.e. the mea-sured surface brightness from cosmological 21 cm emission) and the Lyman-α forestrelative flux transmission δα is defined at a given redshift asξ(r, r′) = 〈T21 (r) δα (r′)〉 , (4.2)where r and r′ are three-dimensional comoving positions and the angle brackets indi-cate an ensemble average. Assuming statistical isotropy, the correlation function onlydepends on the magnitude of the separation in comoving space, ξ(r, r′) = ξ(|r− r′|),and is related to the power spectrum P (k) by a Fourier transform. The power spec-trum is the summary statistic that is most often used to constrain cosmology, as itcontains a wealth of information about the large-scale structure and its evolution,and it is a convenient tool for modelling.However, we measure these fields at angular positions on the sky and as a functionof wavelength, which maps directly to redshift z, not comoving distance r. Every red-shift also corresponds to a different instant in cosmic time. In other words, we observethe Universe on our past lightcone – in redshift space and on the sphere – denotedby (z, nˆ). Transforming from these observer coordinates to comoving coordinates re-quires knowledge of how the expansion of space proceeded – the Hubble rate H(z).Still, any contributions to the observed redshift that are not due to the expansion, likepeculiar velocities of the measured structures, will alter the inferred distance, an effectknown as redshift-space distortion (Kaiser, 1987). On large scales, the velocities thatare responsible for redshift-space distortions are related to the density fluctuationsthemselves, and so they can actually be used to improve constraints on the cosmology(e.g. Percival & White (2009)). In any case, these effects can straightforwardly beincluded in models of the observations, and for the purpose of this work – which isto detect a signal – it will be sufficient to evaluate the correlation function in redshiftspace directly: ξ (∆z,∆θ).In aggregate, the line-of-sight measurements from eBOSS and CHIME containthree-dimensional information on the two fields, as stated in the dependence of thecorrelation function on separations in redshift (∆z) and angular (∆θ) coordinates.However, CHIME’s resolving power is much greater in frequency than it is on the sky,81so we expect the signal-to-noise to be stronger along this axis. Thus there is not muchto gain from the added complexity of a 3-d analysis. In CHIME Collaboration et al.(2022b) we tried evaluating a 3-d correlation and found no significant improvement tothe detection by including angular separations. In any case, including 3-d informationwould add significant complexity to the analysis, so for a first detection attempt Ichose to only consider line-of-sight correlations.In this section, I will describe the procedure used to compute an estimate of theradial cross-correlation function in redshift space, corresponding toξ (∆z) ∼〈∫dz T i21 (z) δiα (z +∆z)〉i, (4.3)where the angle brackets denote the average over individual line-of-sight observa-tions labelled by i. I describe how redshifts are determined for each observable inSection 4.3.1, and how the cross-correlation is computed in Section 4.3.2.4.3.1 Frequency RemappingThe eBOSS Lyman-α forest spectra measurements are arranged in wavelength binsuniformly spaced in lnλ and spanning ∼360–700 nm, whereas the CHIME data areuniformly spaced in frequency in the 400–800MHz band. In order to correlate thetwo, it will be useful to map these onto a common grid in redshift. Given that theresolution of the CHIME data is coarser by a factor of ∼ 1.2−1.5 in the redshift rangewhere they overlap, it will be convenient to work with the CHIME data directly, Ichoose to map the Lyman-α data onto the CHIME band. The mapping is defined bymatching the frequency of Lyman-α absorption to the observed frequency of 21 cmemission for the corresponding redshift, a relationship given by the ratio of their restwavelengthsν ′ =λαλ21ν. (4.4)I use the reverse Lanczos interpolation method described in detail in Appendix Ato regrid the Lyman-α forest spectra onto the CHIME band. The regions of the spec-tra that are empty (redshifts larger than that of the background quasar, or otherwisemasked) are set to zero in the regridded data.824.3.2 Correlation Function EstimationI will write the discrete samples in frequency as T21[n] = T21(n∆ν + 400 MHz), withintegers n labelling samples at the CHIME frequency resolution, ∆ν = 390.625 kHz.Both the Lyman-α and CHIME spectra have noise estimates associated with everysample so I use these to weight the observations and improve signal-to-noise. Withinverse variance weights, w = σ−2, the cross-correlation function is estimated asξ[n] = W−1M∑iN∑jwi21[j]wiα[j + n]Ti21[j]δiα[j + n], (4.5)where i labels lines of sight, n indexes separation in frequency, and W−1 is a normal-isation term given byW =M∑iN∑jwi21[j]wiα[j + n]. (4.6)M = 210005 is the total number of lines of sight and N ≤ 255 is the number ofsamples along the frequency axis being averaged together at that separation.4.4 Foreground Delay FilterThe defining challenge of 21 cm cosmology is the presence of extremely bright fore-ground emission that buries the signal, by a factor ∼ 105 in the CHIME band,(Section 1.1.4). Separating the signal from foregrounds relies on their very differ-ent spectral properties, but is complicated by instrumental effects imprinted on thedata by the measurement process. Many methods have been proposed, and in recentyears have been put to the test as data from 21 cm experiments like CHIME hasbecome available. One common method is to filter the data in delay-space, i.e. theFourier-conjugate of frequency, where the slowly-varying-in-frequency foregrounds areconfined to low delays. This technique has the advantage of necessitating relativelyminimal knowledge of the instrumental response, at the expense of considerable signalloss. The delay-filter used in the analysis will be described in this section, followinga discussion of the delay representation of the data and how the foregrounds aredistributed in that space.834.4.1 Delay Power SpectraThe delay transform of a signal is simply its Fourier transform along the frequencyaxis, i.e.V˜ (τ) =N∑nV (νn)e−2piiνnτ , (4.7)for a signal V sampled discretely at N frequencies νn. It is often useful to considerthe average power at each delay, the power spectrumP (τ) =〈∣∣∣V˜ (τ)∣∣∣2〉 , (4.8)where the angle brackets denote the average along some axis, e.g. time. This quantityis particularly relevant in the context of interferometric observations of the 21 cmline for cosmology, because it is closely related to the power spectrum of large-scalestructure fluctuations P (k), as we shall see.Visibilities in the Flat-sky ApproximationRecall the visibility equationVij(ν) =∫d2nˆ |A(ν, nˆ)|2 S(ν, nˆ)e−2piibij ·nˆν/c , (4.9)where the integral is over the sky brightness S as seen through the telescope beam Afor a baseline bij. The baseline in units of wavelength uij = bij/λ defines a particularmode of the sky being measured. The more baselines that are recorded (i.e. the morecoverage in u-space), the more faithfully an image of the sky S(nˆ) can be recovered.If we restrict the integral to a small area of sky such that curvature can be ne-glected, which is appropriate if the extent of the beam A is small (as is the case forlarge aperture telescopes), then these modes can be simply interpreted as rectangularFourier modesVij(ν) ∝∫dθxdθyS(ν,θ)e−2piiuij ·θ, (4.10)where θ is a Cartesian vector spanning the small patch of sky being integrated over.This is known as the flat-sky approximation. These small angular separations canfurther be related to physical sizes if one knows the comoving distance to the observedfield, χ, which means that the visibility is sensitive to a specific scale of large-scale84structure, for which the baseline vector is a proxy:uij =χ2pik⊥ij . (4.11)Thus the visibilities are directly measuring cosmological k-modes, specifically of struc-tures on a plane transverse to the line of sight. The transverse component of thewavenumber is labelled k⊥.In addition to structures in the transverse plane, the visibilities measure 21 cmsignal along the line of sight via the relationship between observed frequency andredshift 1 + z = ν21/ν, which is in turn related to the comoving distance along theline of sight (Section 1.1.1). Note that peculiar velocities of the radiating objects alsocontribute to the measured redshift, leading to redshift space distortions, as discussedearlier. So there exists a mapping between observed frequency and distance along theline of sight, such that a visibility in the flat sky approximation constitutes a 3Dmeasurement of the large-scale structure where two of the directions are sampled inFourier space. We can complete the Fourier-representation by taking the transformalong the line of sight, i.e. frequency, as a function of delay as in Equation (4.9) above.The delay τ is conjugate to frequency ν, which is a proxy for distance along the lineof sight, so it corresponds to longitudinal k-modes: k‖. The details of the mappingbetween these quantities depends on the cosmological model and is not uniform acrossa wide band of frequencies, but within a narrow region of frequency (redshift) centredon ν0 it can be approximated ask‖ ≈ 2piH(z)cν20ν21τ , (4.12)where H(z) is the Hubble rate at the corresponding redshift.In the flat sky approximation, the delay-transformed visibility can therefore beinterpreted as a measurement of the Fourier modes of the large-scale structure, sub-ject to a coordinate transformation from the space where observations are performed,baseline and delay (uij, τ), to comoving wavenumbers along the transverse and lon-gitudinal directions (k⊥, k‖). Similarly, the delay power spectrum of visibilities canbe thought of as a measurement of the cosmological power spectrum P (k⊥, k‖). Thisis why in the field of 21 cm cosmology, delay power spectra are a common way tovisualise and analyse the observations, typically taking the form of a 2D image whereone axis spans baseline length (as a proxy for∣∣k⊥∣∣) and the other delay (as a proxyfor k‖). Figure 4.4 shows a delay spectrum evaluated from CHIME data.85Of course, given the large field of view of CHIME, the flat sky approximation isnot appropriate in most contexts and the visibilities are better understood as relatedto modes in a spherical harmonics decomposition. But conceptually they can beinterpreted in a similar manner, as measurements of Fourier modes of the sky. Onthe other hand, the delay transform along the line of sight can be applied directly tothe CHIME data and interpreted as discussed.4.4.2 Foreground WedgeThe characteristic feature of astrophysical foregrounds at the frequencies relevant to21 cm cosmology is their spectral smoothness, typically varying in frequency accordingto a power law to a high degree of accuracy (Santos et al., 2005). In delay space, thisproperty corresponds to being confined to small delays, and can be used to separatethem from signal present at high delays. The region where foreground contaminationis concentrated depends linearly on baseline length, known as the foreground wedge(Morales et al., 2012).For illustration, consider a foreground source with a completely flat spectrum andsome arbitrary distribution on the sky S(ν, nˆ) = S(nˆ), and ignore the beam response,such that the visiblity equation (4.9) is simplyVij(ν) =∫d2nˆS(nˆ)e−2piiuij ·nˆ , (4.13)and its delay transformV˜ij(τ) =∫d2nˆS(nˆ)∫dνe2pii(τν−uij ·nˆ) . (4.14)Because the foreground power and instrumental response are frequency independentin this example, the only contribution to power at non-zero delays comes from thegeometric delay of the baseline-sky projection uij ·nˆ. For a given baseline, this quantityvaries with nˆ and is maximal when they are in the same plane. For a baseline tangentto the surface of the Earth, this occurs when nˆ is on the horizon. The correspondingmaximum delay isτ ijhzn =|uij|ν=|bij|c, (4.15)scaling linearly with the length of the baseline, forming the “wedge” that can be seenin Figure 4.4. In that figure, however, there is significant contamination outside ofthe wedge, and this is because for a real instrument like CHIME, the instrumen-86Figure 4.4: Delay power spectrum of CHIME visibilities evaluated over the RAaxis of a sidereal stack of 88 days of observations. From left to right, eachpanel shows the visibilities corresponding to 0–3 cylinder separations,with the NS baseline length varying along the horizontal axis. The dottedlines bound the foreground wedge as the maximum delay expected for flatspectrum foregrounds at the horizon (Section 4.4.2). The X polarisationis shown here, displayed in uncalibrated correlator units.tal response imprints spectral features on the foreground power, which complicatesthe simple wedge picture sketched out here. Correcting for this imprint via precisecalibration can in principle rein the foreground power back to within the wedge.4.4.3 DAYENU FilterThe confinement of foreground power to the wedge makes it possible to excise the bulkof this contamination while preserving the signal at higher delays. This is achievedby constructing a high-pass delay filter and applying it to the visibilities (or to mapsgenerated from the visibilites, as in Section 4.5.4). Designing an adequate filter iscomplicated by the fact that – as will be discussed in Section 4.5.5 – 65% of theband has been masked in order to remove RFI-like features and other narrowband,instrumental artifacts. The DAYENU technique (Ewall-Wice et al., 2021) is used toconstruct a linear filter for the irregularly-sampled spectra that achieves the required87suppression at large spectral scales. In what follows, a brief summary of this techniquethat was adapted from the one published in CHIME Collaboration et al. (2022b) isgiven.The following simple model is assumed for the covariance of the visibilities as afunction of delay τ :C˜(τ, τ ′) = −12τcutδD(τ, τ ′), |τ | ≤ τcut,∆ν δD(τ, τ ′), |τ | > τcut,, (4.16)where the region below τcut is the region of delay space contaminated by bright fore-grounds,  is a small number that corresponds to the ratio of the radiometric noise toforeground variance, ∆ν = 0.390 625MHz is the width of the frequency channel, andδD(τ, τ ′) is the Dirac delta function. This model results in the following analyticalformula for the covariance between frequency channel νm and νn:Cmn = −1 sinc [2piτcut(νm − νn)] + δmn , (4.17)where sinc(x) ≡ sin (x)/x and δmn is the Kronecker delta.To construct the filter, the delay cut τcut and stop-band rejection  are specified.Equation (4.17) is then evaluated for each pair of frequency channels in the frequencyband. Rows and columns of the covariance matrix that correspond to masked frequen-cies are zeroed and the Moore-Penrose pseudo-inverse (Penrose, 1955) is calculatedR =[mTCm]+, (4.18)where m is a vector that is 1 for valid frequencies and 0 for masked frequencies. Thefilter can then be applied to a given visibility,Vhpf(νm) =∑nRmn V (νn) . (4.19)The weights are also propagated through the filtering operation according towhpf(νm) =[∑n(Rmn)2 (σ(νn))2]−1, (4.20)where σ(ν) is the noise variance associated with the visibility V .To help understand how the filter operates, first consider the case without the88frequency mask, and where the covariance C is assumed to be the real data covarianceand is non-singular: R = C−1. C can be decomposed into a set of eigenvectors {ux}with eigenvalues {λx}, some of which will be dominated by foregrounds, {uf}, andothers by signal, {us}. Applying R to the data v is equivalent to projecting it ontothis basis weighted by the inverse of the eigenvalues:Rv =∑sλ−1s us(us · v) +∑fλ−1f uf (uf · v) . (4.21)The key to this method is the realisation that it does not matter what the precisevalues of the eigenvalues are: as long as λs λf with a sufficient dynamic range, andthat the eigenvectors {uf} span the foregrounds, the filter will effectively suppress theforeground-dominated modes. Thus, the assumed covariance C only has to have theseproperties, not represent the data exactly. The DAYENU filter consists in choosing abasis for C that is well-suited to describing foregrounds that are confined to a top-hat region in delay, Equation (4.17) (these turn out to be discrete prolate spheroidalsequence (DPSS)s, as explained in Ewall-Wice et al. (2021)).The application of this filter to the data for the cross-correlation analysis will bedescribed in Section 4.5.5.4.5 Data Processing PipelineTo produce measurements along the line-of-sight to the eBOSS Lyman-α forest quasarsfrom CHIME data, multiple stages of processing were involved. These processingpipelines are summarised in Figure 4.5 and will be described in detail in this section.Software for CHIME data analysis and processing was developed by membersof the collaboration within a modular pipeline framework designed for CHIME, butwhich is also publicly available. The caput (Shaw et al., 2020c) software packagedefines the core components of this framework, which is essentially composed of “con-tainers” that represent structured data and “tasks” that operate on or generate data,as well as a syntax for defining pipelines composed of a sequence of tasks in a “config”file. The pipeline illustrated in the diagram of Figure 4.5 is a visual representationof such a config. Tasks with increasing specificity to the CHIME analysis, along withsupporting analysis code, are organised in the cora, driftscan, draco (Shaw et al.,2020b,a,d), and ch_pipeline (not publicly released) code repositories. These arewritten in Python for the most part, with lower-level code in C where necessaryfor performance, and support multiprocessing as well as distributed computing with89MPI1.These packages amount to a large codebase, to which I have contributed ex-tensively, both with general infrastructure-level code and particular analysis tasks,including all those described below that are specific to the Lyman-α analysis. Detailsof the CHIME analysis pipeline have been published in the CHIME overview and firstcross-correlation detection papers (CHIME Collaboration et al., 2022a,b), and I willonly summarise those steps here to focus instead on the tasks specific to this analysis.4.5.1 Daily ProcessingThe data from the CHIME telescope are saved to disk continuously in time, and thefirst step of the analysis pipeline is to collect them into individual sidereal days andregrid them onto a common axis in RA to produce a set of repeated observations ofthe sky. At this stage, RFI that was not caught by the real-time flagging is identifiedand masked, along with daytime and any other region of time and frequency that wasflagged for contamination. The list of tasks is as follows.• Apply a timing correction. The ADC boards that digitize the CHIME signalssample at a rate derived from a common 10MHz clock, which is distributedto them via a network of coaxial cables. The F-engine is divided physicallybetween two shielded containers near the telescope, and again within them intoeight “crates” where the ADC boards are located. Each of these experiences aslightly different thermal environment, which leads to differential drift betweenthe clock signals. Corrections are derived for this effect from a noise sourcesignal that is distributed to each crate. A thorough description of this systemis given in CHIME Collaboration et al. (2022a).• Correct telescope rotation. The CHIME cylinders are not oriented perfectlyalong geodetic North, offset by about 0.07◦. This introduces a small shift tothe nominal positions of the antennas, and thus to the gain calibration, whichis accounted for in this task.• Mask decorrelated cylinders. On rare occassions, communication errors withinthe F-engine can cause certain data streams in the correlator to become mis-aligned with all other data streams. This results in one quarter of the feedsdecorrelating with all other feeds for 64 of the 1024 frequency channels, whichpersists until the data streams are re-aligned by restarting the correlator. The1https://en.wikipedia.org/wiki/Message_Passing_Interface90Real-time RFI excisionN2 visibilitiesCalibrate gainsAvg. redundant baselinesReal-time processingAveraging sidereal daysAverage within seasonsApply timing correctionAverage over cylinder pairsFlag bad daysAverage over seasonsSidereal stackDaily processingCompute noise weightsConstruct RFI maskSmooth noise weightsApply thermal calibrationRegrid in Earth rotation angleApply flagsLegend Used to define operationDirect operationFigure 4.5: A summary of the pipeline that processes CHIME data from N2correlation products streaming from the telescope to an average sky mea-surement. The box labelled Real-time processing highlights relevant op-erations performed on the data as they are acquired (Section 2). Dailyprocessing contains the stages of the analysis (offline) pipeline that groupand regrid the data onto a common grid in sidereal time (Section 4.5.1),and Averaging sidereal days those that combine these daily measure-ments into a sidereal stack (Section 4.5.2). This figure was reproducedfrom CHIME Collaboration et al. (2022b).9164 decorrelated frequency channels are uniformly spaced across the band andcan cause significant foreground leakage. We perform an automated search forthese decorrelation events and exclude any frequency channel and time affectedfrom further analysis.• Collate baselines. As part of the real-time processing, the N2 correlation prod-ucts are averaged across nominally redundant baselines to reduce the storagerequirements of the archive, but only within cylinder pairs. This stage completesthe averaging across cylinders, such that only unique baselines are passed on.• Group by sidereal day. Boundaries between sidereal days are identified and thedata are grouped accordingly.• Identify and mask RFI. A three-stage scheme to detect RFI is applied to thedata and produces a mask in frequency and time.• Apply thermal gain correction. Gains for every CHIME feed are derived onceper day as part of the real-time pipeline using a bright astrophysical sourceas a reference; however, these drift between calibration events, mostly as aresult of temperature changes throughout the day. A model for the thermalsusceptibility of the analogu signal chain (CHIME Collaboration et al., 2022a)is used to correct this drift according to the recorded air temperature on site.• Smooth noise weights. The radiometric noise is not expected to change appre-ciably on short timescales. In order to reduce the uncertainty on our estimateof the variance of the radiometric noise and also make our estimate less sen-sitive to transient RFI events, a rolling median filter with a 5min window isapplied to the time axis of the inverse variance weight dataset. Any time thatwas masked is ignored when calculating the median and also remains maskedafter the filtering is applied.• Regrid onto sidereal time. The data are interpolated onto a common siderealaxis, spanning 0–360◦ over 4096 samples, using the reverse Lanczos scheme,described in Appendix A.• Mask daytime and flagged periods. The data are masked during times when theSun was above the horizon, the Moon was within 5◦ of meridian, or a period ofpoor data quality was flagged (for a variety of reasons, e.g. when rain or snowwas affecting the performance of the instrument, or maintenance tasks werebeing performed).924.5.2 Sidereal AverageAt this stage, the daily observations are combined to produce a single measurement ofthe sky with full coverage of the sidereal day. Before a day is added to the stack, it isvisually inspected to flag any residual corruption that was not identified algorithmi-cally. This daily validation task was distributed between members of the collaborationwho volunteered time to inspect the output of the daily processing pipeline. A stan-dard set of diagnostic plots are generated for every day, including maps at a referencefrequency, an estimate of the sensitivity in time and frequency, and power spectraof the delay-transformed (i.e. Fourier transform of the frequency axis) visibilities.(CHIME Collaboration et al., 2022b). At least two “good” assessments were requiredfor a day to be included in the stack. In the end, 88 days out of the 166 processedfrom the year 2019 passed these tests and were included in the stack.Correlated noise between antennas appears in the visibilities as a bias added tothe sky signal. This “noise crosstalk” arises from common sources of thermal emissionthat are picked up by multiple antennas. Likely sources include noise from the am-plifiers that is rebroadcast out of the antennas and ground radiation. The crosstalkis observed to be stable over a period of days, so it can be suppressed by removing aconstant offset from each individual day of observations. The bias identified in thisway, however, is completely degenerate with the average sky contribution over thatperiod, so care must be taken to constrain this procedure so that it is consistent acrossall days. The simplest way to do this is to choose a reference period common to alldays within which to estimate the crosstalk, which ensures that any sky signal that issubtracted is identical for every day. However, this is not possible because we consideronly nighttime data, which over the course of a year covers disjoint sections of thesidereal day. The solution that was devised is to proceed with the averaging in twosteps. First, the year is divided into quarters and cross-talk is estimated separatelyin each within a region that is shared by every day in the quarter. Each quarter canthen be stacked consistently. The second step is to combine the quarter-stacks into afull sidereal average, requiring consistency in the regions where they overlap to finda solution. A consequence of this crosstalk subtraction procedure is that the averagesky signal over the sidereal day is lost, due to it being degenerate with the crosstalk.The implementation of this procedure is discussed in detail in CHIME Collaborationet al. (2022b).To allow for jack-knife tests, the available days are also divided into two parti-tions, with alternating days in time assigned to each (“even” and “odd” days). Thesetwo partitions are separately averaged in the same way as described above to pro-93duce independent but nominally redundant measurements of the sky. The differencebetween these will be used to test for systematic effects later.4.5.3 Point-Source SubtractionThe CHIME feeds have non-negligible off-axis response. As a result, the signal fromthe four brightest point sources – Cygnus A, Cassiopeia A, Taurus A, and VirgoA – is significant whenever these sources are above the horizon. The first stage ofprocessing following the sidereal average is a targeted removal of these four sources.The description of this task was published in section 4.1 of CHIME Collaborationet al. (2022b).The following model is assumed for the contribution of the four brightest pointsources to the visibility measured by baseline b at local Earth rotation angle φ andfrequency ν:Vpsrc(b, ν, φ) =4∑s=1as(ν, φ) ej2piνb·nˆ(θs,φ−φs)/c , (4.22)where as(ν, φ), θs, and φs denote the primary-beam-modulated amplitude, declina-tion, and right ascension of source s, nˆ is the unit vector pointing towards the source’slocation, and c is the speed of light. At every frequency and local Earth rotation anglewe estimate the set of source amplitudes a using a weighted linear regression. Themodel for the four brightest sources is then computed using (4.22) and subtractedfrom the data. Note that only visibilities measured by baselines consisting of feeds ondifferent cylinders are used to solve for the source amplitudes – because contaminationfrom diffuse Galactic emission and noise crosstalk is significantly reduced for theseinter-cylinder baselines – but the resulting model is subtracted from all baselines.4.5.4 Map-MakingIn order to evaluate the cross-correlation function along the lines-of-sight to theeBOSS backlight quasars as laid out in Section 4.3.2, the filtered CHIME visibili-ties must first be made into maps of the sky, from which pixels at the locations of thequasars can be extracted. The map-making method used in this analysis is discussedin great detail in Section 4.3 of CHIME Collaboration et al. (2022b), from which thedescription below has been adapted and summarised. This method is very similar tothe ringmaps described in Section 2.3, but with some additional steps.The first step is the same, and involves taking a linear combination of the NS94baselines to form beams along the meridian. Since the CHIME baselines are arrangedin a regular grid of NS and EW separations, let us introduce the notation Vxy(ν, φ) torepresent a CHIME visibility for the (x, y) baseline, where x indicates EW separationin grid units and y along NS, at frequency ν and RA φ. The unit grid separationsare dx ≈ 22m and dy ≈ 0.3m. The three cylinder separations and 256 antennasper cylinder mean that the baseline grid spans ±3dx and ±255dy. The beamformedvisibilities are thenVx(ν, φ, θ) =∑yWxy(ν, φ)Vxy(ν, φ)e−2piiydyν sin(θ−θL)/c (4.23)for a declination θ, such that θ− θL is the zenith angle on meridian at the telescope’slatitude θL. The relative weightsWxy(ν, φ) = wxy(ν, φ)[∑y′wxy′(ν, φ)]−1(4.24)are computed from the inverse variance of the visibilities, wxy = [Var (Vxy)]−1, andnormalised such that the flux of a point source will be preserved in the map.In the next step, we attempt to deconvolve the effect of the telescope beam fromthe map, and to understand this it will be useful to write down a model for how thesky is observed by CHIME. The measured visibilities can be decomposed into skyand noise contributionsVxy(ν, φ) = Vxy(ν, φ) + nxy(ν, φ) , (4.25)where Vxy, the sky visibility, is an integral over the sky brightness S(ν, θ, φ) as seenthrough the telescope beam A(ν, θ, φ):Vxy(ν, φ) =∫|A(ν, θ′, φ− φ′)|2 e2piib·nˆν/cS(ν, θ′, φ′) cos θ′dθ′dφ′ . (4.26)Here c is the speed of light and nˆ(θ′, φ − φ′) is the unit vector pointing towardsdeclination θ′ and hour angle ha ≡ φ− φ′, given bynˆ(θ′, φ− φ′) =− cos θ′ sin (φ− φ′) xˆ+ (cos θL sin θ′ − sin θL cos θ′ cos (φ− φ′)) yˆ+ (sin θL sin θ′ + cos θL cos θ′ cos (φ− φ′)) zˆ .(4.27)95The phase in the exponential corresponds to the geometric path-length due to theprojection of the sky direction onto that baseline. The CHIME beam is compact inha (.2.5◦ FWHM) so it is appropriate to restrict the integral to small ha, and writethis phase asb · nˆ = −xdx cos θ′ (φ− φ′) + ydy sin (θ′ − θL) . (4.28)(Note that this expression neglects the rotation of the cylinders from geodetic norththat adds a small correction to the phase, but this is included in the map-makingcode.)Plugging this model for the visibility into Equation (4.23) for the beamformedvisibilities, we can define functions that describe the beamVx(ν, φ, θ) =∫bθˆsynth(ν, θ, θ′, φ) Bx(ν, θ′, φ− φ′)S(ν, θ′, φ′) cos θ′dθ′dφ′ (4.29)whereBx(ν, θ′, φ− φ′) = |A(ν, θ′, φ− φ′)|2 e−2piixdx cos θ′(φ−φ′)ν/c (4.30)will be referred to as the beam transfer function, which captures the response of anindividual feed, andbθˆsynth(ν, θ, θ′, φ) =∑yWxy(ν, φ) e2piiydy [sin (θ′−θL)−sin (θ−θL)]ν/c (4.31)is the synthesized beam in the θˆ direction, the shape of which is determined bythe geometry of the baseline array and the weights used in the linear combination.Although this quantity depends on sky coordinate φ through the inverse varianceweights, in practice this is a very weak dependence and the synthesised beam can bethought of as stationary.Since the beam transfer function does not change appreciably on scales less thanthe FWHM of the synthesized beam, it can be brought outside of the integral over θ′in (4.29), resulting in:Vx(ν, φ, θ) =∫Bx(ν, θ, φ− φ′) cos θ dφ′∫bθˆsynth(ν, θ, θ′)S(ν, θ′, φ′) dθ′ . (4.32)Given a model for the primary beam A, the beam transfer function is computedusing (4.30) and then deconvolved from the hybrid beamformed visibilities at eachdeclination to recover the flux density of the sky S convolved with bθˆsynth. The con-struction of the primary beam model will be described in Section 4.5.4. The fast96Fourier transform of the hybrid beamformed visibility is taken along the φ axis,V˜xm(ν, θ) =1Nφ∑nVx(ν, φn, θ)e−imφn . (4.33)Here the sum runs over the Nφ = 4096 samples on the grid in local Earth-rotationangle and the m-modes range over [−2047, 2048]. We will henceforth refer to thisoperation as the m-mode transform. The same operation is performed on the beamtransfer function. Figure 4.6 provides an example of the m-mode transform of boththe hybrid beamformed visibility and the beam transfer function.The beam transfer function is then deconvolved from the data in m-space usinga Tikhonov regularization scheme:M˜m(ν, θ) =∑xWx(ν) B˜∗xm(ν, θ) V˜xm(ν, θ)η +∑xWx(ν) |B˜xm(ν, θ)|2, (4.34)whereWx(ν) = wx(ν)/∑xwx(ν) (4.35)denotes the relative weight given to each east-west baseline and η is a regularizationparameter. The different east-west baselines measure a largely disjoint set of m-modes, with each baseline primarily sensitive to the range of ms centred onmcenter,x(ν, θ) = −2piν cos (θ)xdx/c , (4.36)with widthmwidth(ν, θ) = 2piν cos (θ)w/c , (4.37)where w = 20m is the width of the cylinder. However, there is some mild overlapthat is dependent on the aperture illumination and accounted for by the m-modetransform of the primary beam pattern. (4.34) first performs a weighted average ofthe measurements made by the different east-west baselines, and then deconvolves theprimary beam by effectively dividing by the corresponding weighted average of them-mode transform of the beam transfer function. The regularization parameter η isthe assumed inverse signal-to-noise. It defines which m-modes are signal dominated,and hence should be divided by the beam, and conversely which m-modes are noise97Figure 4.6: (Top) m-mode transform of the hybrid beamformed visibilities, for0–3 cylinder separations, at a single frequency. The areas with signifi-cant amplitude indicate the scales that the visibilities are sensitive to.(Bottom) m-mode transform of the beam model at these separations andfrequency. Note that the declination axis has been extended beyond 90◦to reflect the sensitivity of the beam to regions of the sky on either sideof the north celestial pole (NCP). At the NCP itself, the beam patternapproaches a δ-function at m = 0 and this causes artefacts in the m-mode transform, which explains the horizontal line there.dominated, and should not be amplified further by dividing by the beam. A value η =10−4 was chosen as a compromise between the effectiveness of the beam deconvolutionand the resulting noise in the map.We setwx(ν) =0, x = 0,[σx(ν)]−2 , |x| > 0, (4.38)where[σx(ν)]2 =1N2φ∑n[σx(ν, φn)]2 (4.39)is the variance of the noise in the m-mode transform of the beamformed visibility,98V˜xm, and[σx(ν, φ)]2 =(∑ywxy(ν, φ))−1(4.40)is the variance of the noise in the beamformed visibility, Vx. This weighting schememasks all intra-cylinder baselines and propagates the inverse-variance weights throughthe beamforming and m-mode transform for inter-cylinder baselines. The redun-dancy of the array results in Wx(ν) ≈ [0.0, 0.5, 0.333, 0.166] for |x| = [0, 1, 2, 3],corresponding to the intra-cylinder auto-correlation that is removed, the three-foldredundancy in the one-cylinder separation, two-fold redundancy in the two-cylinderseparation, and single appearance of the three-cylinder separation, as a fraction ofthe total number.The intra-cylinder baselines are masked for this analysis because they contain twosources of contamination that are significantly reduced in the inter-cylinder baselines:(1) large-scale diffuse Galactic emission; and (2) noise crosstalk (see discussion in4.5.2). Since the noise crosstalk changes slowly with time, it contaminates only lowm, which is where the signal from the sky resides for intra-cylinder baselines. Notethat the signal from the sky at declinations near the NCP will also appear at low m,even for inter-cylinder baselines.Finally, the deconvolved map is obtained by taking the inverse m-mode transformM(ν, θ, φ) =∑mM˜m(ν, θ)eimφ . (4.41)A map generated in this way is presented in the top panel of Figure 4.9.Beam ModelOur primary beam model is obtained by deconvolving a model for the radio skythat consists only of extragalactic point sources from the visibilities measured withbaselines that have a large east-west component. The long baselines resolve outthe diffuse Galactic emission, making a point-source-only sky model a reasonabledescription of the data. There are several high-resolution, large-area sky surveysthat can be interpolated to the 400−800MHz CHIME band to construct this skymodel. At lower frequencies we rely on the VLA Low-frequency Sky Survey (VLSS;Cohen et al. (2007)) at 74MHz and the Westerbork Northern Sky Survey (WENSS;Rengelink et al. (1997)) at 326MHz. At higher frequencies we rely on the NRAO VLA99Figure 4.7: Slices of the 3D beam model – from left to right – at constant hourangle, declination, and frequency, for the X polarisation. The top rowshows the beam model that was derived from data by fitting to a cata-logue of point sources with known fluxes, which can be compared to asmooth analytic model in the bottom row (discussed in Section 4.6.2).In the last column, the beam pattern has been projected onto sky coor-dinates HA and declination, which is why it broadens significantly nearthe NCP. Note that the declination axis has been extended beyond 90◦to reflect the sensitivity of the beam to regions of the sky on either sideof the NCP.Sky Survey (NVSS; Condon et al. (1998)) at 1400MHz and the Green Bank survey(GB6; Gregory et al. (1996)) at 4850MHz. The method used to deconvolve the skymodel from the data is similar to the method used to construct a map, which wasdescribed in the previous section. Whereas the map maker deconvolves a model forthe primary beam from the hybrid beamformed visibilities to estimate the intensity ofthe sky, the beam calibration deconvolves a model for the sky intensity from the samehybrid beamformed visibilities in order to estimate the primary beam. Appendix Bof CHIME Collaboration et al. (2022b) describes this method in detail. Slices of theresulting beam model are displayed in Figure 4.74.5.5 Foreground FilterThe DAYENU high-pass delay filter that was described in Section 4.4 is applied atthis stage in the pipeline to each pixel independently in the maps. The lower edge100Figure 4.8: Delay power spectrum for a ringmap of a sidereal stack of 88 daysof CHIME observations prior to foreground filtering, evaluated acrossright ascensions between 110 and 270◦, at fixed declination. The colourscale represents power in units of noise power. Vertical dashed lines showthe lower bounds of the DAYENU filter pass-band.of the delay pass-band is set at 200 ns, which was chosen to encompass the bulk offoreground power at all declinations in the delay power spectrum, as illustrated inFigure 4.8. This power spectrum was evaluated using the Gibbs sampling methoddescribed in appendix A of CHIME Collaboration et al. (2022b) and averaging acrossright ascensions 110–270◦, a region corresponding roughly to the eBOSS north galacticcap (NGC) field, where the majority of the Lyman-α quasars are found. Filteringdelays below 200 ns corresponds to removing large scales with k‖ . 0.13 Mpc−1, atthe average redshift z¯ = 2.3 (∼430MHz). Figure 4.9 shows maps at a single frequencyand polarisation before and after the delay filter.Frequency FlaggingAlthough RFI is flagged algorithmically at multiple stages in the real-time and dailyprocessing, low-level events that are fainter than the foregrounds are very hard todetect. This contamination becomes apparent after the high-pass delay filter hasbeen applied to remove most of the foreground power, and is strong enough to obscurethe cosmological signal. I identified affected frequencies by iteratively filtering and101flagging outlier frequencies until none remained. The flagging was performed byvisually inspecting maps generated from the filtered data at every frequency, wherethe spatial information helps identify bad cases. After this final stage of frequencyflagging, 65% of frequencies in the 400–500MHz band are masked out.In addition to RFI, instrumental effects that have not yet been identified or cor-rected for could have a similar effect. Even small calibration errors can be sufficientto cause significant foreground leakage once the bulk of the power in the maps hasbeen filtered, especially when these errors affect isolated frequencies, as those willpropagate widely in delay. As the calibration and flagging parts of the daily pipelinecontinue to improve, it may be possible to recover some of these flagged frequencies.Outlier MaskPixels near foreground features such as bright point sources tend to show an excessin power even after filtering (see Figure 4.9), so we also derive a mask that excisesany pixel above a 6σ threshold determined by the expected thermal noise. 5–10%of unflagged pixels were masked at this stage, reflecting the significant non-Gaussiantails of the pixel value distribution. These will be discussed in Section 4.8. The6σ threshold was chosen during the analysis of the first CHIME cross-correlationdetection (CHIME Collaboration et al., 2022b) based on a series of signal-injectionsimulations that determined the bias on the 21 cm signal amplitude introduced bythis choice to be small. In this analysis, I do not attempt to interpret the signalamplitude for cosmology, and the chosen threshold is sufficient to detect a cross-correlation signal, so I will leave further tuning of this parameter in the lower bandfor future work.Apart from some minor changes, the pipeline stages outlined so far are the sameas were described in detail in CHIME Collaboration et al. (2022b).4.5.6 Extraction of Line-of-sight SpectraThe catalogue of eBOSS Lyman-α forest spectra is provided in structured data filesusing the FITS format2, which contain the fractional flux transmission spectra (Equa-tion 4.1) on the line of sight to every background quasar. The CHIME pipeline oper-ates on data containers that require a specific structure and are packaged in HDF53files, so I implemented pipeline tasks to load the eBOSS FITS files and prepare the2The structure is defined here: https://data.sdss.org/datamodel/files/EBOSS_LYA/DELTA_LYAB/delta.html3https://www.hdfgroup.org/solutions/hdf5/102Figure 4.9: (Top) Ringmap generated from inter-cylinder baselines (&22m EWdistance) at 476.56 MHz for the XX polarisation. (Middle) Ringmapgenerated from visibilities filtered using DAYENU with a 200 ns delay cut.(Bottom) The same as the middle panel but with an outlier mask derivedfrom a 6σ threshold based on the expected noise. In all panels, verticaldashed lines indicate the region that was used to estimate the delaypower spectrum in Figure 4.8.data for this analysis. The following operations are performed.• Load the FITS files; concatenate and interpret the contents into a single pipelinecontainer.• Combine the Lyman-α and Lyman-β regions of the spectra into continuousobservations for every quasar (these are provided as separate catalogues byeBOSS, as discussed in Section 4.2.1).• Interpolate the spectra from their native uniform sampling in lnλ to theirredshift-equivalent on the CHIME frequency grid (as defined in Section 4.3.1,Equation 4.4). This is implemented with the reverse Lanczos interpolationscheme described in Appendix A.1031.9 2.0 2.1 2.2 2.3 2.4 2.50.010.00RA 133.2 DEC 52.0CHIME1.9 2.0 2.1 2.2 2.3 2.4 2.5101eBOSS0.0050.0000.005RA 354.8 DEC -1.3 100.00250.00000.0025RA 358.7 DEC 18.8 101.9 2.0 2.1 2.2 2.3 2.4 2.5z0.0020.0000.002RA 15.6 DEC 8.11.9 2.0 2.1 2.2 2.3 2.4 2.5z10400 420 440 460 480 500 (MHz)3600 3800 4000 4200(Å)Figure 4.10: (right) Randomly selected spectra from the eBOSS catalogue ofthe fractional transmission to background quasars. (left) Correspond-ing spectra extracted from the maps generated from filtered CHIMEdata, in relative units of surface brightness. Some regions of the spec-tra were set to zero according to the mask derived as explained inSection 4.5.5. The dual horizontal axes label the redshift of emission(21 cm) or absorption (Lyman-α) by HI on the bottom, and the ob-served frequency/wavelength on the top.The result is a single data container that holds a catalogue of background quasarpositions on the sky and in redshift, and their associated Lyman-α forest spectraon a redshift grid that corresponds to the frequency of 21 cm emission observed byCHIME. The counterparts to these spectra in the CHIME data are extracted from thefiltered maps by finding for every quasar the CHIME pixel that covers its position andstoring the frequency spectrum of that pixel in the catalogue. Figure 4.10 shows a fewrandomly selected spectra from the eBOSS catalogue and the corresponding filteredCHIME data. The mapping between redshift and the observed frequency/wavelengthis also shown.1044.5.7 Correlation FunctionThe final task in the pipeline implements the cross-correlation method that was de-scribed in Section 4.3.2. This is essentially evaluating a weighted convolution forevery eBOSS and CHIME spectrum, an operation with a computational cost thatscales with O(N2) for spectra of length N (N = 256 in this case). Because this steprequires a significant amount of computation and no appropriate methods exist inthe numpy library that provides fast computational tools for Python, I implementedthe correlation function estimation in a C code that is wrapped into the higher-levelpipeline task.I also explored an alternative method for evaluating the cross-correlation functionin Fourier-space, where the convolution becomes a product. The advantage of thismethod is that it can use the fast fourier transform (FFT) method to reduce thecomputational scaling to O(N lnN). Both methods have been integrated into thepipeline, and they produce similar results. However, the direct convolution is easierto interpret given non-uniform weights, and the computational cost is not prohibitive,so it was chosen over the FFT. This will be discussed in more detail in Section 4.7.2.Nevertheless, the FFT method was useful in developing the analysis, as its speedfacilitated running through many iterations.4.6 Signal Model and SimulationsTo interpret the results of this analysis and assess detection significance, I generatedsimulations of correlated 21 cm and Lyman-α forest measurements. These are runthrough the same processing pipeline and analysis as the data in order to forwardmodel the cross-correlation result. Although the models described below are para-metric, I did not attempt to fit them to the data. Instead, the purpose of these modelsis to generate a reasonably realistic template to compare to the measurement.4.6.1 Large-scale Structure SimulationThe common basis for the simulated 21 cm and Lyman-α fields is a map of the dis-tribution of matter on large scales, the matter overdensity δm(x, z) ≡ [ρm(x, z) −ρ¯m(z)]/ρ¯m(z), where an overbar denotes a spatial average. A detailed descriptionof the model used to generate maps of large-scale structure in the CHIME redshiftrange is given in Section 5 of CHIME Collaboration et al. (2022b), which is brieflysummarised here.As input to our simulations, we use the halo model prediction for the nonlinear105matter power spectrum from Mead et al. (2021), as implemented in the CAMB code(Lewis et al., 2000). This model is evaluated at z = 1 and transformed to a 3dcorrelation function using the hankl Python package (Karamanis & Beutler, 2021) viathe FFTlog method. We additionally employ Richardson extrapolation to repeatedcomputations of the correlation function with increasingly fine k sampling in order toreduce numerical errors. We then transform this to a multi-frequency angular powerspectrum C`(ν, ν ′), and perform further frequency integrals over top-hats with width0.390 625MHz in order to mimic the effect of CHIME’s frequency channelization. Weform a set of Nν HEALPix maps (Górski et al., 2005) from a Gaussian realization ofthis angular (matter) power spectrum and use the linear growth factor for our fiducialcosmology to scale each map to the redshift corresponding to its frequency.4.6.2 CHIME Maps SimulationFrom the underlying large-scale structure simulation, mock CHIME observations arederived by modelling the HI density and its associated 21 cm emission that traces thematter field δm, then simulating the measurement process leading to visibilities. Thedescription given here is adapted from what was published in CHIME Collaborationet al. (2022b).21 cm Signal ModelThe model for the HI overdensity δHI assumes that it is a linearly biased tracer of thetotal matter density δm:δHI(k; z) =[bHI(z) + f(z)µ2]D˜FoGHI (kµ, z)δm(k, z) , (4.42)where bHI is the bias factor (assumed to be scale-independent), and the f(z)µ2 termencodes the effect of redshift-space distortions at linear order (Kaiser, 1987), withf as the logarithmic growth rate and µ ≡ k‖/k. We aim to capture the key non-linear contributions to the two-point statistics of the fields: we include real-spacenon-linear clustering in δm itself, and the impact of small-scale velocities on redshift-space observations (“Fingers of God”; Jackson 1972) is modelled with the dampingfunction D˜FoGHI . Each of these ingredients is described in the following paragraphs.For the linear bias of HI, we follow Cosmic Visions 21 cm Collaboration et al.(2018) in smoothly interpolating between measurements from the IllustrisTNG sim-ulations (Villaescusa-Navarro et al., 2018) at z < 2 and the analytical model from1061.0 1.5 2.0 2.51.52.0b HI1.0 1.5 2.0 2.51.01.52.0P (h1 Mpc)1.0 1.5 2.0 2.5z0.60.8103HI1.0 1.5 2.0 2.5z0.20.3T b (mK)Figure 4.11: Modelled parameter values described in this section as a functionof redshift.Castorina & Villaescusa-Navarro (2017) at z > 2.4 We show our bias model in Fig-ure 4.11.We model the Finger of God damping in Fourier space as a Lorentzian:D˜FoGHI (kµ, z) =11 + k2µ2σP(z)2/2, (4.43)where the damping scale σP can approximately be associated with the pairwise ve-locity dispersion of HI emitters on nonlinear scales. The (constant-redshift) Fourierconjugate of this function is an exponential in comoving distance (e.g. Scoccimarro2004),DFoGHI (x‖, z) =e−∣∣x‖∣∣√2/σP(z)√2σP(z), (4.44)and we implement the Finger of God effect by convolving our simulated maps withthis kernel along the line-of-sight axis. This is equivalent to multiplying the 3d powerspectrum by D˜FoGHI (kµ, z)2. We choose the damping scale based on simulations fromSarkar & Bharadwaj (2019), who attempt to account for the motion of HI withingalaxies in addition to the contribution from the velocity dispersion within darkmatter halos. They assume that the Finger of God damping of the 3d HI powerspectrum is given by a Lorentzian, and fit a σP(z) relation to their simulations. Weuse these results, multiplied by a factor of 2−1/2 to translate to the damping given4This bias model has been implemented in the PUMANoise code, available from https://github.com/slosar/PUMANoise.107by a squared Lorentzian (as implied by Equation 4.43). The adopted σP(z) model iswell fit by a quadratic function of redshift, given byσP(z)h−1Mpc= 1.93− 1.48(z − 1) + 0.81(z − 1)2 . (4.45)Over the redshift range of interest, this σP(z) model is within 20% of the valuesobtained in Villaescusa-Navarro et al. (2018) from fits of a squared Lorentzian tomeasurements from the IllustrisTNG simulations. We plot our model for the dampingscale in Figure 4.11.We convert simulated maps of δHI into brightness temperature fluctuations bymultiplying by the mean 21 cm brightness temperature T¯b(z). After the end of reion-ization, the 21 cm brightness temperature can be written as (Bull et al., 2015)Tb(x, z) =3h¯c3A1016kBν221(1 + z)2H(z)nHI(x, z) , (4.46)where nHI is the comoving HI number density and A10 is the Einstein coefficient forspontaneous emission in the 21 cm line. Using nHI(x, z) = n¯HI(z)[1+ δHI(x, z)], we canwriteδTb(x, z) ≡ Tb(x, z)− T¯b(z) = T¯b(z)δHI(x, z) , (4.47)which justifies our method of converting maps of δHI into Tb. Using n¯HI(z) = ΩHI(z)ρc,0/(mp+me), where ρc,0 = 3H20/8piG is the critical density today5, we can write T¯b asT¯b(z) =[9h¯c3A10H100128pikBν221G(mp +me)]× h H0H(z)ΩHI(z) (1 + z)2 , (4.48)where H100 = 100 km s−1Mpc−1 and h = H0/100. The prefactor in square brackets isindependent of cosmology, consisting only of fundamental constants and A10. UsingA10 = 2.8843× 10−15 s−1 (Gould, 1994), (4.48) can be written more compactly as6T¯b(z) ≈ 191.06[hH0H(z)ΩHI(z) (1 + z)2]mK . (4.49)5Following the common convention in the 21 cm intensity mapping literature, we have definedΩHI(z) in terms of the comoving HI number density at redshift z and the critical density at z = 0.6Other versions of (4.49) in the literature have prefactors that vary significantly from 180–190mK,most of which is accounted for by using values of A10 from older calculations. The value quoted inthe main text is taken from a recent review of atomic transition properties (Wiese & Fuhr, 2009),which takes its A10 value for hydrogen from Gould (1994).108For ΩHI(z), we use the fitting function from Crighton et al. (2015), which wasdetermined from a compilation of ΩHI estimates over 0 < z < 5:ΩHI(z) = 4× 10−4(1 + z)0.6 . (4.50)We plot Eqs. (4.49) and (4.50) in Figure 4.11.21 cm Brightness MapsThe model components just described are combined with the simulated maps of large-scale structure to produce maps of 21 cm brightness. For each map of δm (previoussection), we generate a map of the gravitational potential φ, to which we apply afinite-difference second derivative in the radial direction and appropriate prefactorsto generate a velocity field that is added to the biased matter to include linear redshift-space distortions. These maps are then convolved with the frequency kernel designedto reproduce the desired form of Finger-of-God damping in Fourier space. Finally, themaps corresponding to δHI are multiplied by the mean 21 cm brightness temperatureT¯b(z).Simulated VisibilitiesWe make use of the m-mode formalism (Shaw et al., 2015) to translate simulated21 cm maps into visibilities. In this formalism, the spherical harmonic coefficientsaP`m(ν) of sky maps for Stokes parameter P ∈ {I,Q, U, V } are related to the sidereal-time Fourier transform of the visibility timestream, V˜ pxy,m(ν), via multiplication by abeam transfer matrix Bp,Pxy;`m(ν):V˜ pxy,m(ν) =∑P∑`Bp,Pxy;`m(ν)aP`m(ν) , (4.51)where p ∈ {XX,XY, Y X, Y Y } labels the polarisation pair for that visibility. Afterperforming this multiplication, we convert the result to a visibility timestream byinverse Fourier transforming in m, applying zero-padding to the Fourier transformsuch that the corresponding time resolution matches that of the observed siderealstacks. Gaussian noise at a level consistent with the CHIME observations is optionallyadded to the simulated visibilities.The beam model used to calculate the beam transfer matrix Bp,Pxy;`m(ν) is an ana-lytic description for the beam of a cylindrical telescope, defined in Section 4 of Shawet al. (2015). This model treats the beam as factorisable into an unfocused NS profile109and a diffracted EW part. Although it is not an accurate description of the CHIMEtelescope, the simple model captures the broad features of the beam that determinethe field of view and its scaling in frequency. A comparison of this analytic model tothe data-derived beam that is used to calibrate the observations is shown in Figure 4.7.Up to this point, the simulated data are in brightness temperature units (K). Totransform into spectral flux density units (Jy), which is what the CHIME visibilitiesare calibrated to, we first compute the beam solid angle for the assumed beam model:Ωp(ν) =∫dθ cos θ∫dφ |Ap(ν, θ, φ)|2 . (4.52)We then multiply the visibilities by the standard Rayleigh-Jeans conversion factorand the beam solid angle, normalized by the power beam evaluated at ha = 0′ and areference declination θref:V pxy(ν, φ)∣∣Jy=2× 1026kBν2c2Ωp(ν)|Ap(ν, θref, 0)|2 Vpxy(ν, φ)∣∣K. (4.53)With this normalization, a visibility corresponding to a point source that transitsat θ = θref has an amplitude equal to the flux of the source. For consistency withCHIME’s beam and complex gain calibration, we set θref to the declination of CygnusA, the bright point source that is used as a daily calibrator.4.6.3 Lyman-α Forest SimulationLyman-α Signal ModelIn order to model the Lyman-α forest measurements, we need a prescription forconverting the matter density fluctuation maps δm(z, nˆ) generated by our large-scalestructure (LSS) simulation into an optical depth to Lyman-α photons. A commonlyused and straightforward approach is to model the mildly non-linear baryon field aslog-normally distributed (Bi & Davidsen, 1997) and the Lyman-α absorption by thefluctuating Gunn-Peterson approximation (FGPA) (Farr et al., 2020).The log-normal transform is performed so as to preserve the variance of the originalfield σ2m1 + δLN(z, nˆ) = exp[δm(z, nˆ)− σ2m(z)2]. (4.54)This step also ensures that the density is strictly non-negative.The optical depth τ is proportional to the number density of HI atoms, which the110FGPA assumes is related to the density by a power lawτ(z, nˆ) = τ0(z) (1 + δLN(z, nˆ))α . (4.55)I adopt the parameter valuesτ0(z) = 0.3(1 + z3.4)4.5, α = 1.6, (4.56)quoted in Cieplak & Slosar (2016) and Seljak (2012). These steps are illustrated inthe top three panels of Figure 4.12 for a realisation of LSS along a given line of sight.Rather than use the FGPA, it might be more self-consistent to derive the opticaldepth from the HI density that was modelled for the 21 cm temperature directly(Equation (4.42)). This is expected to be part of a future effort to more carefullymodel the Lyman-α signal and perhaps derive constraints on the model parameters.I will emphasize that for the purpose of this work, the simplest approach was taken.However, the analysis of the cross-correlation result (Section 4.7) will show that asimple FGPA-based model is sufficient to describe the measurement at the currentsignal-to-noise level, so I leave more detailed modelling to future analyses.Mock CatalogueFrom the Lyman-α fractional flux transmission simulation I extract individual spectraalong the lines of sight at a number of sky locations, chosen to be those of the quasarsin the eBOSS catalogue. Each spectrum is also masked so that regions of missing data– redshifts larger than the backlight quasar and instrumental masks – match thosein the corresponding eBOSS data. Gaussian random noise is optionally added tothe simulated spectra at a level consistent with the noise variance recorded alongsidethe data. This results in a mock catalogue with exactly the same quasar sample asthe data, and thus the same sensitivity to cross-correlation with CHIME, but witha synthetic Lyman-α forest signal. The bottom panel of Figure 4.12 shows a mockspectrum generated in this way.It must be noted that these simulations of the Lyman-α flux transmission neglecta number of effects that would be necessary to include in order to produce realisticsimulations of the eBOSS spectra (and which are included in more sophisticatedsimulations, e.g. Farr et al. (2020)). Among the most significant are:• redshift-space distortions.1112.50.0m0.02.5em101.8 2.0 2.2 2.4z101+nFigure 4.12: Fluctuations along a given line of sight in redshift, at subsequentsteps of the procedure for generating simulated Lyman-α forest spectra.From top to bottom: the matter fluctuations from the LSS simulation;the result of the log-normal transform; the Lyman-alpha fractional fluxtransmission; and finally masked and with noise added to emulate thereal spectrum.• metal and Lyman-β line absorption.• high column density (HCD) absorbers.• biases due to the modelling of the quasar spectra continuum.A careful accounting of these will be necessary for future work that will aim toconstrain the physics of the Lyman-α forest.4.6.4 Cross-correlation TemplateThe simulated CHIME visibilities are run through a processing pipeline very similar tothe one that is used on the data, as was described in Section 4.5, with a few differences.This is mainly because the simulations contain signal only (i.e. no foreground sourceswere modelled and included). The difference in processing are listed here.• No bright point-source subtraction is performed, since no point sources areincluded in the simulation.• The beam model that is deconvolved from the map is the analytic one, not thedata-derived model.112Generate mapsApply HI bias Apply linear RSD andFoG damping Convert Generate visibilities Convert K JyApplyfrequency mask Make mapsHigh-passdelay filter Mask outlier pixelsLog-normaltransform FGPAExtract lines-of-sightto eBOSS quasars Apply spectral masks Add gaussian noiseCalculate cross-correlationSimulated CHIME 21 cm measurementsSimulated Lyman- catalogueExtract lines-of-sightto eBOSS quasarsFigure 4.13: Graphical representation of the simulation pipeline, showing theparallel paths generating simulated Lyman-α forest and 21 cm obser-vations. The input to the pipeline is a model for the multi-frequencyangular power spectrum C`(ν, ν ′), from which a realisation of large-scale structure at the relevant redshifts is generated, consisting of mapsof the matter density contrast δm and gravitational potential φ. Thefinal step is to evaluate the cross-correlation function of the correlatedsimulations.• An outlier pixel mask is applied, but it is not derived from the simulations.Instead, the same mask as was used on the data is employed.The resulting catalogue of simulated CHIME spectra thus contains no foregroundresiduals, but does reflect the masking (of frequencies and pixels), as well as theeffect of the delay filter that is present in the real data. This is exactly what isneeded to produce a signal template for the measurement that captures the effect ofthe filtering performed on the data. Figure 4.13 summarises the simulation pipelinesteps in a flowchart diagram.The template is generated by evaluating the cross-correlation function of the simu-lated CHIME data and Lyman-α forest catalogue, ξ(∆ν), as defined in Equation (4.5).Figure 4.14 presents templates derived for a few different delay cuts, showing how thesignal is suppressed as more data are excised. The shape of the correlation function issharply peaked at zero separation, reflecting the fact that the larger scales are thrown11330 20 10 0 10 20 30 (MHz)1.00.50.0Surface brightness (arbitrary units)1e 6cut=100 nscut=200 nscut=300 nsFigure 4.14: Cross-correlation functions evaluated on the simulated 21 cm andLyman-α forest observations, for the X polarisation. These show thecorrelation as a function of frequency separation ∆ν in the CHIMEband averaged across all lines-of-sight in the eBOSS catalogue. Dif-ferent choices of delay cut τcut for the high-pass delay filter are shownas separate curves. As the delay cut is increased, so is the amount of21 cm signal excised, resulting in a suppression of the cross-correlationamplitude.away by the high-pass delay filter. It also demonstrates the characteristic negativecorrelation of the 21 cm and Lyman-α fields.4.7 ResultsThe cross-correlation functions of CHIME and the eBOSS Lyman-α forest (evaluatedas defined in Equation 4.5) are presented in Figure 4.15, where the two polarisationshave been combined to produce total intensity (they are shown separately in Fig-ure 4.16). An estimate of the background from noise and residual foregrounds wasderived by cross-correlating the CHIME data against permutations of the Lyman-αforest spectra, i.e. exchanging the spectrum measurements between lines of sight sothat they are uncorrelated with the CHIME observations but maintain exactly thesame distribution on the sky (this is described in more detail in Section 4.8). Themeasurements show a clear excess over the background around the zero-lag bin. Alsoshown in this figure is a signal template derived from the result of cross-correlating the11430 20 10 0 10 20 304202Surface brightness (arbitrary units) 1e 7Stokes ITemplate fit30 20 10 0 10 20 30 (MHz)2.50.02.5resid/Figure 4.15: (Top) Cross-correlation function of CHIME and eBOSS Lyman-αforest data, for the combined X and Y polarisations. An estimate ofthe standard deviation of the background is plotted as a shaded region.The black dashed line is a template derived from simulations with anamplitude fit to the data. Note that the full range of measured ∆νhas been cropped to focus on the region where signal is present abovethe noise. (Bottom) Residuals normalised by the estimated backgroundlevel.simulated datasets described in the previous section. Only the amplitude of the tem-plate was fit to the data, as explained below, and it appears to be broadly consistentwith the measurement.4.7.1 Detection SignificanceTo assess the significance of the detection, I derived a signal template from the simu-lations. As discussed in the previous section, synthetic CHIME and Lyman-α forestdata were generated without adding noise and run through the same filtering andcross-correlation analysis. The resulting correlation functions, the dashed lines inFigure 4.16, were then fit to the cross-correlation measurements by varying a sin-gle amplitude parameter and minimising the χ2. Finding the ∆χ2 = 1 bounds inparameter space provides an estimate for the standard deviation of the amplitudeconstraint. The fitted template and residuals are presented in Figure 4.15 for thecombined polarisations. The significance of the detection is quoted as the amplitude115of the template in units of its standard deviation from the fit in Table 4.1.Table 4.1: Detection significance of correlation function template fits. Thecolumns list the χ2 for the best-fit amplitude parameter; the reduced χ2;the χ2 difference from the null hypothesis; and the implied significance.The number of degrees of freedom used to compute χν is 254.measurement χ2 χ2ν ∆χ2null significanceX pol 217.1 0.85 40.2 6.4σY pol 295.8 1.16 42.8 6.5σcombined 238.5 0.94 82.7 9.1σ4.7.2 Comparison to FFT Cross-correlation MethodAs mentioned in Section 4.5.7, I considered an alternate method for evaluating thecross-correlation function of the CHIME and eBOSS spectra, using the FFT to com-pute the power spectrum in Fourier space, then transforming back to obtain thecorrelation function. In detail, the first step is to compute the transform of the fieldsδx (where x labels either the Lyman-α or HI field)δ˜x[n] =N∑je−2piijn/Nδx[j] , (4.57)followed by estimating the power spectrum averaged across all spectraP [n] =1MM∑iδiHI[n]δiα[n] . (4.58)Finally the correlation function is simply the inverse transform of the power spectrumξ[j] =1NN∑ne2piijn/NP [n] . (4.59)The correlation functions for each polarisation evaluated using this method are pre-sented in Figure 4.16, along with their associated detection significance levels.Although this procedure is nearly equivalent to the direct method presented inSection 4.3.2, there are a few important differences. First, the FFT method implicitlyassumes that the signal is periodic, and thus is equivalent to a direct correlationfunction computation where the sum wraps at the edges. This results in an equalnumber of samples contributing to every correlation function bin, which is why the11640 20 0 20 405.02.50.02.51e 75.0 , 2=1.05X polarisation (FFT x-corr)Measurement from CHIME and eBOSSSimulated signal40 20 0 20 405.02.50.02.51e 74.3 , 2=1.12Y polarisation (FFT x-corr)40 20 0 20 405.02.50.02.51e 76.3 , 2=0.85X polarisation (real space x-corr)40 20 0 20 40 (MHz)5.02.50.02.51e 76.6 , 2=1.16Y polarisation (real space x-corr)Surface brightness (arbitrary units)Figure 4.16: Measured cross-correlation functions for each polarisation, calcu-lated using either the FFT method or directly in real-space. The dashedcurves are the corresponding template fits, and the reduced χ2ν as wellas the detection significance are printed near the centre of each panel.Also note that the full range of frequency separations ∆ν has been in-cluded in this figure.variance plotted as a shaded region in Figure 4.16 is approximately uniform in theFFT cases. However, when I compute the direct method, I do not wrap the spectra(we do not expect the physical correlation to be periodic) so the variance in theestimate grows as it nears the edges of the range, where there are fewer samplescontributing to the sum.Second, in Equation 4.58 the inverse-variance noise weights are not used in theaverage. Propagating these weights (and the associated mask) to Fourier space is notas straightforward as is using them in the direct correlation function computation.This is the main disadvantage of the FFT method, and why I chose the direct methodinstead, despite it being more computationally expensive. The use of the noise weightsin the direct estimate is reflected in the better signal-to-noise obtained from the117template fit to those correlation functions.4.8 ValidationIn this section I describe a few checks that were performed to validate the detec-tion and estimate the level of background fluctuations (due to noise and residualcontamination).4.8.1 Lyman-α Forest PermutationsIn order to characterise the background fluctuations in the cross-correlation estimates,it is useful to generate a catalogue of Lyman-α forest measurements with the sameproperties as the eBOSS one, but which is uncorrelated with the CHIME data. Noise,foreground residuals, and any other systematic effects that are not correlated with thesignal will remain present in the resulting correlation functions, and their power can bemeasured. I do this by drawing random permutations of the Lyman-α forest spectra,while keeping the sample of line-of-sight directions fixed. 1000 such permutationswere generated and correlation functions evaluated for all of them. The standarddeviation over permutations at every spectral separation is shown as the shaded regionin Figure 4.15.By construction, this method preserves the distribution of background quasarson the sky and in redshift separately, but if there exist correlations between thesetwo axes in the selection function, they will be lost. I check for the presence of sucha correlation by binning the background quasars by their angular position onto aHEALPix (Górski et al., 2005) grid of Nside = 32 and computing the average redshiftin each bin. The standard deviation of the average redshifts across non-empty binsis ∆z < 0.08, or about 3%. For comparison, I shuffled the positions of the quasars,and repeated the binning calculation, resulting in a typical standard deviation of justunder 3%. There does not appear to be more variation in the redshifts with locationon the sky than would be expected for an uncorrelated distribution and I concludethat spatial selection effects are unlikely to be important.4.8.2 Even/Odd DaysThe full set of days of CHIME data that went into the sidereal stack were dividedinto two sets according to an even/odd split in chronological order, and each wasstacked in the same way as the full set (see Section 4.5.2). The 21 cm and foregroundsignals are the same on every sidereal day, so they are expected to be common to the11830 20 10 0 10 20 302021e 7 X polarisationeven - oddnoise realisations30 20 10 0 10 20 30 (MHz)2021e 7 Y polarisationSurface brightness (arbitrary units)Figure 4.17: Correlation function evaluated on the difference of even and odddays, with the standard deviation derived from noise realisations as ashaded regioneven and odd stacks, but thermal noise and RFI should be uncorrelated between thetwo. By taking the difference (and dividing by 2 to match the normalisation of thefull average), the common foreground residuals should cancel and what is left can becompared to our expectation for thermal noise. Any remaining excess over the noisecan be attributed to day-to-day variations, e.g. RFI or changes to the instrumentalresponse that were not entirely captured in the daily calibration.The correlation functions evaluated on the even/odd difference are shown in Fig-ure 4.17, along with the expected noise level. Not surprisingly, there is no signalaround the zero lag bin, but the background fluctuations appear to have excessnoise. To study the source of this excess noise, I generated 1000 permutations ofthe Lyman-α forest spectra and used them to evaluate the variance of the even-odddifference correlation function, using both the real CHIME data and a realisation ofCHIME noise. In all cases, the same set of permutations was used. Figure 4.18 showsthe variance across permutations for these two cases along with the regular data. Theratio between the variance in the regular analysis and the even-odd difference givesan idea for the contribution of residual foregrounds to the background fluctuationsin the correlation function. These account for about twice as much variance as theeven-odd difference, which in turn is about 3–5 times as large as the noise, indicatingthat there is significant contamination from sources that vary day-to-day. Note how-ever that the variance is skewed by significant non-Gaussian tails in the distributions11940 20 0 20 4010 1510 1410 13X polfulleven - oddnoise40 20 0 20 40 (MHz)10 1510 1410 13Y polvariance across 1000 permutationsFigure 4.18: Variance computed across correlation functions estimated on 1000permutations of the Lyman-α forest spectra. In blue is the full CHIMEdata, orange is the difference of even and odd days and green is a noiserealisation. Note that the minimum does not occur at zero separationbecause the maximal overlap of the highest weighted regions in theCHIME and eBOSS spectra is not achieved there.of the pixel values.Figure 4.19 shows a histogram of the pixels in the filtered and masked maps acrossall frequencies that were used in the 400–500MHz band, within a region of the skyrestricted roughly to the eBOSS NGC field. In the even-odd difference, the distribu-tion is closer to Gaussian for values . 2σ, but large tails remain. For comparison,similar histograms were generated for ringmaps in the upper band, 585–800MHz,displayed in Figure 4.20. These are much closer to Gaussian for . 3σ, especially inthe even-odd difference.The analysis of the difference of even and odd days in the lower band suggeststhat:1. there is a large excess of variance in the maps above the thermal noise (by afactor of 6− 10);2. this excess is not fully accounted for by features that are repeated day-to-day likecelestial foregrounds (the variance is only ∼halved in the even-odd difference);120Figure 4.19: Histogram of pixel values in filtered and masked ringmaps acrossall frequencies used in the lower band, 400–500MHz. The pixels arenormalised to the expected noise standard deviation, and a Gaussiancurve with a matching amplitude is plotted in orange for comparison.Rows show the two polarisations, and columns show the difference ofeven and odd days and the average of all days. The effect of the outliermask can be seen as the sharp cutoff at 6σ.3. the excess above the noise is non-Gaussian, with large tails in the distributionof map pixels;4. the distribution is significantly less Gaussian in the lower band than it was inthe upper band, especially in the even-odd difference.A possible explanation for this behaviour is that the lower band is more contaminatedthan the upper band by unflagged RFI – which does not repeat day-to-day andcontributes in a non-Gaussian-distributed way. Another contributing effect couldalso be errors in the daily calibration, which would have a similar signature. Bothof these factors have been reported in other internal CHIME analyses of the lowerband. Identifying the source of this contamination and removing it will likely require121Figure 4.20: Histogram of pixel values in filtered and masked ringmaps acrossall frequencies used in the upper band, 585–800MHz. The pixels arenormalised to the expected noise standard deviation, and a Gaussiancurve with a matching amplitude is plotted in orange for comparison.Rows show the two polarisations, and columns show the difference ofeven and odd days and the average of all days. The effect of the outliermask can be seen as the sharp cutoff at 6σ.another look at the individual days, before they are averaged into a sidereal stack.4.9 ConclusionThis chapter presents the first reported detection of 21 cm emission at redshiftsz > 1.5, obtained by cross-correlating CHIME data from 88 days in the 400–500 MHzband with the DR16 eBOSS measurements of the Lyman-α forest. To do so, I ex-tended the analysis methods described in the first CHIME detection paper (CHIMECollaboration et al., 2022b) to the lower quarter of the CHIME frequency band, anddeveloped a cross-correlation method to combine this data with the Lyman-α forestas a tracer of large-scale structure.122This demonstrates that CHIME is able to detect emission from cosmological HIin the entire redshift range available to it, but it also highlights the difficulty ofseparating this signal from the extremely bright foregrounds. Even in cross-correlationand after aggressive delay filtering, the excess variance from foregrounds was foundto be about 6–10 times above thermal noise. This level of contamination neverthelessallowed for a 9σ detection. The shape of the measured correlation function appearsto be dominated by the effect of the foreground filter, and I have not attempted tointerpret it in this work. The amplitude of the cross-correlation should in principletell us something about the physics of the Lyman-α forest, 21 cm brightness and large-scale structure, but connecting this measurement with the relevant theory will requiremore sophisticated modelling of the Lyman-α and 21 cm signals. Additions to themodelling of the Lyman-α field should include redshift-space distortions, metal lines,high-density absorber systems, and biases due to the quasar spectra fits. Anotherimprovement would be to model the Lyman-α and 21 cm signals jointly, rather thanindependently as was done here. Including more realistic instrumental effects in theCHIME simulations should also be investigated. This will be the subject of futurework.Improvements to the CHIME calibration and analysis pipelines, as well as theinclusion of vastly more data that has yet to be processed, will eventually lead to anincrease in the constraining power of the analysis. In the short term, these mightreduce the number of frequencies that need to be flagged and enable a less restrictivedelay cut, allowing for information from larger scales to be recovered. The Lyman-αforest is known to be correlated up to and beyond the BAO scale, so we can expectthat measuring its cross-correlation with 21 cm emission beyond the very small scalesmanaged here will yield useful information. Nevertheless, the cross-correlation signalthat was detected is already a valuable benchmark for the effectiveness of the 21 cmcosmology analysis of CHIME data, and in the longer term has the potential to be apowerful probe of cosmology.123Chapter 5ConclusionIn this thesis I presented research that I performed towards commissioning and cal-ibrating CHIME, as well as an analysis of the data that resulted in a detection ofcosmological 21 cm signal in cross-correlation with another tracer. Although com-posed of many disparate tasks, the work I performed over the years I spent as aPhD student at University of British Columbia (UBC) was directed towards makingCHIME successful, and amounts to a significant contribution to the experiment, aswell as enabling other science cases being investigated using the telescope we built.The two detections in cross-correlation with external surveys that have now beenreported at the lower and upper ends of the redshift range accessible to CHIME havedemonstrated that intensity mapping with the 21 cm line is a viable probe of cos-mology and have established CHIME as a pathfinder in this nascent field. The nextstep for the analysis is to attempt a detection using CHIME data alone, which wouldconstitute another milestone for the field. Achieving this will require further progressin the precise calibration of the instrumental response, as well as new methods ofanalysis in order to identify and remove contamination from the bright foregrounds.Current efforts within the collaboration are focused on the analysis of the high-delayregion of the data, where foreground contamination is less significant. These mea-surements of the power spectrum are beginning to approach the expected level of the21 cm signal.These challenges are shared by the other experiments in the 21 cm intensity map-ping field, including those targeting higher redshifts and the epoch of reionisation,which all stand to benefit from the continued development and sharing of methods ofanalysis. Those being reported in the literature show steady progress (as summarisedin the Introduction; Section 1.1.4). Notably, the recent publication of a detectionfrom 21 cm observations alone (Paul et al., 2023) is very encouraging.124For CHIME, continued development of the processing and analysis for the holo-graphic beam measurements will be important to refine and expand the beam modelof the telescope that is essential to effectively remove foregrounds. The data productsproduced by the processing pipeline that I described in this thesis will facilitate suchanalyses, and extensions to the pipeline will refine the beam measurements thus de-rived. Gated holography of pulsars, especially the polarised ones, is another valuablesource of information about the beam, and the only avenue currently available forcalibrating its polarisation response. Using the gated acquisition mode that I devel-oped, observations for many pulsars have been taken, but they have yet to be putto use. The analysis method that I presented in the thesis explains how to extractpolarised beam measurements from this dataset, and it should be put to the test inthe future.The cross-correlation I performed between CHIME data and measurements of theLyman-α forest yielded the highest-redshift detection of 21 cm emission published todate, at an average redshift of z¯ = 2.3. Future work on this analysis will be directedtowards interpreting the measurement and its implications for cosmology and thephysics of the inter-galactic medium. For instance, using the well-studied Lyman-αabsorption signal as a reference has the potential to shed light on the distribution andbias of the 21 cm field at these redshifts, for which measurements are presently lacking.However, this will require a more careful modelling approach than was presented inthis thesis, in order to capture more of the physical and instrumental effects that areat play. There is also a lot of potential to improve the quality of the measurement withrefinements to the data processing pipeline, especially with regard to the flagging ofcontaminated data. There is evidence that the lower part of the CHIME band used inthis analysis is more heavily contaminated than the upper part that was consideredin the lower-redshift cross-correlation result, in particular by sporadic and fainterRFI. This was suggested by the analysis of a difference of even and odd partitionsof the data that showed excess variance at the high delays relevant to this analysis.Correlations up to larger scales are expected to be present in the cross-correlation ofthese two tracers, and reining in the aggressive filtering that had to be applied couldunlock information that would be valuable for a cosmological interpretation.CHIME has also made an impact outside of cosmology. Fast radio bursts (FRBs)– short, apparently random, and extremely bright bursts of radiation at radio wave-lengths – were first discovered in 2007, and up to a few years ago very little wasknown about their properties due to the impossibility of predicting where and whento look for them. That is until CHIME was outfitted with a parallel backend to125search for fast radio burst (FRB)s, which has resulted in the publication of thousandsof FRB detections (CHIME/FRB Collaboration et al., 2021), upending the field andprompting a flurry of interest in their origin and emission mechanism. CHIME is alsoequipped with a separate backend to monitor known pulsars and search for new ones(CHIME/Pulsar Collaboration et al., 2020), and is an element in the NANOGrav Pul-sar Timing Array. Finally, CHIME is one of the instruments that will participate inthe Global Magneto-ionic Medium Survey (GMIMS), which aims to observe the MilkyWay over the entire sky and a broad band in frequency in order to perform Faradaytomography and map the magnetic fields coursing through it (Wolleben et al., 2021).My contributions and those of the other CHIME team members to the constructionand operation of the telescope enabled the realisation of these additional projects.126BibliographyAbdalla, E., Abellán, G. F., Aboubrahim, A., et al. 2022, Journal of High EnergyAstrophysics, 34, 49, doi: 10.1016/j.jheap.2022.04.002 → page 11Alam, S., Aubert, M., Avila, S., et al. 2021, Phys. Rev. D, 103, 083533,doi: 10.1103/PhysRevD.103.083533 → pages 11, 13Amiri, M., Bandura, K., Boskovic, A., et al. 2022, ApJ, 932, 100,doi: 10.3847/1538-4357/ac6b9f → page 46Anderson, C. J., Luciw, N. J., Li, Y. C., et al. 2018, MNRAS, 476, 3382,doi: 10.1093/mnras/sty346 → page 16Baars, J. W. 2007, The Paraboloidal Reflector Antenna in Radio Astronomy andCommunication (Springer) → page 46Bandura, K., Cliche, J. F., Dobbs, M. A., et al. 2016, Journal of AstronomicalInstrumentation, 5, 1641004, doi: 10.1142/S225117171641004X → page 28Baumann, D. 2022, Cosmology (Cambridge University Press),doi: 10.1017/9781108937092 → pages 2, 6, 9Bennett, J. C., Anderson, A. P., McInnes, P. A., & Whitaker, A. J. T. 1976, IEEETransactions on Antennas and Propagation, 24, 295,doi: 10.1109/TAP.1976.1141354 → page 46Berger, P., Newburgh, L. B., Amiri, M., et al. 2016, in Society of Photo-OpticalInstrumentation Engineers (SPIE) Conference Series, Vol. 9906, Ground-basedand Airborne Telescopes VI, ed. H. J. Hall, R. Gilmozzi, & H. K. Marshall,99060D, doi: 10.1117/12.2233782 → page 47Bi, H., & Davidsen, A. F. 1997, 479, 523, doi: 10.1086/303908 → page 110Bull, P., Ferreira, P. G., Patel, P., & Santos, M. G. 2015, ApJ, 803, 21,doi: 10.1088/0004-637X/803/1/21 → page 108Carucci, I. P., Villaescusa-Navarro, F., & Viel, M. 2017, Journal of Cosmology andAstroparticle Physics, 2017, 001, doi: 10.1088/1475-7516/2017/04/001 → page 77Castorina, E., & Villaescusa-Navarro, F. 2017, MNRAS, 471, 1788,doi: 10.1093/mnras/stx1599 → page 107127Chang, T.-C., Pen, U.-L., Bandura, K., & Peterson, J. B. 2010, Nature, 466, 463,doi: 10.1038/nature09187 → page 16CHIME Collaboration, Amiri, M., Bandura, K., et al. 2022a, The AstrophysicalJournal Supplement Series, 261, 29, doi: 10.3847/1538-4365/ac6fd9 → pagesvi, 22, 23, 24, 25, 26, 28, 30, 31, 32, 33, 34, 48, 80, 90, 92—. 2022b, doi: 10.48550/arXiv.2202.01242 → pagesvi, 17, 24, 31, 46, 76, 79, 80, 82, 88, 90, 91, 93, 94, 100, 101, 102, 105, 106, 122, 133CHIME/FRB Collaboration, Amiri, M., Andersen, B. C., et al. 2021, ApJS, 257, 59,doi: 10.3847/1538-4365/ac33ab → page 126CHIME/Pulsar Collaboration, Amiri, M., Bandura, K. M., et al. 2020, arXive-prints, arXiv:2008.05681. https://arxiv.org/abs/2008.05681 → pages 38, 40, 126Cieplak, A. M., & Slosar, A. 2016, 2016, 016, doi: 10.1088/1475-7516/2016/03/016→ page 111Cohen, A. S., Lane, W. M., Cotton, W. D., et al. 2007, The Astronomical Journal,134, 1245, doi: 10.1086/520719 → page 99Condon, J. J., Cotton, W. D., Greisen, E. W., et al. 1998, AJ, 115, 1693,doi: 10.1086/300337 → page 100Cosmic Visions 21 cm Collaboration, Ansari, R., Arena, E. J., et al. 2018, arXive-prints, arXiv:1810.09572. https://arxiv.org/abs/1810.09572 → page 106Crighton, N. H. M., Murphy, M. T., Prochaska, J. X., et al. 2015, MNRAS, 452,217, doi: 10.1093/mnras/stv1182 → page 109Croft, R. A. C., Weinberg, D. H., Katz, N., & Hernquist, L. 1998, ApJ, 495, 44,doi: 10.1086/305289 → page 77Cunnington, S., Li, Y., Santos, M. G., et al. 2023, MNRAS, 518, 6262,doi: 10.1093/mnras/stac3060 → page 16Deng, M. 2020, PhD thesis, University of British Columbia,doi: http://dx.doi.org/10.14288/1.0395494 → page 28Deng, M., & Campbell-Wilson, D. 2014, 2014 16th International Symposium onAntenna Technology and Applied Electromagnetics (ANTEM),doi: 10.1109/antem.2014.6887670 → page 28Dodelson, S., & Schmidt, F. 2020, Modern Cosmology,doi: 10.1016/C2017-0-01943-2 → pages 2, 4du Mas des Bourboux, H., Rich, J., Font-Ribera, A., et al. 2020, The AstrophysicalJournal, 901, 153, doi: 10.3847/1538-4357/abb085 → pages 77, 78128Edwards, R. T., Hobbs, G. B., & Manchester, R. N. 2006, Monthly Notices of theRoyal Astronomical Society, 372, 1549, doi: 10.1111/j.1365-2966.2006.10870.x →page 40Ewall-Wice, A., Kern, N., Dillon, J. S., et al. 2021, MNRAS, 500, 5195,doi: 10.1093/mnras/staa3293 → pages 87, 89Farr, J., Font-Ribera, A., Bourboux, H. d. M. d., et al. 2020, 2020, 068,doi: 10.1088/1475-7516/2020/03/068 → pages 110, 111Górski, K. M., Hivon, E., Banday, A. J., et al. 2005, ApJ, 622, 759,doi: 10.1086/427976 → pages 106, 118Gould, R. J. 1994, ApJ, 423, 522, doi: 10.1086/173830 → page 108Gregory, P. C., Scott, W. K., Douglas, K., & Condon, J. J. 1996, ApJS, 103, 427,doi: 10.1086/192282 → page 100HERA Collaboration, Abdurashidova, Z., Adams, T., et al. 2023, ApJ, 945, 124,doi: 10.3847/1538-4357/acaf50 → page 17Hu, W. 1995, arXiv e-prints, astro, doi: 10.48550/arXiv.astro-ph/9508126 → page 12IEEE. 2019, IEEE Std 754-2019 (Revision of IEEE 754-2008), 1,doi: 10.1109/IEEESTD.2019.8766229 → page 33Jackson, J. C. 1972, Mon. Not. Roy. Astron. Soc., 156, 1P,doi: 10.1093/mnras/156.1.1P → page 106Kaiser, N. 1987, MNRAS, 227, 1, doi: 10.1093/mnras/227.1.1 → page 81Kaiser, N. 1987, Mon. Not. Roy. Astron. Soc., 227, 1 → page 106Karamanis, M., & Beutler, F. 2021, arXiv e-prints, arXiv:2106.06331.https://arxiv.org/abs/2106.06331 → page 106Lewis, A., Challinor, A., & Lasenby, A. 2000, ApJ, 538, 473, doi: 10.1086/309179 →page 106Li, L.-C., Staveley-Smith, L., & Rhee, J. 2021, Research in Astronomy andAstrophysics, 21, 030, doi: 10.1088/1674-4527/21/2/30 → page 16Liu, A., & Shaw, J. R. 2020, PASP, 132, 062001, doi: 10.1088/1538-3873/ab5bfd →pages 17, 20, 21Liu, A., & Tegmark, M. 2011, in American Astronomical Society Meeting Abstracts,Vol. 218, American Astronomical Society Meeting Abstracts #218, 220.03 →page 16Masui, K., Amiri, M., Connor, L., et al. 2015, Astronomy and Computing, 12, 181,doi: https://doi.org/10.1016/j.ascom.2015.07.002 → page 33129Masui, K. W., Switzer, E. R., Banavar, N., et al. 2013, ApJ, 763, L20,doi: 10.1088/2041-8205/763/1/L20 → page 16Mead, A. J., Brieden, S., Tröster, T., & Heymans, C. 2021, MNRAS, 502, 1401,doi: 10.1093/mnras/stab082 → page 106Morales, M. F., Hazelton, B., Sullivan, I., & Beardsley, A. 2012, ApJ, 752, 137,doi: 10.1088/0004-637X/752/2/137 → page 86Narayan, R. 1992, Philosophical Transactions of the Royal Society of London SeriesA, 341, 151, doi: 10.1098/rsta.1992.0090 → page 65Paciga, G. 2013, PhD thesis, University of Toronto, Canada → page 72Particle Data Group, Zyla, P. A., Barnett, R. M., et al. 2020, Progress of Theoreticaland Experimental Physics, 2020, doi: 10.1093/ptep/ptaa104 → page 10Paul, S., Santos, M. G., Chen, Z., & Wolz, L. 2023, arXiv e-prints,arXiv:2301.11943, doi: 10.48550/arXiv.2301.11943 → pages 17, 124Pen, U.-L., Staveley-Smith, L., Peterson, J. B., & Chang, T.-C. 2009, MNRAS, 394,L6, doi: 10.1111/j.1745-3933.2008.00581.x → page 16Penrose, R. 1955, Mathematical Proceedings of the Cambridge PhilosophicalSociety, 51, 406–413, doi: 10.1017/S0305004100030401 → page 88Percival, W. J., & White, M. 2009, MNRAS, 393, 297,doi: 10.1111/j.1365-2966.2008.14211.x → page 81Perley, R. A., & Butler, B. J. 2017, ApJS, 230, 7, doi: 10.3847/1538-4365/aa6df9 →pages 30, 53Perlmutter, S., Aldering, G., Goldhaber, G., et al. 1999, ApJ, 517, 565,doi: 10.1086/307221 → page 9Peterson, J. B., Aleksan, R., Ansari, R., et al. 2009, in astro2010: The Astronomyand Astrophysics Decadal Survey, Vol. 2010, 234, doi: 10.48550/arXiv.0902.3091→ page 15Planck Collaboration, Ade, P. A. R., Aghanim, N., et al. 2014, A&A, 571, A1,doi: 10.1051/0004-6361/201321529 → pages vi, 12Planck Collaboration, Aghanim, N., Akrami, Y., et al. 2020, A&A, 641, A6,doi: 10.1051/0004-6361/201833910 → pages vi, 3, 9, 10, 12Renard, A., Shaw, R., Ng, C., et al. 2021, Kotekan: A framework forhigh-performance radiometric data pipelines, 2021.11, Zenodo,doi: 10.5281/zenodo.5842660 → pages vi, 32130Rengelink, R. B., Tang, Y., de Bruyn, A. G., et al. 1997, A&AS, 124, 259,doi: 10.1051/aas:1997358 → page 99Riess, A. G., Filippenko, A. V., Challis, P., et al. 1998, AJ, 116, 1009,doi: 10.1086/300499 → page 9Santos, M. G., Cooray, A., & Knox, L. 2005, ApJ, 625, 575, doi: 10.1086/429857 →page 86Sarkar, D., & Bharadwaj, S. 2019, MNRAS, 487, 5666, doi: 10.1093/mnras/stz1691→ page 107Scoccimarro, R. 2004, Phys. Rev. D, 70, 083007, doi: 10.1103/PhysRevD.70.083007→ page 107Scott, D., & Rees, M. J. 1990, MNRAS, 247, 510 → page 15Scott, P. F., & Ryle, M. 1977, MNRAS, 178, 539, doi: 10.1093/mnras/178.4.539 →page 46Seljak, U. 2012, 2012, 004, doi: 10.1088/1475-7516/2012/03/004 → page 111Shaw, J. R., Foreman, S., Nitsche, R., et al. 2020a, radiocosmology/driftscan:20.5.0, v20.5.0, Zenodo, Zenodo, doi: 10.5281/zenodo.5949193 → pages vi, 89Shaw, J. R., Nitsche, R., Foreman, S., & Kefala, A. 2020b, radiocosmology/cora:20.5.0, v20.5.0, Zenodo, Zenodo, doi: 10.5281/zenodo.5949195 → page 89Shaw, J. R., Sigurdson, K., Sitwell, M., Stebbins, A., & Pen, U.-L. 2015, PhysicalReview D, 91, 083514, doi: 10.1103/PhysRevD.91.083514 → page 109Shaw, J. R., Masui, K., Nitsche, R., et al. 2020c, radiocosmology/caput: 20.10.0,v20.10.0, Zenodo, Zenodo, doi: 10.5281/zenodo.5846375 → pages vi, 89Shaw, R., Nitsche, R., Siegel, S. R., et al. 2020d, radiocosmology/draco: 20.10.0,v20.10.0, Zenodo, Zenodo, doi: 10.5281/zenodo.5828373 → pages vi, 89Stairs, I. H., Thorsett, S. E., & Camilo, F. 1999, ApJS, 123, 627,doi: 10.1086/313245 → page 65Taylor, J., Denman, N., Bandura, K., et al. 2018, Journal of AstronomicalInstrumentation, 08, doi: 10.1142/S225117171940004X → page 29Tramonte, D., & Ma, Y.-Z. 2020, MNRAS, 498, 5916, doi: 10.1093/mnras/staa2727→ page 16van Straten, W. 2004, ApJS, 152, 129, doi: 10.1086/383187 → pages 70, 73Villaescusa-Navarro, F., Genel, S., Castorina, E., et al. 2018, ApJ, 866, 135,doi: 10.3847/1538-4357/aadba0 → pages 15, 106, 108131Weiland, J. L., Odegard, N., Hill, R. S., et al. 2011, ApJS, 192, 19,doi: 10.1088/0067-0049/192/2/19 → page 65Wheeler, J. A., & Ford, K. 1998, Geons, black holes, and quantum foam: A life inphysics → page 2Wiener, N. 1949, Extrapolation, Interpolation, and Smoothing of Stationary TimeSeries: With Engineering Applications (The MIT Press),doi: 10.7551/mitpress/2946.001.0001 → page 70Wiese, W. L., & Fuhr, J. R. 2009, Journal of Physical and Chemical ReferenceData, 38, 565, doi: 10.1063/1.3077727 → page 108Wolleben, M., Landecker, T. L., Douglas, K. A., et al. 2021, The AstronomicalJournal, 162, 35, doi: 10.3847/1538-3881/abf7c1 → page 126Wolz, L., Pourtsidou, A., Masui, K. W., et al. 2022, MNRAS, 510, 3495,doi: 10.1093/mnras/stab3621 → page 16132Appendix AReverse Lanczos InterpolationIn this appendix, I describe the reverse Lanczos method for interpolating a noisysignal onto a new coordinate grid with coarser sampling. This is the same methodthat is described in Section 3.2.6 of CHIME Collaboration et al. (2022b). First,consider the direct Lanczos resampling scheme, which uses a truncated sinc kernel tointerpolate between regularly spaced discrete samples, with another sinc envelope forsmooth apodization:La(xi, yj) =sinc(xi − yj) sinc(xi−yja), −a < xi − yj < a,0, otherwise .(A.1)Here a is the width of the kernel, the xi are the positions of the interpolated pointsand yj are the coordinates of the original sequence. To interpolate onto a new grid,this kernel is convolved with the data, an operation that can be represented as amatrix multiplication, as the i, j indices in the equation suggest:v′ = Lv. (A.2)In the limit where a → ∞, the Lanczos kernel becomes a sinc function over aninfinite domain resulting in the exact recovery of a band-limited signal as long asit was sampled at least at the Nyquist rate. Increasing the width parameter betterapproximates this limit at the expense of extra computation.To interpolate a noisy dataset like a CHIME visibility onto a downsampled grid,the direct Lanczos scheme may not be the best choice since it will not combine datapoints in such a way as to beat down the noise. For instance, if we just chose tocut the sample cadence in half, keeping the same grid, applying the Lanczos method133would throw away half the data points, not average them, since every adjacent pointwould fall onto a node in the sinc function. Instead, we can cast this as an inverseproblem, accounting for the noise covariance.In this scheme, the measured signal d is modelled in terms of the true signal vdefined at the desired coordinates and interpolated to the measured points by theLanczos kernel L, and noised = Lv + n. (A.3)The problem is to find an estimate vˆ of the true signal v that minimises the squarederror(d− Lvˆ)† (d− Lvˆ) . (A.4)Assuming a solution in the form of a linear filter – vˆ = Wd – and that we havean expectation for the signal and noise covariances S =〈vv†〉and N =〈nn†〉, thesolution is the Wiener filter,W = SL†(LSL† +N)−1, (A.5)andC =〈vˆvˆ†〉=(S−1 + L†N−1L)−1(A.6)is the covariance for the interpolated data.The Wiener filter defined above assumes knowledge of the signal covariance, whichwe may not have a priori. In practice however, the filter is not very sensitive to thischoice, and using a covariance proportional to the identity works well. The amplitudeof the signal relative to the noise is a parameter that can be tuned to modulate theaggressiveness with which it suppresses noisy modes depending on the application.134Appendix BSpectral IncoherenceA brief derivation of the equivalence of stationarity and spectral incoherence for arandom process is given here. Consider a random signal in time, and its conjugate infrequency,s(t) =∫dνe2piiνts(ν) . (B.1)We say it is statistically stationary if its auto-correlation function is independent oftime:C(t, t′) = 〈s(t)s∗(t′)〉 = C(t′ − t) , (B.2)where the angle brackets indicate an ensemble average. Similarly, the correlationbetween frequencies is on averageP (ν, ν ′) = 〈s(ν)s∗(ν ′)〉 (B.3)=∫dtdt′e−2piiνte2piiν′t′ 〈s(t)s∗(t′)〉 (B.4)=∫dtdt′e−2piiνte2piiν′t′C(t′ − t) , (B.5)where we have used the definition of C(t, t′) and stationarity C(t, t′) = C(t′ − t).Changing variables to τ = t′ − t for the inner integral,P (ν, ν ′) =∫dtdτe2piiν′τe2pii(ν′−ν)tC(τ) (B.6)=∫dτe2piiν′τδ(ν ′ − ν)C(τ) (B.7)= δ(ν ′ − ν)∫dτe2piiν′τC(τ) (B.8)≡ δ(ν − ν ′)P (ν) . (B.9)135This reflects the Wiener-Khinchin theorem, while emphasizing the lack of correlationsbetween frequencies that is a corollary.136Appendix CSpatial Aliasing in InterferometricMeasurementsFigure 2.3 in Chapter 2 shows a map of the sky generated from CHIME data, wherestructures near the top and bottom edges show up again as aliased copies on the otherside of the map. This phenomenon arises due to the interferometric measurementprocess, in which the direction to a source of radiation is inferred within the broadfield of view of an individual antenna from the delay in its reception by anotherantenna separated from the first by a known distance, as illustrated in Figure C.1.At an observing wavelength λ and a baseline separation b, this delay isντ =bλcos θ , (C.1)where the delay τ is multiplied by the frequency ν to give the number of cycles, and θis the direction to the source above the horizon. From horizon to horizon, delays canspan −b/λ < ντ < b/λ. The periodicity of the radiation means that delays cannotbe distinguished if they are separated by an integer number of cycles: ντ ∼ ντ + n.Aliasing occurs when b/λ > 1/2, in which case a measured ντ > 1/2 is confused witha source from a direction corresponding to ντ − 1 > −1/2, which is also above thehorizon. This is a manifestation of the Nyquist-Shannon sampling theorem, appliedto visibilities, which are measurements in aperture space at separations b/λ, conjugateto sky position, cos θ in this example. In order to measure the sky up to cos θ = 1,baselines need to be sampled with a shortest separation ∆b/λ < 1/2. For example,for CHIME the shortest baseline separation is about 30 cm, which is insufficient tosample the sky from horizon to horizon without aliasing for frequencies ν & 500MHz.This is apparent in Figure 2.3.137θbb cosθb cosθ - λFigure C.1: Diagram illustrating the situation described in this section. Radi-ation from a celestial source located an angle θ above the horizon (inred) is received by two feeds separated by a baseline distance b. Thedelay between the two signals is the projection of the source directiononto the baseline divided by the speed of light, b cos θ/c. If b < λ/2,it is possible for radiation from another direction above the horizon (inblue) to arrive with a delay different from the first by an integer numberof wavelengths, λ/c. These two thus cannot be distinguished and arealiased onto each other.138
Thesis/Dissertation
eng
Physics
Vancouver : University of British Columbia Library
Attribution-NonCommercial-NoDerivatives 4.0 International
Graduate
Hinshaw, Gary
A detection of cosmological 21 cm emission from CHIME in cross-correlation with the eBOSS Lyman-𝜶 forest
Text