TY - JOUR
AU - Chen, T
AB - ABSTRACT Cosmology has entered an era where the experimental limitations are not due to instrumental sensitivity but instead due to inherent systematic uncertainties in the instrumentation and data analysis methods. The field of neutral hydrogen (HI) intensity mapping (IM) is still maturing, however early attempts are already systematics limited. One such systematic limitation is 1/f noise, which largely originates within the instrumentation and manifests as multiplicative gain fluctuations. To date, there has been little discussion about the possible impact of 1/f noise on upcoming single-dish HI IM experiments such as BINGO (Baryonic Acoustic Oscillations from Integrated Neutral Gas Observations), FAST (Five hundred metre Aperture Spherical Telescope), or SKA (Square Kilometre Array). Presented in this work are Monte Carlo end-to-end simulations of a 30 d HI IM survey using the SKA–MID array covering a bandwidth of 950 and 1410 MHz. These simulations extend 1/f noise models to include not just temporal fluctuations, but also correlated gain fluctuations across the receiver bandpass. The power spectral density of the spectral gain fluctuations are modelled as a power law, and characterized by a parameter β. It is found that the degree of 1/f noise frequency correlation will be critical to the success of HI IM experiments. Small values of β (β < 0.25) or high correlation is preferred as this is more easily removed using current component separation techniques. Spectral index of temporal fluctuations (α) is also found to have a large impact on signal-to-noise ratio. Telescope slew speed has a smaller impact, and a scan speed of 1 deg s−1 should be sufficient for an HI IM survey with the SKA. instrumentation: spectrographs, methods: data analysis, methods: statistical, cosmology: observations, large-scale structure of Universe, radio lines: galaxies 1 INTRODUCTION Cosmology has entered an era where most cosmological parameters are known to the one per cent level using a combination of probes including the cosmic microwave background (CMB, Planck Collaboration et al. 2016a), galaxy redshift surveys (e.g. Eisenstein et al. 2005; Alam et al. 2015; Aubourg et al. 2015; Hinton et al. 2017), and Type Ia supernovae (e.g. Riess et al. 1998; Betoule et al. 2014). Another possible cosmological probe that is yet to be exploited is the hyperfine 21 cm spin-flip transition of neutral hydrogen (HI), which can be used to, for example, trace large-scale structure and the ionization history of the Universe (Pritchard & Loeb 2012). Given the abundance of hydrogen within the Universe one might expect HI emission to be a bright signal. However, given the expected lifetime of 107 yr for the upper transition state of the 21 cm line, HI emission in reality is extremely faint. To date, the detection of individual galaxies up to modest redshifts of z < 0.4 has proved extremely challenging using the current generation telescopes (Zwaan, van Dokkum & Verheijen 2001; Zwaan et al. 2004; Giovanelli et al. 2005; Verheijen et al. 2010; Catinella & Cortese 2015; Fernández et al. 2016). By the stacking of HI emission from many galaxies (Lah et al. 2009; Rhee et al. 2016), exploiting strongly lensed systems (e.g. Hunt, Pisano & Edel 2016), and measurements of damped Lyman-α systems (e.g. Noterdaeme et al. 2012; Neeleman et al. 2016) it is known that there is a vast quantity of HI present throughout the cosmic history of the Universe. It is this extremely faint HI signal and the potential cosmological information within that was the original science driver behind the Square Kilometre Array (SKA) proposal (Wilkinson 1991). HI emission as a probe of cosmological structure is considered extremely promising as it is expected to have far moreinformation content than the CMB. This is because tracing HI emission allows for the Universe to be divided in to redshift slices and the evolution of structure to be studied. HI cosmology is expected to yield both vastly improved constraints on parameters that describe the evolution of dark energy, neutrino masses, HI bias, HI mass density, and the astrophysical processes occurring within the Epoch of Reionization (EoR, Loeb & Zaldarriaga 2004; Santos, Cooray & Knox 2005; Sethi 2005; Fan, Carilli & Keating 2006; McQuinn et al. 2008; Loeb & Wyithe 2008; Pritchard & Loeb 2008). However, contemporary HI instrumentation lacks the sensitivity to conduct a traditional galaxy survey using the 21 cm line. Such a survey today would suffer from high levels of incompleteness, uncertainties in redshifts, and long-survey times. For these reasons, there has been a drive to find an alternative method of measuring the cosmic HI signal. The proposed alternative to HI galaxy surveys is known as the HI intensity mapping (IM) method. HI IM foregoes the detection of individual galaxies in favour of the total integral HI flux density over large cosmological volumes. The benefits of HI IM are no bias in galaxy selection, intrinsic knowledge of redshift, and a higher surface brightness relative to galaxy surveys. HI IM has the potential to be one of the most important tools available to understanding the evolution of large-scale structure, as it contains several orders of magnitude more information than the CMB (Battye, Davies & Weller 2004; Wyithe & Loeb 2008; Masui, McDonald & Pen 2010; Battye et al. 2013; Santos et al. 2014, 2015; Padmanabhan, Choudhury & Refregier 2015; Bull et al. 2015; Ma & Scott 2016). There have been several attempts at probing the post-reionization HI signal using the IM technique with existing telescopes at the Parkes (Pen et al. 2009; Anderson et al. 2017) and Green Bank (Chang et al. 2010; Masui et al. 2010, 2013; Shaw et al. 2014; Wolz et al. 2017) observatories. However, due to systematic limitations from the instrumentation and radio-frequency interference (RFI), the HI power spectrum was only detected in cross-correlation with optical surveys. There are a number of bespoke post-reionization epoch experiments planned such as the Baryonic Acoustic Oscillations from Integrated Neutral Gas Observations (Battye et al. 2013), Canadian Hydrogen Intensity Mapping Experiment (Bandura et al. 2014), Tianlai Telescope (Chen 2012), Five hundred metre Aperture Spherical Telescope (FAST, Smoot & Debono 2017; Bigot-Sazy et al. 2016), Hydrogen Intensity and Real-Time Analysis experiment (Newburgh et al. 2016), and a potential survey with the SKA–MID array (Santos et al. 2015). There has also been a recent proposal to perform an HI IM survey with MeerKAT (Santos et al. 2017). At lower frequencies, there are numerous phased array telescopes planning to probe the EoR such as the Precision Array for Probing the Epoch of Reionization (Parsons et al. 2010), Low Frequency Array (van Haarlem et al. 2013), the Giant Metrewave Radio Telescope (Paciga et al. 2011), Murchison Widefield Array (Bowman et al. 2013), and Hydrogen Epoch of Reionization Array (DeBoer et al. 2017). Similarly there are plans to perform IM using other emission lines such as CO (Lidz et al. 2011; Keating et al. 2015; Li et al. 2016; Keating et al. 2016), Lyman-α (Croft et al. 2016), and CII (Crites et al. 2014). The focus of this paper will be on the HI IM survey discussed in Bull et al. (2015) that proposes the use of the SKA–MID array as a collection of single-dish (SD) radio telescopes as opposed to a single interferometer. This is sometimes referred to as autocorrelation or total-power observations. The forecast for such a survey as described in Bull et al. (2015) suggests that the cosmological constraints on dark energy parameters at redshifts less than z < 1 will be comparable to upcoming Stage IV dark energy experiments such as Euclid (Albrecht et al. 2006). The claims in Bull et al. (2015) were further supported by simulations of such an SKA survey by Alonso et al. (2015), where it was demonstrated that existing methods for component separation are sufficient to remove astrophysical foregrounds in intensity. However, for real HI IM data, there will be challenges due to instrumental systematics such as standing waves, frequency-dependent asymmetric beams, or polarizations leakage. The current work will investigate the impact on the recovery of the HI angular power spectrum in the presence of a single instrumental systematic that is commonly referred to as 1/f noise. As many instrumental systematics are time dependent, determining the exact impact each has on the underlying HI signal requires end-to-end simulations of signal on the sky, through a receiver model and resulting observed sky map. This work presents such end-to-end simulations with the intent of determining the impact of 1/f noise on the recovery of the HI angular power spectrum. The paper proceeds as follows: In Section 2, an overview of SD HI IM will be given as well as a brief review of 1/f noise; Section 3 will describe the simulation pipeline that was used to generate the mock SKA HI IM data and the data analysis methods; Section 4 will describe how 1/f noise will impact SKA observations; and Section 5 will derive an empirical model of the impact of 1/f noise and potential avenues of future investigation into the impact of 1/f noise on HI IM. Finally, in Section 6, the key results and insights of the paper are summarized. 2 SINGLE-DISH INTENSITY MAPPING AND 1/f NOISE The SKA array will be a revolutionary interferometric radio observatory. It has been proposed in Santos et al. (2015) and Bull et al. (2015) that the approximately 200 15 m dishes in the SKA–MID array could each be used as independent SD radio telescopes (i.e. autocorrelation or total-power observations). The advantage of SD HI IM over interferometric surveys is increased HI surface brightness and sensitivity to HI on large angular scales. The above papers suggest that an HI IM survey using the SKA could yield cosmological constraints that are competitive with Stage IV dark-energy surveys on comparable time-scales. There are disadvantages to using the SKA as an array of SD instruments. The first is the poor resolution of SD observations. The limited resolution of an SKA HI IM survey means that spatial scales measured in the angular power spectrum will only be useful for late-time cosmology or z< 1 (assuming approximately ∼1 deg resolution). However, at low frequencies, there will still be sensitivity to HI fluctuations in the frequency direction that could be utilized. The limited resolution is especially important if the objective is to measure the characteristic wiggles of baryonic acoustic oscillations (BAOs) within the HI power spectrum on scales of 0.1 Mpc−1. This was discussed extensively in Bull et al. (2015), and is easily summarized by looking at the extent of the beam resolution in k-space as \begin{equation*} k_{{\rm beam}} \approx \frac{2 \pi D_{{\rm dish}}}{ r \lambda }, \end{equation*} (1) where Ddish is the diameter of the telescope dish (accounting for underillumination), λ is the wavelength of observation, and r is the comoving distance at a given redshift. For even the largest SD radio telescopes, such as the 500 m diameter FAST telescope, the maximum redshift at which BAOs could still be measured is z < 4. Another disadvantage of using the SD approach is the challenge of separating spurious contaminants from the HI signal. When observing the sky with an interferometric array, there are several ways to utilize the combined information of all the telescopes to suppress external sources of RFI (e.g. Hellbourg et al. 2014), or instrumental effects such as standing waves and 1/f noise that are specific to individual receivers within the array. However, this is not possible with SD radio telescopes and suppressing systematics or spurious signals requires a combination of careful modelling and a well-designed observing strategy. This work will focus on the contamination of SD observations by 1/f noise, which is a form of correlated noise that is ubiquitous to radio receiver systems and manifests as small gain fluctuations. When binned into a map of the sky 1/f noise manifests as large-scale spatial fluctuations that are not trivial to separate from the true underlying sky signal. Exactly how the fluctuations manifest in sky maps depends on the noise properties (for example the Gaussianity of the noise, the characteristic time-scale of the fluctuations or whether it is stationary) and the details of the observing strategy. For these reasons, 1/f noise has for a long time has been of interest to SD observers in both radio (e.g. Seiffert et al. 2002) and submillimetre (e.g. Emerson & Graeve 1988; Traficante et al. 2011). The removal of 1/f noise has also become an area of research that has resulted in several advanced map-making methods (e.g. Natoli et al. 2001; Sutton et al. 2010). Before discussing 1/f noise, it should be made clear that it is a phenomenon separate to the well understood concept of thermal noise. Thermal noise is caused by the thermal motion of charge carriers and has an rms related to the noise temperature of the receiver through the well-known radiometer equation \begin{equation*} \sigma = A T_{{\rm sys}}\sqrt{ \frac{ f_{{\rm sr}}}{\delta \nu } }, \end{equation*} (2) where A is constant that is dependent on the receiver system of order unity (Wilson, Rohlfs & Hüttemeister 2009), Tsys is the system temperature, fsr is the sample rate, and δν is the system channel width. Thermal noise can be accurately modelled as a Gaussian white noise distribution (Wilson et al. 2009). There have been several proposed physical explanations for 1/f noise within electronic circuits stemming from its initial discovery (Johnson 1925; Nyquist 1928). However, even with its prevalence in a range of physical systems including cognition (Gilden, Thornton & Mallon 1995), biomechanics (Kobayashi & Musha 1982), geological records (Mandelbrot & Wallis 1969), music (Voss & Clarke 1978), and others (e.g. Press 1978), a fundamental description of 1/f noise has yet to be found. Fortunately, 1/f noise can still be modelled and phenomenologically described by a small number of statistical properties. It is common in astronomy to define the power spectral density (PSD) of a receiver contaminated with thermal and 1/f noise as (e.g. Seiffert et al. 2002; Bigot-Sazy et al. 2015) \begin{equation*} {\rm PSD}(f) = \frac{T_{{\rm sys}}^{2}}{\delta \nu } \left[1 + \left(\frac{f_{\mathrm{ k}}}{f}\right)^{\alpha} \right], \end{equation*} (3) where Tsys is the system temperature of the receiver, δν is the channel bandwidth, fk is known as the knee frequency, and α is the spectral index of the noise. The unity term inside the brackets of equation (3) describes the thermal noise contribution and the 1/f noise is described by the reciprocal power law on the right. The reciprocal power law, from which 1/f noise derives its name, is its key property and implies that for α > 0 long time-scale fluctuations will have more power than short time-scale fluctuations. It is common to find named types of 1/f noise that describe specific spectral indices, such as pink noise (α = 1), and brown noise (α = 2), or sometimes generally referred to as red noise for any α > 0. Another property of 1/f noise is that the fluctuations need not be Gaussian distributed, which could impact how it averages down over time. However, in this work, it is assumed that 1/f noise does have Gaussian properties, this should be confirmed through measurements from real receiver systems. A third important property of 1/f noise is that equation (3) only defines it for a finite period, as the 1/f noise term in equation (3) is unbounded and tends to infinity at the zeroth frequency. This of course cannot be true because it would be in violation of both Parseval’s theorem and conservation of energy. One simple solution to this paradox would be to suppose that on some sufficiently long time-scale the 1/f PSD flattens. This must be true, but intriguingly, for semiconductors at least, no such turn-off has been observed even after months of continuous observation (Caloyannides 1974; Mandal, Arfin & Sarpeshkar 2009). For real astronomical observations, this is also not true as no observation is taken over infinite time, and data calibration effectively acts as a high-pass filter on the lowest frequency 1/f noise modes. For many past CMB and other SD experiments, each receiver has a single output that is integrated over a wide bandwidth. Therefore, the 1/f noise fluctuations in each receiver are entirely independent and as such can be sufficiently characterized by the parameters fk and α from equation (3). For the spectroscopic receivers used in IM experiments, each receiver will have many outputs. In this instance, the 1/f noise fluctuations in two output channels from one receiver are likely to be highly correlated. Similarly, the width of each channel is arbitrary with wider channels having lower thermal noise levels. However, if the 1/f noise is highly correlated then it will not average down with wider channel widths and so will effectively increase the fk defined in equation (3). For this reason, the use of fk to characterize 1/f noise properties of a spectroscopic system should be used with care, as it depends on both the channel bandwidth and the correlations of the 1/f noise in frequency. At present, there is very little known about the frequency correlations of 1/f noise in spectroscopic receivers. This work will explore the impact of frequency correlated 1/f noise on a simulated SKA HI IM survey. Due to lack of knowledge as to what the true functional form for the PSD of frequency correlated 1/f noise may be, a simple power-law model is adopted. Modifying equation (3) to include frequency correlations results in \begin{equation*} {\rm PSD}(f, \omega ) = \frac{T_{{\rm sys}}^{2}}{\delta \nu } \left[1 + C(\beta , N_{\nu} ) \left(\frac{f_\mathrm{ k}}{f}\right)^{\alpha} \left( \frac{1}{\omega \Delta \nu } \right)^{\frac{1-\beta }{\beta }} \right], \end{equation*} (4) where ω is the inverse spectroscopic frequency wavenumber, Δν is the total receiver bandwidth, β is used to parametrize the spectral index of the PSD, and C(β, Nν) is a constant that is discussed in more detail in Section 3.3. Equation (4) describes a 2D PSD for which the temporal and spectroscopic 1/f noise fluctuations are separable. The spectral index of the frequency correlations is defined in the limits 0 < β < 1, where β = 0 describes 1/f noise fluctuations that are identical in every frequency channel, whereas β = 1 would describe 1/f noise that is independent in every channel. In Section 3.3, the 1/f noise model and simulations will be discussed in detail. Fig. 1 shows a toy example of equation (4), taking slices of the PSD along the smallest wavenumber mode in time and frequency for an arbitrary receiver system with a β = 0.25 and α = 1. It is expected that 1/f noise will be highly correlated, and as such the β value will be small. Fig. 1 shows that in this case the spectroscopic 1/f noise slope will be very steep, and may only dominate on the very largest spectroscopic scales making measurements of β challenging. Fig. 1 also defines a knee for the spectroscopic 1/f noise fluctuations (ωk) that is not defined in equation (4). This knee is intrinsically linked to the temporal knee frequency fk by \begin{equation*} \omega_{\mathrm{ k}}= \left(f_{\mathrm{ k}} T_{\rm obs}\right)^{\frac{\alpha \beta }{1-\beta }} \Delta \nu ^{-1}, \end{equation*} (5) where Tobs is the total observing time per receiver and Δν is the total bandwidth of the receiver. Figure 1. View largeDownload slide The PSD of the temporal (top) and spectroscopic (bottom) 1/f noise as described by equation (4). Both plots have the same thermal noise level ($$\omega _w^2$$), and an α = 1 and a β = 0.25 for the temporal and spectroscopic PSD, respectively. The black line is the combined spectral power of the 1/f noise slope and the thermal noise. Both the thermal noise and the 1/f noise have equal spectral power at fk and ωk marked by the red dotted line in the top and bottom plots, respectively. Figure 1. View largeDownload slide The PSD of the temporal (top) and spectroscopic (bottom) 1/f noise as described by equation (4). Both plots have the same thermal noise level ($$\omega _w^2$$), and an α = 1 and a β = 0.25 for the temporal and spectroscopic PSD, respectively. The black line is the combined spectral power of the 1/f noise slope and the thermal noise. Both the thermal noise and the 1/f noise have equal spectral power at fk and ωk marked by the red dotted line in the top and bottom plots, respectively. 3 PIPELINE AND SIMULATIONS Mock data sets are a crucial component of any precision cosmology experiment. For systematics in radiometers that are time dependent, the only way to study how these contaminants interact with the signal of interest is through end-to-end simulations. Within this section, a description of an end-to-end simulation pipeline for SD HI IM experiments is given. The pipeline simulates the expected total-power outputs from each radiometer within an arbitrary array. The simulated receiver time-ordered data (TOD) model is described by \begin{equation*} d(t, \nu ) = \left[1 + \delta G(t,\nu )\right]\left[T_{\rm rx} + T_{\rm CMB} + T_{\rm sky}\right], \end{equation*} (6) where d is the TOD vector output from the simulations, Trx + TCMB is the combined receiver and CMB temperatures (see Table 1), Tsky is the sum of the sky model (including the HI signal) as described in Section 3.2, and δG describes the multiplicative 1/f noise contribution discussed in Section 3.3. The pipeline has been designed to be modular such that in the future it can be expanded to include other systematics expected to impact HI IM experiments such as RFI, standing waves, calibration errors, and more. The pipeline is written in a combination of python and fortran, and makes use of mpi functionality though the mpi4py module (Dalcin et al. 2011). 3.1 Instrument design The simulated SKA array described in this section is based on the design outlined in the SKA 2016 baseline document SKA Collaboration et al. (2016), and several recent SKA HI IM forecast and modelling papers (Santos et al. 2015; Bull et al. 2015; Alonso et al. 2015; Villaescusa-Navarro, Alonso & Viel 2017). The simulated survey will use a subset of the Band 2 frequencies spanning the frequency range 950 < ν < 1410 MHz with Nν= 23 frequency channels, each with a channel width of δν= 20 MHz. The choice of channel width was motivated to optimize between the effectiveness of the component separation, which performs better with more channels (Olivari, Remazeilles & Dickinson 2016), and computational efficiency. The corresponding upper redshift limit will be z ≈ 0.5. The choice of Band 2 over Band 1 is based on the limited resolution of the SKA1–MID dishes at redshifts greater than z > 0.5 resulting in an insensitivity to the BAO signal as discussed in Section 2. The number of dishes used is Ndish = 200, which is several more than the current baseline SKA description, but parametrization of the results presented later will allow this number to be scaled to that of the resulting final array. The precise telescope positions do not have a significant impact on the final results. Table 1 provides a description of all the other parameters. Table 1. Input parameters describing the simulated SKA1–MID array and instrumentation. Description Parameter Value Dish diameter Ddish 15 ma Number of dishes Ndish 200 Receiver + CMB TCMB + Trx 20 Kb Number of polarimeters Npol 2 Number of channels Nν 23 Bandwidth Δν 950 < ν < 1410 MHz Channel width δν 20 MHz Sample rate fsr 4 Hz Integration time Tobs 30 d Elevation E 55 deg Slew speed vt 0.5 < vt < 2.0 deg s−1 Description Parameter Value Dish diameter Ddish 15 ma Number of dishes Ndish 200 Receiver + CMB TCMB + Trx 20 Kb Number of polarimeters Npol 2 Number of channels Nν 23 Bandwidth Δν 950 < ν < 1410 MHz Channel width δν 20 MHz Sample rate fsr 4 Hz Integration time Tobs 30 d Elevation E 55 deg Slew speed vt 0.5 < vt < 2.0 deg s−1 aThe slightly smaller dish sizes used in the MeerKAT array (Ddish = 13.5 m) is ignored for simplicity. bThe total system temperature also includes a frequency- and position-dependent contribution from the sky, however the receiver temperature (Trx) and CMB temperature (TCMB = 2.73 K) are assumed to be constant. View Large Table 1. Input parameters describing the simulated SKA1–MID array and instrumentation. Description Parameter Value Dish diameter Ddish 15 ma Number of dishes Ndish 200 Receiver + CMB TCMB + Trx 20 Kb Number of polarimeters Npol 2 Number of channels Nν 23 Bandwidth Δν 950 < ν < 1410 MHz Channel width δν 20 MHz Sample rate fsr 4 Hz Integration time Tobs 30 d Elevation E 55 deg Slew speed vt 0.5 < vt < 2.0 deg s−1 Description Parameter Value Dish diameter Ddish 15 ma Number of dishes Ndish 200 Receiver + CMB TCMB + Trx 20 Kb Number of polarimeters Npol 2 Number of channels Nν 23 Bandwidth Δν 950 < ν < 1410 MHz Channel width δν 20 MHz Sample rate fsr 4 Hz Integration time Tobs 30 d Elevation E 55 deg Slew speed vt 0.5 < vt < 2.0 deg s−1 aThe slightly smaller dish sizes used in the MeerKAT array (Ddish = 13.5 m) is ignored for simplicity. bThe total system temperature also includes a frequency- and position-dependent contribution from the sky, however the receiver temperature (Trx) and CMB temperature (TCMB = 2.73 K) are assumed to be constant. View Large The beam of each telescope is assumed to be Gaussian and the full width at half-maximum (FWHM) of the beam scales withwavelength (λ) as \begin{equation*} \theta _{{\rm FWHM}} = 1.1 \frac{\lambda }{D_{{\rm dish}}}. \end{equation*} (7) The scaling factor of 1.1 comes from measurements of the SKA–MID primary beam models (Lehmensiek, private communication). For these simulations, it is assumed that all observations are taken at the same resolution equal to the resolution at 950 MHz, which corresponds to θFWHM = 1.33 deg. This approximation was made to simplify the noise estimates in data analysis stage of the simulations. The chosen observing strategy is to continuously slew the SKA dishes 360 deg at a constant elevation of 55 deg. Observing at a constant elevation is commonly used by ground-based SD survey instruments such as CBASS (Irfan et al. 2015) and QUIJOTE (Rubiño-Martín et al. 2010). The advantage of scanning at a constant elevation is to control systematics such as the local horizon or ground pick-up. In practice, several elevations would be used to more evenly distribute observing time over the entire observed field, however to simplify the data analysis of the simulations only one elevation is used here. The choice of 55 deg elevation results in approximately a survey area of 20500 sq. deg, which was chosen to be comparable with the survey areas proposed in Santos et al. (2015) and Bull et al. (2015). Simulations were generated for surveys with slews speeds of vt = 0.5, 1, and 2 deg s−1. These speeds were chosen as the SKA–MID dishes are designed to slew at speeds <1 deg s−1, while maintaining sufficient pointing accuracy (SKA Collaboration et al. 2016). The simulations at 2 deg s−1 are provided to determine whether there is sufficient motivation to push to higher slew speeds. The chosen receiver sample rate of 4 Hz is sufficient to Nyquist sample the sky up to 2 deg s−1 at an elevation of 55 deg. 3.2 Sky model 3.2.1 Synchrotron Synchrotron radiation originates from the interaction of relativistic cosmic ray electrons (CRE) interacting within the Galactic magnetic field. The dominant source of CRE are generated in supernovae and have a lifetime of approximately 105–106 yr, this means that synchrotron emission is a tracer for relatively recent star formation within the Galaxy (Condon 1992). At low radio frequencies, synchrotron emission is vastly brighter than any other emission from the sky. The best measurement of the all sky distribution of Galactic synchrotron intensity comes from a reprocessing of the Haslam et al. (1982) 408  MHz all-sky map by Remazeilles et al. (2015). It is well known that Galactic synchrotron emission varies smoothly with frequency and the ensemble average of the emission along a given line of sight can be approximated by a power-law (e.g. Scheuer & Williams 1968; de Oliveira-Costa et al. 2008; Planck Collaboration et al. 2016b). The synchrotron model used in this work is \begin{equation*} T(\nu , \hat{n}) = T_{408{\rm \,MHz}}(\nu , \hat{n}) \left(\frac{\nu }{408{\rm \,MHz}}\right)^{\alpha _s(\nu , \hat{n})}, \end{equation*} (8) where ν is frequency, $$\hat{n}$$ is a line of sight on the sky, T408 MHz is the pixel amplitude of the 408 MHz map, and αs is the spectral index for a given line of sight. The spectral index from pixel-to-pixel was estimated using the all-sky spectral index map by Platania et al. (2003). The spectral indices do not include spectral curvature. The top image of Fig. 2 shows the simulated Galactic synchrotron emission at 1190 MHz for the simulated SKA strip of sky. The bright features to the left and right of the image are saturated slices through the Galactic plane emission, which were masked before performing the component separation analysis. Looking at the total simulated emission map, which combines all foregrounds and the HI signal, it is clear that the synchrotron emission is the dominant foreground by at least an order of magnitude. It is also clear that, even after masking out the brightest low Galactic latitude emission, the residual foreground components are still far brighter than the underlying HI signal. Figure 2. View largeDownload slide Cartesian projections in celestial coordinates of the simulated SKA HI IM survey strip for (from top to bottom): Galactic synchrotron emission (K), Galactic free–free emission (K), the cosmic HI intensity field (mK), the combined sky emission (K), and the integration time distribution per pixel (seconds) for the adopted scanning strategy. All images use a healpix grid with Nside= 512 and centred at right ascension α = 0 and declination δ = −28.5 deg. All colour scales are linear with high brightness regions saturated to highlight the low brightness emission. Due to the projection of spherical data on to a Cartesian plane low declination region will look stretched in the right ascension direction. Figure 2. View largeDownload slide Cartesian projections in celestial coordinates of the simulated SKA HI IM survey strip for (from top to bottom): Galactic synchrotron emission (K), Galactic free–free emission (K), the cosmic HI intensity field (mK), the combined sky emission (K), and the integration time distribution per pixel (seconds) for the adopted scanning strategy. All images use a healpix grid with Nside= 512 and centred at right ascension α = 0 and declination δ = −28.5 deg. All colour scales are linear with high brightness regions saturated to highlight the low brightness emission. Due to the projection of spherical data on to a Cartesian plane low declination region will look stretched in the right ascension direction. 3.2.2 Free-Free Free–free radiation is generated from unbound interactions within regions of the interstellar medium that have been ionized by nearby OB stars. Free–free emission is a very well understood Galactic foreground component and for diffuse, optically thin HII regions have a spectral index of αs = −2.1 at frequencies around 1 GHz (e.g. Draine 2011). For low Galactic latitudes, the radio free–free emission was simulated using estimates of emission measure from the Planck Galactic component maps (Planck Collaboration et al. 2016a). However, for these particular simulations such low Galactic latitudes are masked, and the remaining intermediate-to-high Galactic latitude free–free emission was modelled using the combined all-sky Hα emission map and Hα-to-radio relations described in Dickinson, Davies & Davis (2003). For both models, the mean electron temperature was fixed at 7000 K for all lines of sight (Alves et al. 2012). Although free–free emission is the subdominant foreground as shown in Fig. 2, it will still act to flatten the overall foreground spectrum at higher frequencies. Therefore, free–free emission effectively adds spectral curvature to the foreground components, making the foreground signal more challenging to characterize and remove. 3.2.3 HI The HI emission at low redshifts of z < 0.5 is assumed to be confined to galaxies. The HI signal was simulated using the model for the mean HI brightness described in Battye et al. (2013) \begin{equation*} \bar{T}_{{\rm obs}}(z) = 44\,\mu {\rm K} \left(\frac{\Omega _{{\rm HI}}(z)h }{2.45 \times 10^{-4}} \right) \frac{(1 + z)^2}{E(z)}, \end{equation*} (9) where $$\bar{T}_{\rm obs}$$ is the mean observed brightness temperature of the HI at redshift z, ΩHI is the neutral HI fraction, and E(z) describes the Hubble expansion. It is assumed that there is no evolution in the neutral HI fraction with redshift and a value of ΩHI = 6.2 × 10−4 is adopted (Switzer et al. 2013). The underlying cold dark matter power spectrum, Pcdm, was generated using the camb software package (Lewis & Bridle 2002). To project the matter spectrum into the angular power spectrum at each redshift, the Limber approximation was assumed (Datta, Choudhury & Bharadwaj 2007), \begin{equation*} C_{\ell } = \frac{H_0 b^2}{c} \int {\rm d}z E(z) \left[ \frac{\bar{T}_{{\rm obs}}(z) D(z)}{r(z)} \right]^2 P_{{\rm cdm}}\left(\frac{\ell + \frac{1}{2}}{r} \right) \end{equation*} (10) where r(z) is the comoving distance out to redshift z, D(z) is the growth factor, and b is the HI bias (assumed here to be constant and unity). Each map realization was generated using the synfast module within the healpix package (Górski et al. 2005). Fig. 2 shows the HI emission for the simulated SKA strip at 1190 MHz. The figure demonstrates the foreground subtraction challenge faced by HI IM experiments as the typical fluctuations of the HI signal in figure are of the order ∼10 μK, while the largely Galactic synchrotron foreground fluctuations are several orders of magnitude higher around ∼100 mK. 3.3 Modelling 1/f noise The 1/f noise model used in these simulations assumes that the source of the 1/f fluctuations in the TOD originate as correlated fluctuations in the gain of the receiver amplifiers. The 1/f noise fluctuations are therefore small multiplicative variations around the system temperature. The 1/f noise can be represented in the TOD as \begin{equation*} \Delta T(t, \nu ) = \delta G(t, \nu ) T_{{\rm sys}}(t, \nu ), \end{equation*} (11) where ΔT(t, ν) is the power of the 1/f fluctuations for a given interval, which is the combination of the instantaneous fluctuation in the gain δG and system temperature Tsys. It is the fluctuations in the gain, δG, that this section will describe. The PSD of δG can be described by four parts. The first is \begin{equation*} F(f) = \frac{1}{\delta \nu }\left( \frac{f_\mathrm{ k}}{f} \right)^\alpha , \end{equation*} (12) which is the power spectrum of the temporal fluctuations as in equation (3), but without the white noise component. δν describes the frequency channel bandwidth, fk is the knee frequency, and α is the spectral index of the noise. The second component of the 1/f power spectrum describes the correlations of the noise in frequency. At present, there is no literature on the possible functional form of the 1/f noise frequency–space PSD, and so a conservative power-law model was adopted (similar to the temporal fluctuations) \begin{equation*} H(\omega ) = \left( \frac{\omega _0}{\omega } \right)^{\frac{1-\beta }{\beta }} \end{equation*} (13) where ω is the Fourier mode of the spectral frequency (i.e. the wavenumber), ω0 is the smallest wavenumber (1/Δν), and 0 ≤ β ≤ 1 describes the correlations in frequency. For example, β = 1 describes 1/f noise that is entirely uncorrelated in frequency and β = 0 would result in the perfect correlation of the 1/f noise across channels. This model assumes no physical interpretation for the origin of the 1/f noise frequency correlations. Fig. 3 shows examples of 1/f noise simulated with different values of β with stationary temporal noise properties (e.g. the total rms, fk, and α are all constant with time). The figure shows the Fourier transform of equations (12) and (13) combined with randomized phases. The mean frequency spectrum and the mean temporal fluctuations are both zero. The simulated gain fluctuations should be interpreted as ripples across a 2D surface. Strictly, there is no white noise component in these simulations, however when β = 1, the frequency spectrum is entirely uncorrelated giving an effective appearance of white noise. This means that the number of modes needed to describe the β = 1 1/f noise is equal to the number of channels, therefore removing the 1/f noise will be very challenging for typical component separation methods. Conversely, fewer modes are required for 1/f noise with β < 1, implying that component separation will find it easier to remove the 1/f noise contamination. Figure 3. View largeDownload slide Time–frequency plots of 1/f noise generated using the simulation pipeline. (a) is example of highly correlated 1/f noise, β = 0.25, while (b) β = 0.5, and (c) β = 0.75 are increasingly uncorrelated in frequency, and finally (d) β = 1 is an example of completely uncorrelated 1/f noise where each channel has a unique noise realization. The images represent waterfall plots, with time along the x-axis and frequency along the y-axis. The parallel and vertical plots associated with each image are slices along either a single channel or time interval, respectively. These figures do not include any white noise contribution. Figure 3. View largeDownload slide Time–frequency plots of 1/f noise generated using the simulation pipeline. (a) is example of highly correlated 1/f noise, β = 0.25, while (b) β = 0.5, and (c) β = 0.75 are increasingly uncorrelated in frequency, and finally (d) β = 1 is an example of completely uncorrelated 1/f noise where each channel has a unique noise realization. The images represent waterfall plots, with time along the x-axis and frequency along the y-axis. The parallel and vertical plots associated with each image are slices along either a single channel or time interval, respectively. These figures do not include any white noise contribution. For each example in Fig. 3, the rms of the temporal fluctuations is fixed, however the amplitude of the fluctuations in frequency depend on the β correlation parameter, with the uncorrelated noise having far larger fluctuations than the correlated fluctuations. Over long integrations, the uncorrelated 1/f noise will average as $$\sqrt{1/T}$$ (where T is integration time), whereas the correlated 1/f noise will average down slower depending on the β parameter. This can result in the reversal of what is observed in Fig. 3, as demonstrated by Fig. 4. The figure shows 1/f noise within the simulated SKA field after 30 d of integration, and that the resulting 1/f noise fluctuations in frequency are higher for β = 0.25 than β = 1. For the analysis in this work, the higher amplitude of the correlated fluctuations after integration is not problematic as component separation is very effective at subtracting smoothly varying correlated signals. However, such an effect could be a concern for experiments interested in the cosmological HI signals frequency–space correlations such as redshift space distortions. Further, the assumption that the 1/f noise integrates down as $$\sqrt{1/T}$$ is dependent on both the Gaussianity and stationarity of the 1/f noise, either of which may not be true for real instruments. Figure 4. View largeDownload slide Spatial and radial slices of 1/f noise with (a) β = 0.25 , and (b) β = 1.0 . The radial slice through redshift in the left-hand plots is taken from the declination highlighted by the red stripe in the right-hand spatial image. The spatial images have central frequencies of 960 MHz. The radial slices are presented as arcs to illustrate that these are cuts through a sphere. Due to the different frequency correlation parameters, the 1/f noise averages down far slower in part (a) than in part (b). Figure 4. View largeDownload slide Spatial and radial slices of 1/f noise with (a) β = 0.25 , and (b) β = 1.0 . The radial slice through redshift in the left-hand plots is taken from the declination highlighted by the red stripe in the right-hand spatial image. The spatial images have central frequencies of 960 MHz. The radial slices are presented as arcs to illustrate that these are cuts through a sphere. Due to the different frequency correlation parameters, the 1/f noise averages down far slower in part (a) than in part (b). For small numbers of channels, β and the amplitude of the 1/f noise frequency fluctuations become highly dependent when using equation (13) due to aliasing. For example, in the extreme case of just one channel, the choice of β no longer matters and the rms of the fluctuations must be zero. However, for these simulations, a sufficient number of frequency channels have been used such that the dependence of β and amplitude of the fluctuations on the number of frequency channels is expected to be minimal. A high-pass cosine filter is applied to the PSD to remove correlations on very long time-scales. The filter is used in these simulations to avoid any additional complexity in the interpretation of the results that may be due using an advanced map-making method. The filter has the form \begin{eqnarray*} W(f) \!=\!\left\lbrace \begin{array}{l@{\quad}l}1, \!&\! \text{if}\quad |f| \gt f_\mathrm{ c} + \frac{f_\mathrm{ w}}{2} \\ 0, \!&\! \text{if}\quad |f| \lt f_\mathrm{ c} - \frac{f_\mathrm{ w}}{2} \\ \frac{1}{2}\left[ 1 + \cos \left(\pi \frac{|f| - f_\mathrm{ c} - f_\mathrm{ w}/2}{f_\mathrm{ w}}\right) \right], \!&\! \text{otherwise} \end{array}\right. \end{eqnarray*} (14) where fc is the temporal frequency for which half the signal power is filtered, and fw is the width of the filter. The choice of the filter cut-off scale was chosen to synchronize with the period of a single 360 deg slew in azimuth as this is both easy to achieve in reality and suppresses spikes in the power spectrum caused by the interaction of the 1/f noise and the scan strategy at the pixel scale. Though scan speed is one variable investigated in this paper, the actual details of the scan strategy are not. The overall power of the 1/f noise angular power spectrum is only weakly dependent on filter scale used as long as the filter period is greater than the slew time. The final component needed to characterize the 1/f noise PSD is a normalization factor that forces the temporal variance of δG to be constant regardless of the β chosen for the correlations in frequency. The correction factor is \begin{equation*} C(\beta , N_{\nu }) = \frac{N_{\nu }}{1 + \frac{\sum \limits _{i=1}^{N_\nu /2}H(\omega _i, \beta )\Delta \omega }{\int H (\omega , 0)\text{d}\omega } \left(N_\nu - 1 \right) }, \end{equation*} (15) where Nν is the number of channels. For a small number of channels, it is important to use a discrete sum and not the continuous integral over equation (13). The resulting 1/f PSD can now be described by \begin{equation*} {\rm PSD}(f,\omega ) = F(f) H(\omega ) W(f) C(\beta , N_\nu ), \end{equation*} (16) where F(f) is from equation (12), H(ω) is equation (13), W(f) is equation (14), and C(β, Nν) is equation (15). Each 1/f noise TOD realization is generated using Gaussian distributed phases with means equal to the square root of the power as \begin{equation*} Z(f, \omega ) = \sqrt{{\rm PSD}(f,\omega )}\left( x + iy \right) \end{equation*} (17) where x and y are Gaussian random variates with zero mean and unit variance. The Fourier transform of Z results in the 1/f gain fluctuations \begin{equation*} \delta G (t, \nu ) = \int \int \text{d}f \text{d}\omega \, Z(f, \omega ) \mathrm{ e}^{2\pi i (t f + \nu \omega ) }. \end{equation*} (18) It is an assumption of these simulations that the 1/f noise is Gaussian distributed. In real data that may not be true. Also, as the 1/f noise is multiplicative with system temperature, which can be time variable, it is non-stationary in terms of its variance (e.g. the rms of the 1/f noise increases when observing bright sources on the sky), but the fk, α, and β of the 1/f noise are all assumed stationary, which also may not be true in real observations. Alone or in combination, these effects will influence how the 1/f noise averages down on long time-scales. 3.4 Sky mask Before performing the component separation and power spectrum estimation discussed in the proceeding subsections, the output of the sky maps from the simulation pipeline were multiplied with a mask. The mask is based on a threshold cut of 40 K in the Haslam 408 MHz sky map (Haslam et al. 1982; Remazeilles et al. 2015), which was chosen to agree with other SKA simulations by Alonso et al. (2015). It removes just the brightest regions of the Galactic plane. The mask also removes edge pixels that, due to the arrangement of the SKA dishes, have very few hits in the highest and lowest observed declinations relative to the rest of the observed sky. Fig. 2 shows the regions that have been masked, the resulting total observed sky fraction after masking is fsky = 0.3. 3.5 Component separation There are currently a wide range of parametric and non-parametric methods for component separation that have been developed during the past decade for CMB experiments, but are now being applied to HI IM. Examples include Fast Independent Component Analysis (FASTICA) (Chapman et al. 2012), correlated component analysis (Bonaldi & Brown 2015), Generalised Morphological Component Analysis (GMCA) (Chapman et al. 2012), and Generalised Needlet Internal Linear Combination (GNILC) (Olivari et al. 2016). For the purposes of these simulations, the method of principal component analysis (PCA) was chosen. PCA is not as sophisticated a method as those mentioned above, and as such is not as robust at removing foregrounds. However, PCA has the advantage of being computationally fast, making it highly suited to Monte Carlo simulations. Also, PCA relies upon eigenvalue decomposition, which is at the heart of most non-parametric component separation methods. The method of PCA works by assuming that a map that is composed of Npix  pixels and Nν frequencies can be decomposed into a set of Ne independent components (where Ne ≤ Nν) via a mixing matrix as \begin{equation*} {\bf m} = {\bf U} {\bf s}, \end{equation*} (19) where m is the set of Nν maps, U is the mixing matrix, and s is the set of Ne independent components. The formalism adopted here is that vectors are lower case, bold characters, matrices are upper case, bold characters, and rank-0 variates are lower case and not bold. The covariance matrix of the maps is defined as \begin{equation*} {\bf m} {\bf m}^{{\rm T}} = {\bf C} = {\bf U} {\bf s} {\bf s}^{{\rm T}} {\bf U}^{{\rm T}} = {\bf U} \Lambda {\bf U}^{{\rm T}}, \end{equation*} (20) where C is the covariance matrix, and it has been assumed that since each component of s is independent, the outer product yields a diagonal matrix Λ. To determine the mixing matrix U and the diagonal matrix Λ, an eigenvalue decomposition can be performed on C. This means that diagonal elements of Λ are the eigenvalues of C and the columns of the mixing matrix U are the corresponding eigenvectors. The mixing matrix U also obeys the condition that \begin{equation*} {\bf U}^{{\rm T}} {\bf U} = {\bf U} {\bf U}^{{\rm T}} = {\bf I}, \end{equation*} (21) where I is the identity matrix. The eigenvalues of C describe the relative variance of each of the components in m. As the Galactic foregrounds are bright relative to the HI signal, and highly correlated it is assumed that they can be described by just a few bright eigenvalues (e.g. Bigot-Sazy et al. 2015). Therefore, to remove the foregrounds, a subspace of the mixing matrix $$\tilde{{\bf U}}$$ is derived by removing the eigenvectors associated with the largest eigenvalues. The foreground subtracted map is then defined as \begin{equation*} \tilde{{\bf U}} {\bf U}^{{\rm T}} {\bf m} = \tilde{{\bf m}} \end{equation*} (22) where $$\tilde{{\bf m}}$$ is the estimated foreground subtracted sky. After applying the PCA method to the simulated data sets with no 1/f noise, it was found that the optimal number of eigenmodes to subtract were Nλ= 6. This minimizes the recovered variance of the HI power spectrum due to foregrounds, but at the cost of the two lowest frequency channels being over subtracted. For the rest of these simulations, these two channels are neglected. Non-parametric methods of foreground subtraction are limited by the available number of frequency channels, Nν, and cannot remove a foreground component that can be decomposed into Nλ > Nν. Similarly, if the sky has large variations in the frequency response of the foregrounds then PCA may work better along some lines of sight more than others. One very important property of PCA is that it will always remove components that are completely correlated in frequency. This means that 1/f noise with β = 0 will always be perfectly removed. This has been verified and therefore is not included in the results as it is equivalent to the white noise only case. A simple proof of this assertion is presented in Appendix A. 3.6 Power spectrum estimation For each realization, frequency and eigenmode the angular power spectrum was measured. Each map was decomposed into spherical harmonic coefficients (alm) (Peebles 1973; Wandelt, Hivon & Górski 2001), \begin{equation*} a_{\ell m} = \int \limits _{\Omega } m(\Omega ) Y_{\ell m}(\Omega ) {\rm d}\Omega , \end{equation*} (23) where m is the observed sky map, Yℓm are spherical harmonics, and Ω describes a given line of sight. The angular power spectrum is then \begin{equation*} C_\ell ^\text{obs} = \frac{1}{2 \ell + 1}\sum \limits _m | a_{\ell m}|^2. \end{equation*} (24) Here, $$C^\text{obs}_\ell$$ describes the observed HI angular power spectrum after foreground subtraction. The observed HI angular power spectrum, in the presence of Gaussian white noise only, is related to the underlying HI spectrum realization by \begin{equation*} C_\ell = \left\langle C_\ell ^\text{obs} - N_\ell \right\rangle = \left\langle \hat{C}_\ell \right\rangle , \end{equation*} (25) where the brackets indicate an average over many noise realizations of the white noise power spectrum, and $$\hat{C}_\ell$$ is the estimate of the Cℓ. Note that Cℓ in this notation is not the underlying theoretical HI angular power spectrum, these are related by \begin{equation*} \left\langle C_\ell \right\rangle = C_\ell ^\text{theory}B_\ell ^2, \end{equation*} (26) where the average now is over many HI sky realizations and Bℓ is a beam window function. The spherical harmonic transform and power spectrum estimation of the simulated maps were performed using the publicly available polspice software (Chon et al. 2004). As only a partial fraction of the sky was observed, neighbouring Cℓ components will be correlated. The correlations can be damped by binning components in ℓ with an approximate analytical bin width of (Hauser & Peebles 1973) \begin{equation*} \Delta \ell \approx \frac{\pi }{\Delta \theta _{0}}, \end{equation*} (27) where Δθ0 is the span of the polar angle over the observed sky area. For these simulations, a binning width of Δℓ = 5 was chosen. Masking of the sky can lead to aliasing of power from one scale to another, which can result in non-physical oscillations within the angular power spectrum. The aliasing was minimized by optimising the apodization and maximum angular scales of the masked sky correlation functions. An apodization of 110 deg, and a maximum angular size of 110 deg were found to minimize aliasing in these simulations. Before power spectrum estimation, each map is smoothed to a common FWHM resolution of 1.3 deg (corresponding to the resolution in the 950 MHz channel). As smoothing is performed after masking, power will be lost from pixels near the mask boundaries, which will result in a biased recovery of the HI power spectrum. To avoid this, polspicewas provided an apodized weight map. The weights were derived by calculating the distance of each unmasked pixel from the nearest masked region using the healpix routine process_mask, and passing these distances through a cosine filter function with a characteristic width of 2 deg. Estimated uncertainties in the angular power spectrum were measured from the standard deviation of all realizations in a given Δℓ bin using \begin{equation*} \Delta F_\ell = \sqrt{ \left\langle \hat{C}_\ell ^2 \right\rangle - \left\langle \hat{C}_\ell \right\rangle ^2}, \end{equation*} (28) where ΔFℓ here denotes the uncertainty in the power spectrum due to a combination of white noise, 1/f noise, residual foregrounds, and any other contaminant within the power spectrum. Note that in most cases discussed later only the additional uncertainty over the expected white noise is of interest, therefore an estimate of the white noise uncertainty is required. The mean angular power spectrum for white noise that is uniformly distributed is well known as (Knox 1995) \begin{equation*} N_\ell = \sigma _\text{pix}^2 \Omega _\text{pix}, \end{equation*} (29) where σpix is the standard deviation of pixels in the map domain, and Ωpix is the solid angle of each pixel. However, for specific realizations, the mean white noise angular power will be different from equation (29) and a transfer function that describes the non-uniform noise distribution and the influence of data processing (e.g. in this case foreground subtraction from Section 3.5) should also be included. An alternative method therefore, which was adopted here, was to simply measure the mean white noise power from simulations of white noise plus signal and foregrounds only, post-foreground subtraction. To calculate the underlying white noise power spectrum uncertainty, it was assumed that the underlying sky signal is non-varying. This means that the additional statistical uncertainty due the white noise is not \begin{equation*} \Delta N_\ell = \sqrt{\frac{2}{(2 \ell + 1)}} \left( C_\ell + N_\ell \right), \end{equation*} (30) as derived in Knox (1995), but instead is given by \begin{equation*} \Delta N_\ell = \sqrt{\frac{2}{(2 \ell + 1)}} N_\ell \sqrt{1 + 2 \frac{C_\ell }{N_\ell }}, \end{equation*} (31) which is derived in Appendix B. 3.7 Summary To conclude this section, the following is a complete summary of each step taken in the process of the simulations. Generate sky foreground and signal maps from models described in Sections 3.2.1 –3.2.3. Simulate TODs for the observing strategy using the model in Section 3.1. For each arbitrary block of TOD produce a correlated 1/f noise signal following the method outlined in Section 3.3. Sample from model sky using observing strategy to produce a noise-free TOD. The noise-free TOD is then multiplied by the 1/f noise gain fluctuations and white noise is added. The TOD is then binned into a single realization of the sky. Separate each sky map realization into eigenmode components using PCA as described in Section 3.5. Calculate the angular power spectrum of every map for each frequency and eigenmode. The parameter space explored included four values of β, four knee frequencies (fk), three telescope slewing speeds (vt), and two 1/f spectral indices (α). For each combination of 1/f noise parameters, 100 realizations were generated for a total of 9600. 4 RESULTS The previous section describes the key components of the HI IM simulation pipeline. The simulation pipeline was used to generate several thousand mock 30 d SKA HI IM surveys with unique 1/f noise properties, but identical astrophysical foreground and HI signals. Using these mocks, HI angular power spectra were recovered. The results presented in this section are predominantly concerned with the impact of 1/f noise correlations in the frequency direction on the statistical and systematic errors of the recovered HI power spectra. Section 4.1 will simply describe the range of 1/f noise parameters that are explored. Section 4.2 will describe the properties of the mean 1/f noise angular power spectra (e.g. how it varies in frequency and harmonic spaces). Then , Sections 4.3–4.5 will describe the impact of 1/f noise on the recovery of the HI angular power spectrum in terms of first the increased statistical uncertainty, second through a bias term, and finally, the combined impact of these two terms, respectively. 4.1 Parameter space The discussions of 1/f noise in Section 3.3 imply that there are four key parameters that must be considered: the spectral index of the temporal 1/f fluctuations (α), the proxy for the spectral index of the 1/f frequency fluctuations (β), the telescope scanning speed (vt), and the 1/f knee frequency (fk) at 20 MHz channel width. The full range of parameters explored by these simulations are listed in Table 2. For each input parameter, there are 100 Monte Carlo realizations. Table 2. The full range of the parameter space explored by this work. Each column represents the input parameters used to generate a given set of 100 Monte Carlo simulations of the 30 d SKA HI IM survey. For each set of realizations with a given β, only one fk was simulated, and other knee frequencies were generated by scaling the noise appropriately before component separation or power spectrum estimation. α 1 vt (deg s−1) 0.5 1 2 β 0.25 0.5 0.75 1 0.25 0.5 0.75 1 0.25 0.5 0.75 1 fk (Hz) 10, 1, 0.1, 0.01 10, 1, 0.1, 0.01 10, 1, 0.1, 0.01 α 2 vt (deg s−1) 0.5 1 2 β 0.25 0.5 0.75 1 0.25 0.5 0.75 1 0.25 0.5 0.75 1 fk (Hz) 10, 1, 0.1, 0.01 10, 1, 0.1, 0.01 10, 1, 0.1, 0.01 α 1 vt (deg s−1) 0.5 1 2 β 0.25 0.5 0.75 1 0.25 0.5 0.75 1 0.25 0.5 0.75 1 fk (Hz) 10, 1, 0.1, 0.01 10, 1, 0.1, 0.01 10, 1, 0.1, 0.01 α 2 vt (deg s−1) 0.5 1 2 β 0.25 0.5 0.75 1 0.25 0.5 0.75 1 0.25 0.5 0.75 1 fk (Hz) 10, 1, 0.1, 0.01 10, 1, 0.1, 0.01 10, 1, 0.1, 0.01 View Large Table 2. The full range of the parameter space explored by this work. Each column represents the input parameters used to generate a given set of 100 Monte Carlo simulations of the 30 d SKA HI IM survey. For each set of realizations with a given β, only one fk was simulated, and other knee frequencies were generated by scaling the noise appropriately before component separation or power spectrum estimation. α 1 vt (deg s−1) 0.5 1 2 β 0.25 0.5 0.75 1 0.25 0.5 0.75 1 0.25 0.5 0.75 1 fk (Hz) 10, 1, 0.1, 0.01 10, 1, 0.1, 0.01 10, 1, 0.1, 0.01 α 2 vt (deg s−1) 0.5 1 2 β 0.25 0.5 0.75 1 0.25 0.5 0.75 1 0.25 0.5 0.75 1 fk (Hz) 10, 1, 0.1, 0.01 10, 1, 0.1, 0.01 10, 1, 0.1, 0.01 α 1 vt (deg s−1) 0.5 1 2 β 0.25 0.5 0.75 1 0.25 0.5 0.75 1 0.25 0.5 0.75 1 fk (Hz) 10, 1, 0.1, 0.01 10, 1, 0.1, 0.01 10, 1, 0.1, 0.01 α 2 vt (deg s−1) 0.5 1 2 β 0.25 0.5 0.75 1 0.25 0.5 0.75 1 0.25 0.5 0.75 1 fk (Hz) 10, 1, 0.1, 0.01 10, 1, 0.1, 0.01 10, 1, 0.1, 0.01 View Large For most of the discussions in this section, a baseline reference simulation will be used with fixed parameters of α = 1 and vt = 1deg s−1. 4.2 Simulated 1/f power spectra To begin this section will examine the properties of the simulated 1/f noise outputs. Fig. 5 shows a comparison between the angular power spectra of the foregrounds, input HI, and 1/f noise for different knee frequencies with a slew speed of 1 deg s−1 and β = 1. For knee frequencies less than 10 mHz, the thermal noise power is greater than the 1/f noise power at all temporal scales for a channel width of 20 MHz. When the knee frequency approaches 1 Hz, then the amplitude of the 1/f noise at low-ℓ is comparable to the amplitude of the underlying HI signal. It is critical to be aware of the knee frequency at which the 1/f noise exceeds the HI signal because it will be shown in the results later in this section that the impact of 1/f noise on the recovered HI angular power spectrum is most significant when the 1/f noise power is comparable to or exceeds the HI angular power. Figure 5. View largeDownload slide Angular power spectra of the foregrounds signal (blue), HI (orange), thermal noise (red), and 1/f noise with different knee frequencies (black lines). These power spectra are for 20 MHz wide channels in the centre of SKA band 2 at 1190 MHz. Figure 5. View largeDownload slide Angular power spectra of the foregrounds signal (blue), HI (orange), thermal noise (red), and 1/f noise with different knee frequencies (black lines). These power spectra are for 20 MHz wide channels in the centre of SKA band 2 at 1190 MHz. Fig. 6 shows the same information as Fig. 5, but in 2D for just 1 Hz knee frequency 1/f noise assuming the baseline model. As the 1/f noise is coupled to the foreground brightness, the 1/f noise angular power spectrum amplitude is slightly larger at low frequencies. Conversely, the amplitude of the HI variations decrease at low frequency. This results in the 1/f noise having significantly more impact on high-redshift observations than low-redshift observations. The impact of this on the low-redshift simulations discussed here is small, but it could potentially become problematic for high-redshift EoR HI IM surveys. Similarly, the 1/f noise impacts large scales more than small scales, leading to a decreased signal-to-noise ratio (S/N) at ℓ < 20. Figure 6. View largeDownload slide Distribution of 1/f noise angular power spectra with frequency and ℓ for the baseline model with fk = 1 Hz. The colour scale shows the 1/f noise angular power, which is largely constant with frequency and a power law in ℓ. The contours show the equivalent input HI angular power spectrum. Both colour scale and contours are in units of μK2. Figure 6. View largeDownload slide Distribution of 1/f noise angular power spectra with frequency and ℓ for the baseline model with fk = 1 Hz. The colour scale shows the 1/f noise angular power, which is largely constant with frequency and a power law in ℓ. The contours show the equivalent input HI angular power spectrum. Both colour scale and contours are in units of μK2. Another informative way to view the 1/f noise is through the reprojection of the spherical harmonic coefficients into a 2D image. Fig. 7 shows the standard deviation of the aℓm’s for 100 different 1/f noise and input HI realiZations. The mapping between the x –y Cartesian grid and the ℓ − m plane is x = ℓ − |m| and y = ℓ, it is assumed that there is a m = −m symmetry and so only the bottom left triangle is shown. The images in Fig. 7 reveal that the combination of the correlated 1/f noise with the observing strategy preferentially gives more power to some aℓm’s over others and has a structure distinct from the underlying HI signal. The difference in the noise and the HI imply that there is a possibility of performing component separation in the harmonic space of the data, which may be complementary to performing component separation in image space. Though not explored further by this work, such differences in the harmonic structure does lead to the possibility of novel power spectrum analysis methods such as m-mode analysis (Shaw et al. 2014, 2015; Berger et al. 2017). Figure 7. View largeDownload slide Standard deviation of spherical harmonic coefficients projected into a 2D image plane for the simulated 1/f noise at 0.5 deg s−1 (left), and the input HI signal (right). The cut-off in HI power at high ℓ is due to the beam of the telescope. For 1/f noise, the beam is not a factor to consider and therefore it still has a contribution at the smallest scales. Figure 7. View largeDownload slide Standard deviation of spherical harmonic coefficients projected into a 2D image plane for the simulated 1/f noise at 0.5 deg s−1 (left), and the input HI signal (right). The cut-off in HI power at high ℓ is due to the beam of the telescope. For 1/f noise, the beam is not a factor to consider and therefore it still has a contribution at the smallest scales. In the above cases, and for the rest of the results presented in this paper, the simulation survey is taken to last 30 d. However, it would be interesting to know whether the results presented here can be simply scaled to different total observing times, number of telescopes or survey sky area. To determine this, SKA simulations for periods of  2–512 d were generated that include only 1/f noise with α = 1 and 2. As may be expected, as the 1/f noise is simulated with Gaussian random variates and stationary values of α and fk, the rms of the maps decrease as $$1/\sqrt{T_\text{obs}}$$. Exploration of how the rms of non-Gaussian and non-stationary 1/f noise is left to futurework. 4.3 Power spectra uncertainty The first question to ask of the 1/f noise in an SKA HI IM survey is how much additional statistical uncertainty does it contribute to the recovered angular power spectrum. One method of determining this is to measure the S/N of the input signal with the combined uncertainties due to noise and cosmic variance. The uncertainty in the HI angular power spectrum is defined as \begin{equation*} \Delta C_\ell = \sqrt{\frac{2}{ \left(2 \ell + 1 \right) \Delta \ell }} C_\ell + \Delta F_\ell \end{equation*} (32) where the first term on the right-hand side is the cosmic variance, calculated analytically from the input HI power spectrum, and the second term ΔFℓ is the uncertainty due to both the thermal and 1/f noise contributions defined in equation (28). The S/N ratio is then defined as \begin{equation*} {\rm S/N}= \frac{C_\ell }{\Delta C_\ell }, \end{equation*} (33) where, as before, Cℓ denotes the angular power spectrum for the input realization of the HI sky signal. Fig. 8 shows the S/N for the frequency range 950 <ν < 1410 MHz for the baseline model of vt = 1 deg s−1, and α = 1. The plot reveals that when the 1/f noise is highly correlated (β = 0.25) in frequency then the contribution to the final HI statistical uncertainty is very small, irrespective of the input 1/f noise knee frequency as expected when using PCA for foreground subtraction and shown in Appendix A. The figure also shows that the survey will be more sensitive to low redshift HI emission, and to scales defined by the ℓ range 10 < ℓ < 100. The effective redshift of the survey can be calculated as \begin{equation*} \bar{z} = \frac{\sum {\rm S/N}^2 z }{\sum {\rm S/N}^2} \end{equation*} (34) where $$\bar{z}$$ is the mean redshift, and z is the redshift of each pixel in the Cℓ − ν plane. The mean redshift is found to be $$\bar{z} \approx 0.2$$. Figure 8. View largeDownload slide The expected S/N when considering only statistical fluctuations in the angular power spectrum due to cosmic, 1/f, and white noise variances as described by equation (33). All plots are using the baseline simulation of α = 1, and vt = 1 deg s−1. Galactic foregrounds have been removed using PCA. The statistical fluctuations in these plots are around 10 per cent as expected for 100 Monte Carlo realizations. Figure 8. View largeDownload slide The expected S/N when considering only statistical fluctuations in the angular power spectrum due to cosmic, 1/f, and white noise variances as described by equation (33). All plots are using the baseline simulation of α = 1, and vt = 1 deg s−1. Galactic foregrounds have been removed using PCA. The statistical fluctuations in these plots are around 10 per cent as expected for 100 Monte Carlo realizations. Summing over each S/N contribution in Fig. 8 results in a combined total S/N \begin{equation*} {\rm S/N}_\text{T} = \sqrt{\sum \limits _\nu \sum \limits _\ell {\rm S/N}^2}, \end{equation*} (35) shown as contour plots in Fig. 9. This figure can be used to numerically determine the impact of the 1/f noise on the HI angular power spectrum recovery. The contour plots show the impact of the 1/f noise when β = 0.25 is comparable to the white noise only when α = 1, with almost no detection when β = 1 and fk ≥ 1 Hz. However, when α = 2, the impact of the 1/f noise is significant even when the noise is highly correlated. This stresses the importance of measuring α of all future SKA receivers that will be used for an HI IM survey. Further the plot shows that small increases in β can very quickly decrease the total S/N of the observation. The impact of α and β on the S/N should be very carefully considered as the observing time is proportional to the square of the S/N, and very quickly the increase in observing time required to achieve the desire sensitivity may become untenable. Figure 9. View largeDownload slide The expected total S/N when considering only statistical fluctuations in the angular power spectrum, after foreground cleaning, due to cosmic, 1/f, and white noise variances as defined by equations (33) and (35). Each value in these plots comes from summing over all the S/N grids, as shown in Fig. 8, for each set of noise parameter inputs. The top row is for α = 1, and from left to right increasing slew speeds of vt = 0.5, 1, and 2 deg s−1. The bottom row is same again, but with 1/f noise with α = 2. The contours show the importance of 1/f noise correlations are when considering the statistical impact of 1/f noise on S/N. Figure 9. View largeDownload slide The expected total S/N when considering only statistical fluctuations in the angular power spectrum, after foreground cleaning, due to cosmic, 1/f, and white noise variances as defined by equations (33) and (35). Each value in these plots comes from summing over all the S/N grids, as shown in Fig. 8, for each set of noise parameter inputs. The top row is for α = 1, and from left to right increasing slew speeds of vt = 0.5, 1, and 2 deg s−1. The bottom row is same again, but with 1/f noise with α = 2. The contours show the importance of 1/f noise correlations are when considering the statistical impact of 1/f noise on S/N. Fig. 10 shows the uncertainty in 1/f noise angular power spectra for α = 1–2, and vt= 0.5, 1, and 2 deg s−1 for fk = 1 Hz and β = 1. The figure shows how the 1/f noise α is a far larger contributor to the overall uncertainty than the slew speed of the telescope. This implies that stable receivers are an order of magnitude more important than choice of scan speed. However, the figure also shows how scan speed is coupled with the spectral index of the 1/f noise as the scan speed is seen to have greater impact for steep spectrum 1/f noise. This coupling between the scan speed and the spectral index is complex but will be strongly dependent on the choice of observing strategy. Figure 10. View largeDownload slide Uncertainty in angular power spectrum averaged over all simulated frequencies for 1/f noise with telescope slew speeds of 0.5, 1, and 2 deg s−1 and spectral indices of α = 1 and 2. All curves are for a knee frequency of 1 Hz and the black-dashed line shows the expected white noise level. There are two principal divisions of the results shown in this figure. The first division is between the solid and dot–dashed lines representing 1/f noise with α = 1 and 2, respectively. The difference in uncertainty within these two divisions is due to the three different scan speeds used, with 2 deg s−1 having the lowest uncertainty in both groups. The key point to take is that the spectral index of the 1/f noise is more important than telescope scan speed. Figure 10. View largeDownload slide Uncertainty in angular power spectrum averaged over all simulated frequencies for 1/f noise with telescope slew speeds of 0.5, 1, and 2 deg s−1 and spectral indices of α = 1 and 2. All curves are for a knee frequency of 1 Hz and the black-dashed line shows the expected white noise level. There are two principal divisions of the results shown in this figure. The first division is between the solid and dot–dashed lines representing 1/f noise with α = 1 and 2, respectively. The difference in uncertainty within these two divisions is due to the three different scan speeds used, with 2 deg s−1 having the lowest uncertainty in both groups. The key point to take is that the spectral index of the 1/f noise is more important than telescope scan speed. As well as the S/N, it is also interesting to know how the ratio of the white noise to the 1/f noise is distributed in ℓ and frequency. The ratio of the 1/f noise to the thermal noise was defined as \begin{equation*} r = \frac{\Delta F_\ell }{\Delta N_\ell }, \end{equation*} (36) where ΔFℓ is the total uncertainty measured from the simulations, and ΔNℓ is the expected uncertainty due to the thermal noise (equation 31). Fig. 11 shows how the ratio r varies with ℓ and frequency for different knee frequencies and β values for the baseline model. The figure shows how 1/f noise can contribute additional uncertainty on all scales, which is dependent only on the ratio of the HI and 1/f noise angular power spectra. At fk = 1 Hz, the residual 1/f noise power exceeds the recovered HI power at all scales and hence there is also additional uncertainty at all scales. The ripples that can be seen in the β = 0.25 column are residual foregrounds from the component separationstep. Figure 11. View largeDownload slide The variation in ratio between the uncertainty in the angular power spectrum, after foreground cleaning, from 1/f noise over white noise calculated using equation (36) a function of ℓ and z. All plots are using the baseline simulation of α = 1, and vt = 1 deg s−1. The statistical fluctuations in these plots are at the 10 per cent level expected for 100 Monte Carlo realizations. Figure 11. View largeDownload slide The variation in ratio between the uncertainty in the angular power spectrum, after foreground cleaning, from 1/f noise over white noise calculated using equation (36) a function of ℓ and z. All plots are using the baseline simulation of α = 1, and vt = 1 deg s−1. The statistical fluctuations in these plots are at the 10 per cent level expected for 100 Monte Carlo realizations. 4.4 Power spectra bias As well as the noise uncertainty, it is also important to quantify the magnitude of the systematics that are introduced by the 1/f noise. The method for characterizing the bias (not to be confused with cosmological bias) is to measure the residual of the angular power spectrum between the input and recovered HI signals defined as \begin{equation*} \epsilon _\ell = \frac{\left\langle \hat{C}_\ell \right\rangle - \left\langle \hat{C}_\ell (f_\mathrm{ k} = {\rm 0\,Hz})\right\rangle }{C_\ell }, \end{equation*} (37) where $$\left\langle \hat{C}_\ell \right\rangle$$ is average over the recovered HI realizations, and $$\left\langle \hat{C}_\ell (f_\mathrm{ k} = {\rm 0\,Hz})\right\rangle$$ is the average recovered HI power spectrum assuming white noise only. Here, realizations of the white noise are used, instead of the analytically derived white noise level, to account for any bias originating from the component separationstep. In Fig. 12, εℓ is presented as a function of frequency and ℓ for the baseline model and a range of fk and β values. The interpretation of the bias in the power spectra is that anything greater than unity implies that the HI power spectrum is dominated by the residual 1/f noise. The bias is only seen to be significant when the 1/f noise power exceeds the HI angular power (i.e. when fk ≥ 1 Hz), however some small-scale contributions to the noise can be seen when fk = 0.1 Hz. Encouragingly, when the 1/f noise is highly correlated (β = 0.25), the bias is close to zero until the 1/f noise greatly exceeds the HI noise power (e.g. fk > 10 Hz). Also note that, similarly to Fig. 11, the ripples that can be seen in the β = 0.25 column are due to poor foreground subtraction. In this case however these ripples are more interesting as, from equation (37), any intrinsic bias due to residual Galactic foregrounds after the component separation step should be subtracted. Therefore, these ripples must be due to residual 1/f noise. Figure 12. View largeDownload slide The relative difference between the mean estimated HI and the input HI angular power spectra, normalized by the input HI spectrum as described by equation (37). All plots are using the baseline simulation of α = 1, and vt = 1 deg s−1. The statistical fluctuations seen at all scales are consistent with 10 per cent, as expected for 100 Monte Carlo with realizations. Figure 12. View largeDownload slide The relative difference between the mean estimated HI and the input HI angular power spectra, normalized by the input HI spectrum as described by equation (37). All plots are using the baseline simulation of α = 1, and vt = 1 deg s−1. The statistical fluctuations seen at all scales are consistent with 10 per cent, as expected for 100 Monte Carlo with realizations. One possible question to ask is whether it would be possible to model this 1/f noise bias and remove it. For these simulations that would be trivial as the input 1/f noise is known, however in real data that may exhibit 1/f noise that is non-Gaussian, non-stationary and coupled with other systematics, accurate modelling may prove challenging. However, should the bias in the real data prove to be a significant difficulty one alternative possibility is to produce cross-correlation HI power spectra from subsets of dishes within the SKA array, but this will come at a cost of increasing the power spectrum uncertainty by at least a factorof 4. 4.5 Combined uncertainties An estimate of the overall uncertainty in the recovered HI $$\hat{C}_\ell$$ is \begin{equation*} \Delta C_\ell = \sqrt{\frac{2}{(2 \ell + 1) \Delta \ell }} C_\ell + \Delta F_\ell + \left\langle \epsilon _\ell \right\rangle C_\ell , \end{equation*} (38) where ΔFℓ is the total uncertainty as measured in the simulations discussed in Section 3.6, and ⟨εℓ⟩ is the mean bias in the recovered power spectra as discussed in Section 4.4. The first term in equation (38) is the cosmic variance term, which has been added to account for the lack of cosmic variance in ΔFℓ. The bias (εℓ) is not a true uncertainty however it is added here to give a tentative value to the impact 1/f noise bias may have on a realobservation. When the bias is added into the uncertainty Fig. 13 shows how the S/N changes. The figure shows that there is a dramatic decrease in the S/N for uncorrelated 1/f noise with β ≥ 0.5. However, even at 10 Hz when β = 0.25 some bins do exhibit a >3σ detection on the HI angular power spectrum, which likely could be improved using a more advanced component separation method. Figure 13. View largeDownload slide The S/N when including not just statistical variations from equation (33) but also includes the systematic residual power from equation (37), as described by equation (38). The overall distribution of the highest S/N regions in this plot are dominated by the residuals induced by poor foreground subtraction and masking. The statistical variations in this plot are consistent with the 10 per cent level expected for 100 Monte Carlo noise realizations. Figure 13. View largeDownload slide The S/N when including not just statistical variations from equation (33) but also includes the systematic residual power from equation (37), as described by equation (38). The overall distribution of the highest S/N regions in this plot are dominated by the residuals induced by poor foreground subtraction and masking. The statistical variations in this plot are consistent with the 10 per cent level expected for 100 Monte Carlo noise realizations. The six plots in Fig. 14 shows the total S/N ratio for the entire parameter space. Comparing this figure to Fig. 9, it can be seen that the residual 1/f noise bias has significant impact for high values of β, decreasing the S/N by factors of 2–5 depending on α and vt. As has previously been discussed, this plot further emphasizes the importance of α on the S/N of the recovered HI angular power spectra. For vt = 0.5 deg s−1 and α = 2, nearly the entire 1/f noise parameter space has an S/N < 1, while for vt = 0.5 deg s−1 and α = 1, there is a reasonable detection of the HI spectra at all but the highest β and fk values. Figure 14. View largeDownload slide The sum in quadrature over all the S/N in each plot of Fig. 13, using equation (35). These figures are similar to those shown in Fig. 9, but include the mean residual power as an uncertainty contribution. As with Fig. 9, each value in these plots comes from summing over all the S/N grids, this time as shown in Fig. 13, for each set of noise parameter inputs. The top row is for α = 1, and from left to right increasing slew speeds of vt = 0.5, 1, and 2 deg s−1. The bottom row is same again but with 1/f noise with α = 2. All plots are post foreground cleaning. Figure 14. View largeDownload slide The sum in quadrature over all the S/N in each plot of Fig. 13, using equation (35). These figures are similar to those shown in Fig. 9, but include the mean residual power as an uncertainty contribution. As with Fig. 9, each value in these plots comes from summing over all the S/N grids, this time as shown in Fig. 13, for each set of noise parameter inputs. The top row is for α = 1, and from left to right increasing slew speeds of vt = 0.5, 1, and 2 deg s−1. The bottom row is same again but with 1/f noise with α = 2. All plots are post foreground cleaning. 5 DISCUSSION This section will present some discussion, informed from the results in Section 4, about 1/f noise in the context of HI IM surveys. To begin in Section 5.1, an empirical model for the 1/f noise angular power spectra measured in these simulations is presented. In Section 5.2, some general comments about how the impact of 1/f noise in an SKA HI IM survey may differ from the results presented here due to interactions with the chosen observing strategy, data analysis methods, or other systematics. Finally, Section 5.3 attempts to provide some practical intuitions about 1/f noise that can be inferred from Section 4. 5.1 An empirical 1/f noise model from power spectrum reconstruction By combining the 1/f noise angular power spectra and the discussions of the previous section, an approximate empirical model of the real spectra can be derived. The first step is to bring together all the constituent parts of the model that are known already. To simplify the analysis, it will be assumed that the 1/f noise has no dependence on frequency (i.e. Trx + Tcmb ≫ Tsky), the noise is Gaussian and therefore integrates down as $$1/\sqrt{T_\text{obs}}$$ as discussed in Section 4.2, and that there is no correlation in frequency such that β = 1 (however an estimate for β < 1 is given later). The ℓ-independent amplitude of the empirical 1/f noise model is defined as \begin{eqnarray*} A &=& \left( \frac{T_{{\rm sys}}}{21\,\mathrm{ K}} \right)^2 \left(\frac{f_\mathrm{ k}}{1\,\text{Hz}}\right)^\alpha \left(\frac{\delta \nu }{20\,\text{MHz}}\right)^{-1} \left(\frac{N_\mathrm{ t}}{200}\right)^{-1} \nonumber\\ &&\times \,\left(\frac{T_\text{obs}}{30\,\text{d}}\right)^{-1} \left(\frac{\Omega _\text{s}}{20500\rm{\,deg}^2}\right) 10^{2} \mu \text{K}^2, \end{eqnarray*} (39) where Tsys is system temperature normalized to 21 K to account for a mean ≈1 K sky contribution within the simulation, the knee frequency fk, the channel width δν, Nt the number of telescopes, Tobs the number of observing days, and Ωs the survey sky coverage. Equation (39) describes the amplitude of the 1/f noise angular power spectrum, but the slope of the 1/f noise angular power spectra shown in Fig. 5 is dependent on a combination of the 1/f noise α and the choice of observing strategy. The impact of the scan speed on the amplitude of the 1/f noise power spectrum is shown in Fig. 10 to be also dependent on the 1/f noise α. The best-fitting model that describes the relationship between α, vt, and the slope of the 1/f noise angular power spectra for the range of parameters explored in this work was determined to be \begin{eqnarray*} \log _{10}\left(\frac{F_\mathrm{ l}}{\mu {\rm K}^2}\right) &=& \log _{10}\left(\frac{A}{\mu {\rm K}^2}\right) + a\left[\alpha - 1\right] \nonumber\\ &&+\, b \sqrt{\alpha } \log _{10}\left(\frac{{v}_{\rm t}}{{\rm deg\,s}^{-1}}\right) - \sqrt{c \alpha } \log _{10}(\ell ),\nonumber\\ \end{eqnarray*} (40) where A is the parameter described by equation (39), α is the 1/f spectral index, vt is the telescope speed in deg s−1, and a, b, and c are fitted constants. The best-fitting values for the three parameters are: a = 1.5, b= −1.5, and c = 0.5. Fig. 15 shows a comparison between the empirical 1/f noise model given by equations (39) and (40) with the output 1/f power spectra from the numerical simulations. The 1/f noise knee frequency was fixed at 10 Hz for each spectrum shown. The bottom three spectra have α = 1, and the top three spectra have α = 2. Smaller variations in the total amplitude of the spectra demonstrate the impact of survey speed on 1/f noise. The empirical model only traces the overall functional form of the numerical simulations, and does not account for the small fluctuations seen in the numerical simulations. The empirical model is accurate to the 10 per cent level over the parameter space explored in this work for spherical harmonics of ℓ < 100. Figure 15. View largeDownload slide Comparison between the empirical 1/f noise angular power spectrum model given by equation (39) and the simulated numerical 1/f noise power spectra. The spectral index of the top three spectra is α = 2 and 1 for the bottom three. The survey speed is varied over 0.5 < vt < 2 deg s−1 and the knee frequency is fixed for all spectra at fk = 10 Hz. Figure 15. View largeDownload slide Comparison between the empirical 1/f noise angular power spectrum model given by equation (39) and the simulated numerical 1/f noise power spectra. The spectral index of the top three spectra is α = 2 and 1 for the bottom three. The survey speed is varied over 0.5 < vt < 2 deg s−1 and the knee frequency is fixed for all spectra at fk = 10 Hz. There are several limitations of this empirical 1/f noise model. The principle limitation of equation (40) model is determining how to account for the 1/f noise frequency correlations described by β. How β impacts the amplitude of the 1/f noise power in the recovered HI angular power spectrum will depend on many aspects of the observations such as the observing strategy, how the noise averages, the choice of component separation, 1/f noise filtering on long time-scales, and more. Further, β will also have an impact on the residual 1/f noise angular power spectrum bias that will add additional systematic uncertainty not accounted for here. Considering these provisos, it is possible to determine the fractional decrease in Fℓ predicted by these particular simulations by measuring the mean fractional change in the angular power spectrum bias with β. Through inspection of the data, the following two parameter model was determined to be a reasonable fit to the fractional change in bias \begin{equation*} \frac{F_\ell (\beta )}{F_\ell (\beta = 1)} = a \sin \left(2 \pi \beta \right) + \beta \end{equation*} (41) where the fitted parameter a = −0.16 for 0 < β < 1. Again, to reiterate, care should be taken when using this model in combination with equation (40) to extrapolate the impact of β on to different HI IM experimental designs. The forms of equations (40) and (41) are not physically motivated, are limited largely to the description of the observations simulated within this work, and give limited insight into the correlations between spectral index and slew speed. Still, one particular use of equations (40) and (41) could be to improve forecasts of future SKA HI IM surveys by adding predictions for the 1/f noise power spectrum into Fisher matrix analyses, such as those presented in Bull et al. (2015) and Pourtsidou (2016). One practical example of using the equation (40) is to predict how much additional observing time may be required for a range of receiver 1/f noise parameters. For simplicity only, the worst case of β = 1 is considered. As S/N is proportional to the total observing time, then the following relationship is true \begin{equation*} T_\text{obs} = T_0 \frac{{\rm S/N}_\mathrm{ w}}{{\rm S/N}}, \end{equation*} (42) where T0 is the observing time to achieve the signal-to-noise ratio of S/Nw given just white noise, Tobs is the observing time required to achieve S/Nw given that S/N was achieved in time T0 due to additional 1/f noise. Fig. 16 shows difference in observing times calculated for 1/f noise with 1 mHz < fk< 10 Hz and 0.5 < α < 2.5. This analysis assumed the full bandwidth (Δν = 450 MHz) of the instrument, the signal to be at the weighted mean redshift of $$\bar{z} = 0.2$$, and includes a sample variance contribution. The plot shows how important the spectral index of the 1/f noise is, as it can easily increase the effective statistical uncertainty in the HI angular power spectrum by an order-of-magnitude between α = 1 and 2. Consideration of the 1/f noise spectral index (and β) should be factored in when designing the instrumentation for an HI IM experiment. However, this may be challenging for multipurpose instruments such as the SKA as competing factors may take precedence. Figure 16. View largeDownload slide Increase in observing time due to 1/f noise calculated using the empirical 1/f noise model of equation (40) assuming β = 1 (solid lines), β = 0.5 (dashed lines), and β = 0.25 (dot–dashed lines). The T0 observing time assumes the same input parameters as the main simulations to calculate the white noise level. Figure 16. View largeDownload slide Increase in observing time due to 1/f noise calculated using the empirical 1/f noise model of equation (40) assuming β = 1 (solid lines), β = 0.5 (dashed lines), and β = 0.25 (dot–dashed lines). The T0 observing time assumes the same input parameters as the main simulations to calculate the white noise level. 5.2 Impact of specific data analysis methods and systematics An interesting question to ask is could the impact of the 1/f noise presented in this work be reduced in a real HI IM survey? The answer is almost certainly yes, the most obvious improvement would be to use a more sophisticated component separation technique in place of PCA. For example, the GNILC (Olivari et al. 2016) method can perform spatially localized component separation for individual needlet scales, which could suppress the large-scale 1/f noise bias discussed in Section 4. Another improvement would be to use a more advanced map-making method, such as Destriping (e.g. Sutton et al. 2010) or optimal map making (e.g. Natoli et al. 2001). However, suchmap-making methods at present are not optimized for HI IM data as they do not consider the spectral covariance of the noise or the sky signal. Without modifying existing map-making codes to consider spectral correlations, there is the possibility that the map-making procedure itself could induce spectral variations in to the 1/f noise spectrum that is effectively similar to increasing β. Further improvements could also be achieved by using a more carefully considered observing strategy. For example, this analysis does not utilize of any of instrument specific advantages of the SKA. The observing strategy presented in this work assumes that all telescopes are all continuously slewing at the same constant elevation of 55 deg. However, each dish or subsets of dishes could be given individual scanning strategies that are designed to maximize cross-linking of scanning tracks, or equalize integration times across the sky. Similarly, no attempt is made to use cross-correlation information between dishes or take advantage of any potential interferometric observations that are taken in parallel. 1/f noise can also be suppressed during the calibration process. This can be achieved by injecting a signal from a calibration diode into the receiver timestream at regular intervals. The use of a diode for calibration purposes can present its own challenges. For example, accurate calibration using a diode requires that over the interval between two diode injections the stability of the diode must be significantly better than the receiver system being calibrated. However, to calibrate the 1/f noise on shorter time-scales would require a brighter diode signal, which can rapidly become impractical. As an example, a simple estimate of the required calibration diode stability implies that calibrating an SKA receiver with a channel width of 50 MHz on relatively long time-scales of 100 s will require the diode stabilities of $$\frac{\Delta T_{\rm cal}}{T_{\rm cal}} \approx 10^{-4}$$ and a diode brightness as 25 K. Roychowdhury (in preparation) intends to explore this problem in greater detail. Unfortunately, for the real observations suppression of 1/f noise will be far more challenging due to the presence of other systematics in the data. Some of these systematics will be intrinsic to the 1/f noise such as different functional forms for the PSD of the frequency correlations, the non-Gaussianity in the distribution of the 1/f noise amplitudes, and non-stationary properties for α or fk. Other systematics will be intrinsic to the instrument or the observations such as the combined interaction of thousands of standing waves unique to each line of sight and each telescope, which in large numbers may exhibit 1/f noise-like properties in frequency and time. Also, how the data are processed, calibrated, and binned into maps from which HI power spectra are extracted are all important, as each step can unintentionally introduce additional frequency structure into the 1/f noise spectrum. In the context of this work, such effects could be seen as increase the effective β value of the 1/f noise, and as such will make it more challenging for component separation methods to recover the underlying HI intensity field. A final point to consider is the impact of the 1/f noise and other systematics on the recovery of the frequency space correlations in the HI signal. Though not the focus of this paper, it should be pointed out that such correlations in the HI intensity field, such as redshift space distortions, will be lost when using the data processing method outlined in this work. Therefore, it may be necessary to have separate analysis pipelines, or perhaps different observations altogether, for spatial and radial science objectives. 5.3 Intuitions for 1/f noise in HI IM surveys As the results and discussions in this work have shown, there is no simple way of determining the exact impact of 1/f noise on any given future HI IM survey. However, several important guidelines can be inferred. First, and most importantly, if 1/f noise is completely correlated across the bandpass and is not interacting with any other systematics then it can be perfectly subtracted from the data. This ideal scenario is not likely to be the case in real data, but the existence of this limit is encouraging as it sets an ultimate goal in terms of receiver frequency stability. Therefore, in the instance that 1/f noise is not fully correlated across the band the following guidelines should be considered: First, and most importantly, attempts to measure the α and β of SD HI IM receivers should be made in order to inform the planning of any future survey. The frequency correlations as described by β (or any other functional form) are critical to determine. Instruments with highly uncorrelated 1/f noise in frequency will find HI IM significantly challenging. Care should be taken to preserve the statistical properties of the 1/f noise frequency spectrum to avoid inadvertently increasing the effective β of the 1/f noise. The spatial power of 1/f noise integrates down as $$1/\sqrt{T_{\rm obs}}$$ (assuming Gaussainity and stationarity), but the initial rms of the 1/f noise fluctuations depends strongly on α, which in turn can significantly impact the observing time required to achieve a desired S/N. Scan speeds should be as fast as feasibly possible, with the aim of at least achieving a speed that matches the period of a slew to the knee frequency of the 1/f noise integrated across the channel width to be used in the final data analysis. The observing time of the experiment should be long enough such that the integrated angular power of the 1/f noise in the map is less than the HI angular power spectrum at all scales of interest. 6 CONCLUSION This work has presented for the first time a set of end-to-end simulations that describe the impact of 1/f noise on a potential future SKA HI IM survey. As HI IM surveys use spectroscopic measurements of the extragalactic HI signal, it was necessary to expand the existing models of 1/f noise to include correlated spectral fluctuations, which were parametrized in this work by β. It was found that in the ideal case that 1/f noise is completely correlated across the bandpass of spectroscopic receivers, then it is trivial for existing component separations methods to remove it. However, should the 1/f noise not be perfectly correlated and have a β > 0 then 1/f noise can have a significant impact on the recovery of the HI angular power spectrum through both increased statistical uncertainty and an intrinsic bias. Exactly how much impact the 1/f noise has depends on the degree of 1/f noise correlation, and a combination of the intrinsic statistical properties of the 1/f noise such as the knee frequency and α, and the choice of telescope observing strategy. As this is the first detailed exploration of 1/f noise in HI IM, numerous assumptions have been made to simplify the 1/f noise signal. However, in real data some of these assumptions are unlikely to be valid and therefore 1/f noise should be expected to be more challenging to remove than implied by these simulations. First the intrinsic properties of the 1/f noise are assumed to have Gaussian distributed amplitudes, stationary values of α and knee frequency, and a simple power-law functional form for the frequency correlations. None of these assumptions are necessarily true, and should be investigated during the pre-commissioning stages of SD IM experiments. Further, the interactions of the 1/f noise with additional systematics is not explored, however in real data, there will be standing waves, RFI, calibration errors, beam deconvolution errors, and many more potential instrumental systematics. Finally, given the above assumptions, this paper does provide some rough numerical guidelines for which a future SKA HI IM survey should attempt to achieve. First, a design consideration for receivers should be to make the 1/f noise as correlated as possible across the band (e.g. β → 0). Following this, the data analysis pipeline should be carefully designed as to not corrupt the receiver 1/f noise properties in order to avoid inadvertent increases in the effective value of β. It was found 1/f noise α can have a major impact on the recovery of the HI angular power spectrum, and ideally α would be minimized during the design of the receivers. However, this may not always be practical due to other design constraints, therefore careful consideration of how to mitigate α should be implemented into any HI IM data analysis pipeline. Finally, the impact of the knee frequency on an HI IM experiment is to increase the potential observing time required to achieve a desired S/N, for a 30 d SKA HI IM survey the knee frequency ideally would be fk < 10 mHz per 20 MHz channel. In terms of scanning strategy a slew speed of 1 deg s−1 should be sufficient. In conclusion, the analysis presented here finds that 1/f noise could be significantly detrimental to future HI IM surveys, but careful choices of observing strategy, data analysis, and component separation method should be sufficient to ensure successful recovery of the cosmological HI signal. ACKNOWLEDGEMENTS SH, CD, and RAB acknowledge support from an STFC Consolidated Grant (ST/L000768/1). SH and CD acknowledges support from an ERC Starting (Consolidator) Grant (no. 307209) under FP7. LCO acknowledges funding from CNPq, Conselho Nacional de Desenvolvimento Científico e Tecnológico − Brazil. YZM acknowledges the support from the National Research Foundation of South Africa with Grant nos 105925 and 104800. YZM acknowledges the NAOC-UKZN Computational Astrophysics Centre (NUCAC), University of KwaZulu-Nata, Durban, 4000, South Africa. The authors would also like to acknowledge Mario Santos for providing useful feedback during the development of this work. REFERENCES Alam S. et al. , 2015 , ApJS , 219 , 12 https://doi.org/10.1088/0067-0049/219/1/12 CrossRef Search ADS Albrecht A. et al. , 2006 , preprint (arXiv) Alonso D. , Bull P. , Ferreira P. G. , Santos M. G. , 2015 , MNRAS , 447 , 400 https://doi.org/10.1093/mnras/stu2474 CrossRef Search ADS Alves M. I. R. , Davies R. D. , Dickinson C. , Calabretta M. , Davis R. , Staveley-Smith L. , 2012 , MNRAS , 422 , 2429 https://doi.org/10.1111/j.1365-2966.2012.20796.x CrossRef Search ADS Anderson C. et al. , 2017 , in American Astronomical Society Meeting Abstracts #230 . p. 314.08 Aubourg É. et al. , 2015 , Phys. Rev. D , 92 , 123516 https://doi.org/10.1103/PhysRevD.92.123516 CrossRef Search ADS Bandura K. et al. , 2014 , Proc. SPIE , 9145 , p. 914522 https://doi.org/10.1117/12.2054950 Battye R. A. , Davies R. D. , Weller J. , 2004 , MNRAS , 355 , 1339 https://doi.org/10.1111/j.1365-2966.2004.08416.x CrossRef Search ADS Battye R. A. , Browne I. W. A. , Dickinson C. , Heron G. , Maffei B. , Pourtsidou A. , 2013 , MNRAS , 434 , 1239 https://doi.org/10.1093/mnras/stt1082 CrossRef Search ADS Berger P. , Oppermann N. , Pen U.-L. , Shaw J. R. , 2017 , MNRAS , 472 , 4928 https://doi.org/10.1093/mnras/stx2329 CrossRef Search ADS Betoule M. et al. , 2014 , A&A , 568 , A22 https://doi.org/10.1051/0004-6361/201423413 CrossRef Search ADS Bigot-Sazy M.-A. et al. , 2015 , MNRAS , 454 , 3240 https://doi.org/10.1093/mnras/stv2153 CrossRef Search ADS Bigot-Sazy M.-A. et al. , 2016 , in Qain L. , Li D. , eds, ASP Conf. Ser. Vol. 502, Frontiers in Radio Astronomy and FAST Early Sciences Symposium 2015 , Astron. Soc. Pac , San Francisco , p. 41 Bonaldi A. , Brown M. L. , 2015 , MNRAS , 447 , 1973 https://doi.org/10.1093/mnras/stu2601 CrossRef Search ADS Bowman J. D. et al. , 2013 , PASA , 30 , e031 https://doi.org/10.1017/pas.2013.009 CrossRef Search ADS Bull P. , Ferreira P. G. , Patel P. , Santos M. G. , 2015 , ApJ , 803 , 21 https://doi.org/10.1088/0004-637X/803/1/21 CrossRef Search ADS Caloyannides M. A. , 1974 , Journal of Applied Physics , 45 , 307 https://doi.org/10.1063/1.1662977 CrossRef Search ADS Catinella B. , Cortese L. , 2015 , MNRAS , 446 , 3526 https://doi.org/10.1093/mnras/stu2241 CrossRef Search ADS Chang T.-C. , Pen U.-L. , Bandura K. , Peterson J. B. , 2010 , Nature , 466 , 463 https://doi.org/10.1038/nature09187 CrossRef Search ADS PubMed Chapman E. et al. , 2012 , MNRAS , 423 , 2518 https://doi.org/10.1111/j.1365-2966.2012.21065.x CrossRef Search ADS Chen X. , 2012 , Int. J. Mod. Phys. Conf. Ser. , 12 , 256 https://doi.org/10.1142/S2010194512006459 Chon G. , Challinor A. , Prunet S. , Hivon E. , Szapudi I. , 2004 , MNRAS , 350 , 914 https://doi.org/10.1111/j.1365-2966.2004.07737.x CrossRef Search ADS Condon J. J. , 1992 , ARA&A , 30 , 575 https://doi.org/10.1146/annurev.aa.30.090192.003043 CrossRef Search ADS Crites A. T. et al. , 2014 , Proc. SPIE , 9153 , 91531W https://doi.org/10.1117/12.2057207 Croft R. A. C. et al. , 2016 , MNRAS , 457 , 3541 https://doi.org/10.1093/mnras/stw204 CrossRef Search ADS Dalcin L. D. , Paz R. R. , Kler P. A. , Cosimo A. , 2011 , Adv. Water Resour. , 34 , 1124 https://doi.org/10.1016/j.advwatres.2011.04.013 CrossRef Search ADS Datta K. K. , Choudhury T. R. , Bharadwaj S. , 2007 , MNRAS , 378 , 119 https://doi.org/10.1111/j.1365-2966.2007.11747.x CrossRef Search ADS de Oliveira-Costa A. , Tegmark M. , Gaensler B. M. , Jonas J. , Landecker T. L. , Reich P. , 2008 , MNRAS , 388 , 247 https://doi.org/10.1111/j.1365-2966.2008.13376.x CrossRef Search ADS DeBoer D. R. et al. , 2017 , PASP , 129 , 045001 https://doi.org/10.1088/1538-3873/129/974/045001 CrossRef Search ADS Dickinson C. , Davies R. D. , Davis R. J. , 2003 , MNRAS , 341 , 369 https://doi.org/10.1046/j.1365-8711.2003.06439.x CrossRef Search ADS Draine B. T. , 2011 , Physics of the Interstellar and Intergalactic Medium , Princeton Univ. Press , Princeton, NJ Eisenstein D. J. et al. , 2005 , ApJ , 633 , 560 https://doi.org/10.1086/466512 CrossRef Search ADS Emerson D. T. , Graeve R. , 1988 , A&A , 190 , 353 Fan X. , Carilli C. L. , Keating B. , 2006 , ARA&A , 44 , 415 https://doi.org/10.1146/annurev.astro.44.051905.092514 CrossRef Search ADS Fernández X. et al. , 2016 , ApJ , 824 , L1 https://doi.org/10.3847/2041-8205/824/1/L1 CrossRef Search ADS Gilden D. L. , Thornton T. , Mallon M. W. , 1995 , Science , 267 , 1837 https://doi.org/10.1126/science.7892611 CrossRef Search ADS PubMed Giovanelli R. et al. , 2005 , AJ , 130 , 2598 https://doi.org/10.1086/497431 CrossRef Search ADS Górski K. M. , Hivon E. , Banday A. J. , Wandelt B. D. , Hansen F. K. , Reinecke M. , Bartelmann M. , 2005 , ApJ , 622 , 759 https://doi.org/10.1086/427976 CrossRef Search ADS Haslam C. G. T. , Salter C. J. , Stoffel H. , Wilson W. E. , 1982 , A&AS , 47 , 1 Hauser M. G. , Peebles P. J. E. , 1973 , ApJ , 185 , 757 https://doi.org/10.1086/152453 CrossRef Search ADS Hellbourg G. , Chippendale A. P. , Kesteven M. J. , Jeffs B. D. , 2014 , in 2014 IEEE Global Conference on Signal and Information Processing (GlobalSIP) , IEEE , Atlanta, GA , p. 1286 https://doi.org/10.1109/GlobalSIP.2014.7032330 Hinton S. R. et al. , 2017 , MNRAS , 464 , 4807 https://doi.org/10.1093/mnras/stw2725 CrossRef Search ADS Hunt L. R. , Pisano D. J. , Edel S. , 2016 , AJ , 152 , 30 https://doi.org/10.3847/0004-6256/152/2/30 CrossRef Search ADS Irfan M. O. et al. , 2015 , MNRAS , 448 , 3572 https://doi.org/10.1093/mnras/stv212 CrossRef Search ADS Johnson J. B. , 1925 , Phys. Rev. , 26 , 71 https://doi.org/10.1103/PhysRev.26.71 CrossRef Search ADS Keating G. K. et al. , 2015 , ApJ , 814 , 140 https://doi.org/10.1088/0004-637X/814/2/140 CrossRef Search ADS Keating G. K. , Marrone D. P. , Bower G. C. , Leitch E. , Carlstrom J. E. , DeBoer D. R. , 2016 , ApJ , 830 , 34 https://doi.org/10.3847/0004-637X/830/1/34 CrossRef Search ADS Knox L. , 1995 , Phys. Rev. D , 52 , 4307 https://doi.org/10.1103/PhysRevD.52.4307 CrossRef Search ADS Kobayashi M. , Musha T. , 1982 , IEEE Trans. Biomed. Eng. , Vol. 6 , p. 456 Lah P. et al. , 2009 , MNRAS , 399 , 1447 https://doi.org/10.1111/j.1365-2966.2009.15368.x CrossRef Search ADS Lewis A. , Bridle S. , 2002 , Phys. Rev. , D66 , 103511 https://doi.org/10.1103/PhysRevD.66.103511 Li T. Y. , Wechsler R. H. , Devaraj K. , Church S. E. , 2016 , ApJ , 817 , 169 https://doi.org/10.3847/0004-637X/817/2/169 CrossRef Search ADS Lidz A. , Furlanetto S. R. , Oh S. P. , Aguirre J. , Chang T.-C. , Doré O. , Pritchard J. R. , 2011 , ApJ , 741 , 70 https://doi.org/10.1088/0004-637X/741/2/70 CrossRef Search ADS Loeb A. , Wyithe J. S. B. , 2008 , Physical Review Letters , 100 , 161301 https://doi.org/10.1103/PhysRevLett.100.161301 CrossRef Search ADS PubMed Loeb A. , Zaldarriaga M. , 2004 , Physical Review Letters , 92 , 211301 https://doi.org/10.1103/PhysRevLett.92.211301 CrossRef Search ADS PubMed Ma Y.-Z. , Scott D. , 2016 , Phys. Rev. D , 93 , 083510 https://doi.org/10.1103/PhysRevD.93.083510 CrossRef Search ADS Mandal S. , Arfin S. K. , Sarpeshkar R. , 2009 , Electronics Letters , 45 , 81 https://doi.org/10.1049/el:20092638 CrossRef Search ADS Mandelbrot B. B. , Wallis J. R. , 1969 , Water Resources Research , 5 , 321 https://doi.org/10.1029/WR005i002p00321 CrossRef Search ADS Masui K. W. , McDonald P. , Pen U.-L. , 2010 , Phys. Rev. D , 81 , 103527 https://doi.org/10.1103/PhysRevD.81.103527 CrossRef Search ADS Masui K. W. et al. , 2013 , ApJ , 763 , L20 https://doi.org/10.1088/2041-8205/763/1/L20 CrossRef Search ADS McQuinn M. , Lidz A. , Zaldarriaga M. , Hernquist L. , Dutta S. , 2008 , MNRAS , 388 , 1101 https://doi.org/10.1111/j.1365-2966.2008.13271.x Natoli P. , de Gasperis G. , Gheller C. , Vittorio N. , 2001 , A&A , 372 , 346 https://doi.org/10.1051/0004-6361:20010393 CrossRef Search ADS Neeleman M. , Prochaska J. X. , Ribaudo J. , Lehner N. , Howk J. C. , Rafelski M. , Kanekar N. , 2016 , ApJ , 818 , 113 https://doi.org/10.3847/0004-637X/818/2/113 CrossRef Search ADS Newburgh L. B. et al. , 2016 , in Proc. SPIE , 9906 , p. 99065X https://doi.org/10.1117/12.2234286 Noterdaeme P. et al. , 2012 , A&A , 547 , L1 https://doi.org/10.1051/0004-6361/201220259 CrossRef Search ADS Nyquist H. , 1928 , Phys. Rev. , 32 , 110 https://doi.org/10.1103/PhysRev.32.110 CrossRef Search ADS Olivari L. C. , Remazeilles M. , Dickinson C. , 2016 , MNRAS , 456 , 2749 https://doi.org/10.1093/mnras/stv2884 CrossRef Search ADS Paciga G. et al. , 2011 , MNRAS , 413 , 1174 https://doi.org/10.1111/j.1365-2966.2011.18208.x CrossRef Search ADS Padmanabhan H. , Choudhury T. R. , Refregier A. , 2015 , MNRAS , 447 , 3745 https://doi.org/10.1093/mnras/stu2702 CrossRef Search ADS Parsons A. R. et al. , 2010 , AJ , 139 , 1468 https://doi.org/10.1088/0004-6256/139/4/1468 CrossRef Search ADS Peebles P. J. E. , 1973 , ApJ , 185 , 413 https://doi.org/10.1086/152431 CrossRef Search ADS Pen U.-L. , Staveley-Smith L. , Peterson J. B. , Chang T.-C. , 2009 , MNRAS , 394 , L6 https://doi.org/10.1111/j.1745-3933.2008.00581.x CrossRef Search ADS Planck Collaboration XIII et al. , 2016a , A&A , 594 , A10 https://doi.org/10.1051/0004-6361/201525967 CrossRef Search ADS Planck Collaboration XXV et al. , 2016b , A&A , 594 , A25 https://doi.org/10.1051/0004-6361/201526803 CrossRef Search ADS Platania P. , Burigana C. , Maino D. , Caserini E. , Bersanelli M. , Cappellini B. , Mennella A. , 2003 , A&A , 410 , 847 https://doi.org/10.1051/0004-6361:20031125 CrossRef Search ADS Pourtsidou A. , 2016 , MNRAS , 461 , 1457 https://doi.org/10.1093/mnras/stw1406 CrossRef Search ADS Press W. H. , 1978 , Comments Astrophys. , 7 , 103 Pritchard J. R. , Loeb A. , 2008 , Phys. Rev. D , 78 , 103511 https://doi.org/10.1103/PhysRevD.78.103511 CrossRef Search ADS Pritchard J. R. , Loeb A. , 2012 , Rep. Prog. Phys. , 75 , 086901 https://doi.org/10.1088/0034-4885/75/8/086901 CrossRef Search ADS PubMed Remazeilles M. , Dickinson C. , Banday A. J. , Bigot-Sazy M.-A. , Ghosh T. , 2015 , MNRAS , 451 , 4311 https://doi.org/10.1093/mnras/stv1274 CrossRef Search ADS Rhee J. , Lah P. , Chengalur J. N. , Briggs F. H. , Colless M. , 2016 , MNRAS , 460 , 2675 https://doi.org/10.1093/mnras/stw1097 CrossRef Search ADS Riess A. G. et al. , 1998 , AJ , 116 , 1009 https://doi.org/10.1086/300499 CrossRef Search ADS Rubiño-Martín J. A. et al. , 2010 , Astrophysics and Space Science Proceedings Vol. 14, Highlights of Spanish Astrophysics V , Springer-Verlag , Berlin Heidelberg , p. 127 https://doi.org/10.1007/978-3-642-11250-8_12 Santos M. G. , Cooray A. , Knox L. , 2005 , ApJ , 625 , 575 https://doi.org/10.1086/429857 CrossRef Search ADS Santos M. G. , , Alonso D. , Bull P. , Camera S. , Ferreira P. G. , 2014 , in Heavens A. , Starck J.-L. , Krone-Martins A. , eds, IAU Symp. Statistical Challenges in 21st Century Cosmology , Cambridge University Press , Cambridge Vol. 306 , p. 165 https://doi.org/10.1017/S174392131401388X Santos M. , Alonso D. , Bull P. , Silva M. B. , Yahya S. , 2015 , Advancing Astrophysics with the Square Kilometre Array (AASKA14) , Proceedings of Science , Trieste , p. 19 Santos M. G. et al. , 2017 , preprint (arXiv:1709.06099) Scheuer P. A. G. , Williams P. J. S. , 1968 , ARA&A , 6 , 321 https://doi.org/10.1146/annurev.aa.06.090168.001541 CrossRef Search ADS Seiffert M. , Mennella A. , Burigana C. , Mandolesi N. , Bersanelli M. , Meinhold P. , Lubin P. , 2002 , A&A , 391 , 1185 https://doi.org/10.1051/0004-6361:20020880 CrossRef Search ADS Sethi S. K. , 2005 , MNRAS , 363 , 818 https://doi.org/10.1111/j.1365-2966.2005.09485.x CrossRef Search ADS Shaw J. R. , Sigurdson K. , Pen U.-L. , Stebbins A. , Sitwell M. , 2014 , ApJ , 781 , 57 https://doi.org/10.1088/0004-637X/781/2/57 CrossRef Search ADS Shaw J. R. , Sigurdson K. , Sitwell M. , Stebbins A. , Pen U.-L. , 2015 , Phys. Rev. D , 91 , 083514 https://doi.org/10.1103/PhysRevD.91.083514 CrossRef Search ADS SKA Collaboration Turner W. , Stringhetti L. , Braun R. , McPherson A. , 2016 , Technical Report SKA-TEL-SKO-0000008, SKA Phase 1 System Requirements Specification . SKA Smoot G. F. , Debono I. , 2017 , A&A , 597 , A136 https://doi.org/10.1051/0004-6361/201526794 CrossRef Search ADS Sutton D. et al. , 2010 , MNRAS , 407 , 1387 https://doi.org/10.1111/j.1365-2966.2010.16954.x CrossRef Search ADS Switzer E. R. et al. , 2013 , MNRAS , 434 , L46 https://doi.org/10.1093/mnrasl/slt074 CrossRef Search ADS Traficante A. et al. , 2011 , MNRAS , 416 , 2932 https://doi.org/10.1111/j.1365-2966.2011.19244.x CrossRef Search ADS van Haarlem M. P. et al. , 2013 , A&A , 556 , A2 https://doi.org/10.1051/0004-6361/201220873 CrossRef Search ADS Verheijen M. A. W. , et al. , 2010 , in ISKAF 2010 Science Meeting , Proceedings of Science , Trieste . p. 54 Villaescusa-Navarro F. , Alonso D. , Viel M. , 2017 , MNRAS , 466 , 2736 https://doi.org/10.1093/mnras/stw3224 CrossRef Search ADS Voss R. F. , Clarke J. , 1978 , J. Acoust. Soc. Am. , 63 , 258 https://doi.org/10.1121/1.381721 CrossRef Search ADS Wandelt B. D. , Hivon E. , Górski K. M. , 2001 , Phys. Rev. D , 64 , 083003 https://doi.org/10.1103/PhysRevD.64.083003 CrossRef Search ADS Wilkinson P. N. , 1991 , in Cornwell T. J. , Perley R. A. , eds, ASP Conf. Ser. Vol. 19, IAU Colloq. 131: Radio Interferometry. Theory, Techniques, and Applications , Astron. Soc. Pac , San Francisco . p. 428 Wilson T. L. , Rohlfs K. , Hüttemeister S. , 2009 , Tools of Radio Astronomy . Springer-Verlag , Berlin https://doi.org/10.1007/978-3-540-85122-6 Wolz L. et al. , 2017 , MNRAS , 464 , 4938 https://doi.org/10.1093/mnras/stw2556 CrossRef Search ADS Wyithe J. S. B. , Loeb A. , 2008 , MNRAS , 383 , 606 https://doi.org/10.1111/j.1365-2966.2007.12568.x CrossRef Search ADS Zwaan M. A. , van Dokkum P. G. , Verheijen M. A. W. , 2001 , Science , 293 , 1800 https://doi.org/10.1126/science.1063034 CrossRef Search ADS PubMed Zwaan M. A. et al. , 2004 , MNRAS , 350 , 1210 https://doi.org/10.1111/j.1365-2966.2004.07782.x CrossRef Search ADS APPENDIX A PCA SUBTRACTION OF PERFECTLY CORRELATED SIGNALS In the main text, it was stated that 1/f noise that is perfectly correlated in frequency (β = 0) was not included in the suite of simulations because in that instance PCA can perfectly subtract all the 1/f noise fluctuations. Recall from equation (13) that the model 1/f noise frequency correlations are described by a power law. The model is designed such that if β = 0, there is zero power in the 2D PSD for all Fourier modes apart from the zeroth mode that describes the temporal fluctuations. This implies that the frequency fluctuations can be described by a single zero-order Fourier mode or constant offset in the TOD. Removing such a constant is trivial. However, the simulations use a receiver model that multiplies the 1/f noise gain fluctuations with the measured Tsys, which imprints a power-law spectral curvature into the 1/f noise due to the sky signal. This implies that the 1/f noise is not perfectly correlated, but critically is described by only one functional form. The following demonstrates a simple example of how PCA can easily remove 1/f noise described by a single functional form. Imagine a simple map with Npix, two frequency channels, and a single sky component of arbitrary complexity in the spatial domain but can be described in frequency by a single function f(ν). Such a map can be written as \begin{equation*} {\bf m}(\Omega , \nu ) = A(\Omega ){\left(\begin{array}{c}f(\nu _1) \\ f(\nu _2) \end{array}\right)}, \end{equation*} (A1) where m(Ω, ν) is a function of direction Ω and frequency ν, where the amplitude of the spatial variations are described by A and the two frequency channels are presented as a vector containing f(ν1) and f(ν2). Next generate the covariance matrix of the map by taking the outer product \begin{equation*} \mathbf {C} \!=\! \mathbf {m}\mathbf {m}^{{\rm T}} \!=\! A^2(\Omega ) {\left(\begin{array}{c{@}{\quad}c}f^2(\nu _1) & f(\nu _1)f(\nu _2) \\ f(\nu _1)f(\nu _2) & f^2(\nu _2) \\ \end{array}\right)} \!=\! \mathbf {U}\Lambda \mathbf {U}^{\rm T}, \end{equation*} (A2) where the covariance matrix has been decomposed into the orthogonal matrix $$\mathbf {U}$$ of eigenvectors and diagonal matrix Λ that contains the corresponding eigenvalues. Solving for the eigenvalues λ can then be done by solving for the condition $${\rm det}(\mathbf {C} - \mathbf {I}\lambda ) = 0$$. The resulting determinant is then \begin{equation*} {\rm det}(\mathbf {C} - \mathbf {I}\lambda ) = \lambda \left(\lambda - f^2(\nu _1) - f^2(\nu _2) \right) \end{equation*} (A3) where the two possible eigenvalues are \begin{equation*} \lambda = \left\lbrace \sum _i^{{\rm N}_\nu } f^2(\nu _i), 0\right\rbrace , \end{equation*} (A4) or simply there is only one non-zero eigenmode that is the sum of squares of each frequency channel $$\lambda = \sum _i^{{\rm N}_\nu } f^2(\nu _i)$$. This statement holds true regardless of the number of frequency channels. Assuming then that the 1/f noise is described by a single power law in all pixels, then it will only need one eigenmode to describe it entirely. Since then the HI signal of interest is uncorrelated in frequency (assuming sufficiently wide redshift bins on the scales of interest), it requires Nchan modes to be fully described, and so can be easily separated from the 1/f noise signal. The effect then of increasing β is to increase the information content of the 1/f noise, with higher values of β requiring more eigenmodes to be fully described. Of course there are further complications such as each pixel may have different spectral indices, or in the real data, the 1/f noise may have a complex relationship with other systematics meaning more than one functional form is required to describe 1/f noise even with β = 0. In this case, 1/f noise will not be perfectly separated from the signal of interest, and the resulting residual 1/f noise power will depend on this effective value of β. APPENDIX B DERIVATION OF WHITE NOISE UNCERTAINTY This section will show how to derive equation (31). It is assumed in this derivation that there are two components of interest, one is the signal of interest and the other is the noise. The spherical harmonic coefficients (aℓm) are a linear combination of the signal and noise such that \begin{equation*} a_{\ell m} = s_{\ell m} + n_{\ell m} \end{equation*} (B1) where sℓm describes the signal coefficients and nℓm the noise coefficients. Taking the absolute square and the fourth power of equation (B1) results in \begin{equation*} \left| a_{\ell m} \right|^2 = \left| s_{\ell m} \right|^2 + \left| n_{\ell m} \right|^2 + 2\left| s_{\ell m} \right| \left| n_{\ell m} \right|, \end{equation*} (B2) and \begin{eqnarray*} | a_{\ell m} |^4 &=& | s_{\ell m} |^4 + | n_{\ell m} |^4 + 6| s_{\ell m} |^2 | n_{\ell m} |^2 \nonumber\\ &&+\, 4 \left[ |s_{\ell m}|^3 |n_{\ell m}| + |s_{\ell m}| |n_{\ell m}|^3 \right]. \end{eqnarray*} (B3) Both the signal and noise are Gaussian random variables sampled from distributions with variances Sℓ and Nℓ, respectively. Gaussianity implies that ⟨|aℓm|2⟩ = Cℓ and $$\left\langle | a_{\ell m} |^4 \right\rangle = 3 C_\ell ^2$$. Therefore, the average over equations (B2) and (B3) results in \begin{equation*} \left\langle | a_{\ell m} |^2 \right\rangle = S_\ell + N_\ell , \end{equation*} (B4) and \begin{equation*} \left\langle | a_{\ell m} |^4 \right\rangle = 3S_\ell ^2 + 3N_\ell ^2 + 6 S_\ell N_\ell , \end{equation*} (B5) where Sℓ and Nℓ are the angular power spectra of the signal and noise, respectively. Taking the square of equation (B3) from equation (B4) gives the uncertainty in the square of the angular coefficients \begin{equation*} \Delta C_\ell ^2 = \frac{2}{2 \ell + 1}(S_\ell + N_\ell )^2, \end{equation*} (B6) where Cℓ = Sℓ + Nℓ and the factor of 2ℓ + 1 comes from remembering that each Cℓ has 2ℓ + 1 m-modes that are independently drawn from the same random distribution. Equation (B6) is the familiar form of power spectrum uncertainty as stated in Knox (1995). To derive equation (31), it must be assumed that there is only one sky realization (as per the simulations described in this work) such that \begin{equation*} \left\langle | a_{\ell m} |^2 \right\rangle = \sqrt{\left\langle | a_{\ell m} |^4 \right\rangle } = C_\ell \end{equation*} (B7) where Cℓ is the measured variance of signal described by aℓm. Substituting equation (B7) for the signal in equations (B2) and (B3) results in \begin{equation*} \left\langle | a_{\ell m} |^2 \right\rangle = S_\ell + N_\ell , \end{equation*} (B8) and \begin{equation*} \left\langle | a_{\ell m} |^4 \right\rangle = S_\ell ^2 + 3N_\ell ^2 + 6 S_\ell N_\ell . \end{equation*} (B9) Differencing the square of equation (B8) from equation (B9) as before, again assuming 2ℓ + 1 independent m-modes, results in \begin{equation*} \Delta C_\ell ^2 = \frac{2}{2 \ell + 1} N_\ell ( 1 + 2\frac{N_\ell }{S_\ell })^2, \end{equation*} (B10) which is equation (31). © 2018 The Author(s) Published by Oxford University Press on behalf of the Royal Astronomical Society This article is published and distributed under the terms of the Oxford University Press, Standard Journals Publication Model (https://academic.oup.com/journals/pages/about_us/legal/notices)
TI - Impact of simulated 1/f noise for HI intensity mapping experiments
JF - Monthly Notices of the Royal Astronomical Society
DO - 10.1093/mnras/sty1238
DA - 2018-05-11
UR - https://www.deepdyve.com/lp/oxford-university-press/impact-of-simulated-1-f-noise-for-hi-intensity-mapping-experiments-QzmlnKqoxg
SP - 1
EP - 2437
VL - Advance Article
IS - 2
DP - DeepDyve
ER -