{'fraction_non_alphanumeric': 0.02989664839533035, 'fraction_numerical': 0.007395706933344491, 'mean_word_length': 4.629313390648923, 'pattern_counts': {'":': 1, '<': 1, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 1, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
The Epoch of Reionization Window: I. Mathematical Formalism

The 21 cm line provides a powerful probe of astrophysics and cosmology at high redshifts, but unlocking the potential of this probe requires the robust mitigation of foreground contaminants that are typically several orders of magnitude brighter than the cosmological signal. Recent simulations and observations have shown that the smooth spectral structure of foregrounds combines with instrument chromaticity to contaminate a"wedge"-shaped region in cylindrical Fourier space. While previous efforts have explored the suppression of foregrounds within this wedge, as well as the avoidance of this highly contaminated region, all such efforts have neglected a rigorous examination of the error statistics associated with the wedge. Using a quadratic estimator formalism applied to the interferometric measurement equation, we provide a framework for such a rigorous analysis (incorporating a fully covariant treatment of errors). Additionally, we find that there are strong error correlations at high spatial wavenumbers that have so far been neglected in sensitivity derivations. These error correlations substantially degrade the sensitivity of arrays relying on contributions from long baselines, compared to what one would estimate assuming uncorrelated errors.


I. INTRODUCTION
Modern cosmological observations have produced exquisite constraints on both the initial and final conditions of structure formation in our Universe. Initial conditions have now been probed to high significance with a large number of cosmic microwave background experiments [1,2], while at low redshifts, a combination of galaxy surveys and traditional astronomical measurements provide the final conditions [3]. Still missing from these direct observations, however, are the intermediate epochs that bridge the gap between early and late times. For example, despite tremendous recent progress in highredshift galaxy observations, details regarding the formation of the first luminous objects and their effects on the intergalactic medium (IGM) during the Epoch of Reionization (EoR) remain uncertain.
In the next few years, direct observations of the EoR will be made possible by measurements of the redshifted 21 cm hyperfine transition of neutral hydrogen (see, e.g., Refs. [4][5][6][7] for reviews). At the relevant redshifts, the intensity field of the 21 cm brightness temperature depends on a rich variety of different astrophysical effects, such as fluctuations in the ionization and spin states of the IGM, as well as cosmological quantities such as the underlying * acliu@berkeley.edu dark matter density field and peculiar velocity gradients. A map of the 21 cm intensity field at redshifts z ∼ 6 and above would therefore be a rich probe of EoR physics, including the nature of the first luminous sources (such as their typical mass and luminosity scales), their ionizing and heating efficiency, and feedback processes on the IGM, among other effects. Such a mapping can be accomplished in three dimensions, since the spectral nature of a 21 cm measurement provides redshift (and therefore line-of-sight distance) information, while the angular directions are mapped using traditional imaging. Because of this, the 21 cm line allows access to a large fraction of our Universe's comoving volume, potentially allowing futuristic measurements to move beyond astrophysics and into the measurement of fundamental cosmological parameters [8][9][10].
There are currently a number of experiments aimed at mapping the fluctuations of the cosmological 21 cm signal, including the Giant Metrewave Radio Telescope Epoch of Reionization experiment (GMRT-EoR [11]), the Low Frequency Array (LOFAR [12]), the Murchison Widefield Array (MWA [13]), and the Donald C. Backer Precision Array for Probing the Epoch of Reionization (PAPER [14]). These interferometer arrays have yet to make a positive detection of the cosmological signal, with the primary challenges being foreground contamination and the high sensitivity requirements. To increase sensitivity, these experiments are primarily targeting binned, statistical measures of the brightness temper-ature field such as the power spectrum. Recent progress has resulted in a number of increasingly stringent upper limits [11,15,16], and proposed next-generation instruments such as the Hydrogen Epoch of Reionization Array (HERA [17]) and the Square Kilometer Array (SKA [18]) promise to yield extremely high significance measurements.
In addition to achieving the required sensitivity, observations targeting the redshifted 21 cm line must also contend with foreground contaminants. In the relevant frequency ranges (roughly ∼ 100 to 200 MHz, corresponding to z ∼ 13 to 6), there exist a large number of noncosmological sources of radio emission that contaminate measurements. These include sources such as the diffuse synchrotron radiation from our own Galaxy, as well as extragalactic point sources, whether they are bright and resolved or part of a dim and unresolved background. The brightness temperature of foregrounds is expected to be 10 5 times greater than theoretical expectations for the amplitude of the cosmological signal. A detection of the reionization power spectrum will therefore be challenging without a robust foreground mitigation strategy.
Historically, cosmic microwave background (CMB) experiments have had to deal with similar problems of foreground contamination. However, strategies for foreground cleaning that have been developed for the CMB cannot be directly applied to 21 cm cosmology for two reasons. First, CMB experiments typically operate at higher frequencies, where foregrounds are not as bright. In fact, microwave-frequency foregrounds are subdominant to the CMB away from the Galactic plane. In addition, CMB experiments measure anisotropies over a twodimensional surface, with different observation frequencies providing consistency-checks and a set of redundant measurements that can be used for foreground isolation. The three-dimensional mapping of the 21 cm line, on the other hand, contains unique cosmological information at every frequency, which makes it more difficult to remove foregrounds in a way that does not result in the loss of cosmological signal [19].
Recently, however, a complication to this simple picture was realized, in what has been colloquially termed the "foreground wedge". Consider a cylindrically-binned FIG. 1. A schematic of the EoR window in the cylindrical k ⊥ k Fourier plane. At the lowest k ⊥ , errors increase because of limits on an instrument's field-of-view. High k ⊥ modes are probed by the longest baselines of an interferometer array, and the sensitivity drops to zero beyond k ⊥ scales corresponding to these baselines. Spectral resolution limits the sensitivity at large k . The lowest k are in principle limited by cosmic variance, but in practice the larger concern is limited bandwidth and the foreground contamination, which intrinsically resides at low k . As one moves towards higher k ⊥ , however, the foregrounds leak out to higher k in a characteristic shape known as the "foreground wedge". The remaining parts of the Fourier plane are thermal-noise dominated, allowing (with a large collecting area or a long integration time) a clean measurement of the power spectrum in this "EoR window". power spectrum measurement, i.e. one where Fourier amplitudes are squared and binned in annuli specified by wavenumbers perpendicular to the line-of-sight k ⊥ and wavenumbers parallel to the line-of-sight k . Because the line-of-sight direction is equivalent to the spectral axis of an interferometer, one might have naively expected smooth spectrum foregrounds to be sequestered to only the lowest k . However, this neglects the fact that interferometers are inherently chromatic instruments, with a given baseline probing finer spatial scales (higher k ⊥ ) at higher frequencies. This coupling of spectral and spatial information is sometimes coined mode-mixing, and results in the leakage of information from low to high k . This effect is particularly pronounced at high k ⊥ , where the modes are typically probed by longer baselines, which are more chromatic. Putting everything together, theoretical studies and simulations [32,[39][40][41][42][43][44] have shown that foregrounds are expected to leak out of the lowest k into a characteristic "wedge" that is schematically shown in Figure 1. Observations with PAPER and MWA have confirmed this basic picture [45], including its evolution with frequency [15].
The foreground wedge is both a blessing and a curse. At the sensitivity levels that have been achieved by current experiments, observations have seen a sharp drop-off in foregrounds beyond the wedge [45]. Theoretical calculations and simulations have shown that such a dropoff is the natural consequence of geometric limitations [32], provided the foregrounds are spectrally smooth. If further integration reveals low-level foregrounds that are spectrally unsmooth, their influence will leak beyond what is typically labeled as the edge of the wedge. However, if foregrounds continue to be reasonably smooth, the fact that physical considerations limit the extent of wedge implies that there must exist an "EoR window": a region in Fourier space that is a priori expected to be foreground-free. The existence of the EoR window thus enables a relatively robust foreground avoidance strategy, where a detection of the power spectrum can be made simply by avoiding measurements within the wedge. On the other hand, such a conservative approach forces one to work at higher k than if the chromatic effects had not caused the wedge in the first place. This is unfortunate because the ratio of the cosmological signal to instrumental noise typically peaks at low k, which means that if it were possible to work within the wedge (or to at least push back its boundaries a little), one would be able to make higher significance detections of the power spectrum. Indeed, in Ref. [17] it was suggested that working within the wedge can increase the detection significance anywhere from a factor of two to six (depending on the interferometer's configuration), with corresponding decreases in the error bars on astrophysical parameters.
Given the potentially high payoff associated with pushing back the influence of the wedge (or equivalently, enlarging the EoR window), it is important to have a statistically rigorous framework for describing the wedge. In this paper, we provide just such a mathematical framework. Since there already exists an extensive literature on the foreground wedge and the EoR window, it is worth summarizing the ways in which this paper builds upon and extends previous results. Works such as Refs. [39][40][41] describe instrumental simulations that take one particular realization of foregrounds and propagate them through a power spectrum estimation pipeline. They therefore only probe the mean power spectrum, and not the scatter (i.e. the errors) about this mean. In Ref. [43], a statistical treatment of point source populations was considered. While error bars were computed, offdiagonal error correlations (i.e. covariances) in the final measurements were neglected. Potential error correlations are important particularly because Ref. [43] considered the application of various tapering functions to their Fourier transforms, and certain choices can result in significant correlations between Fourier bins. We consider full foreground covariances, a full treatment of instrumental effects (such as having a non-tophat beam), a full treatment of data analysis choices such as tapering functions, and a fully covariant propagation of errors. We build upon Ref. [42], which made use of Monte Carlo methods to propagate errors. Our treatment is more analytic, allowing us to capture the large dynamic range needed to accurately compute the error statistics in a measurement where the foregrounds are many orders of magnitude brighter than the signal. This is made computationally feasible by our use of the delay spectrum approach introduced in Ref. [32], where input frequency spectra are sorted into a set of time-delay τ modes via a per-baseline Fourier transform. However, unlike in some works where the delays are used as an approximation for line-of-sight Fourier modes (an assumption that is only valid for short baselines, as we will discuss in Section II), we use delays strictly as a convenient choice of basis. This basis makes it computationally possible for us to deal directly with visibilities in our formalism (bypassing any mapmaking steps), which avoids gridding artifacts in our numerical results. We also take into full account the correlations between partially overlapping baselines, and therefore rigorously treat the possible complications that were highlighted in Ref. [44]. A related treatment pertaining to lower-redshift 21 cm intensity mapping experiments (though focusing less on the details of the wedge) can be found in Ref. [35]. While our fiducial calculations are centered around instruments targeting the EoR, the techniques developed in this paper are equally applicable to cosmological 21 cm at lower (or higher) redshifts.
We accomplish our goals by making use of the quadratic estimator formalism, which was adapted for 21 cm power spectrum measurements in Refs. [36,37], and applied to real data in Ref. [15]. However, appropriate "wedge effects" were not incorporated into the formalism, an omission that we rectify in this paper. Placing everything in the quadratic estimator formalism enables a systematic computation of the aforementioned error statistics, as well as a systematic study of the optimality (or lack thereof) of various power spectrum estimators. In a sequel paper (Ref. [46], henceforth "Paper II"), we will take advantage of this to examine the extent to which statistical methods can enlarge the EoR window.
With our fully covariant treatment, we find that the wedge is not simply a region of large foreground errors and biases, but also as a marker for error correlations: at k ⊥ values where the wedge is a dominant effect, the errors tend to be strongly correlated. With strongly correlated errors, the number of independently measurable Fourier modes is reduced, suggesting that previous sensitivity estimates (such as those in Refs. [17,47,48]) may be overly optimistic, particularly for arrays that make use of long baselines (such as LOFAR or GMRT). In fact, the rewards for working within the wedge may be overrated as a result of this, but of course this cannot be quantified without a rigorous way to compute the error statistics of the wedge-hence the present paper.
The rest of this work is organized as follows. In Section II we examine the measurement equation of an interferometer in detail, paying special attention to chromatic effects. This provides a first non-covariant preview of the foreground wedge, which we generalize to an approx-imate, but fully covariant description in Section IV, following a review of the quadratic estimator formalism in Section III. In Section V we discard the approximations made in Section IV in a full numerical implementation of our formalism. We summarize our conclusions in Section VII. Because a large number of mathematical quantities are defined in this paper, we provide dictionaries in Tables I and II for the reader's convenience.

II. A NON-COVARIANT PREVIEW OF THE FOREGROUND WEDGE
In this section, we will derive the Fourier-space foreground wedge from first principles. We begin with the visibility V measured by a baseline b at frequency ν: where θ is the angular sky position, I(θ, ν) is the sky temperature and A(θ/θ 0 ) is the primary beam, 1 with θ 0 denoting its characteristic width. For notational simplicity in this section, we will omit the instrumental noise contribution to the visibility, but of course it is always implicitly present. Our (arbitrary) convention for A is that it is dimensionless and is normalized so that A(0) = 1. In what follows, we will see that the foreground wedge arises from the fact that the product of ν and θ appears in the complex exponential. Fourier transforms in θ are therefore coupled to ν and vice versa, leading to the "modemixing" phenomena coined in Ref. [41] and ultimately the wedge.
To mimic the discreteness of frequency channels in a real instrument, we introduce a function γ that describes the response of a single frequency channel. Our measurement equation therefore becomes 1 In general, the primary beam will depend on frequency, although for some instruments (such as PAPER) the antennas are intentionally designed to minimize the frequency-dependence of the beam [14]. In this paper, we will neglect the frequencydependence, because our goal is not to provide results pertaining to a particular instrument, but instead to provide a rigorous understanding of how foregrounds enter an interferometric power spectrum measurement. Including a frequency-dependent beam makes many of our analytic manipulations more difficult, which obscures the key physical effects that give rise to the foreground wedge. We note, however, that our general strategy of incorporating an interferometer's measurement equation into the quadratic estimator formalism is one that is capable of including frequency-dependent beams, albeit at a slightly greater computational cost. Taking a Fourier transform of the frequency spectrum of data from a single baseline essentially amounts to taking a Fourier transform along a solid line in this figure. This is a good approximation to taking a Fourier transform along the "true" frequency axis for short baselines.
with B chan denoting the width of a frequency channel, and γ normalized such that ∞ −∞ γ(x)dx = 1. So far, we have no choice in the matter, in that the results of Eq. (2) are handed to us by the instrument. Moving onto data analysis, however, there is considerable freedom as to how one proceeds. For example, we may choose to move into delay-space, which is accomplished by taking the Fourier transform in frequency of the spectrum measured by a single baseline (a "delay transform"): where τ is the delay (with units of time), B band is the bandwidth over which we wish to compute a power spectrum, and ν 0 is the central frequency of our band. The function φ is normalized so that φ(0) = 1, and captures both the bandpass of our instrument and any tapering that one may wish to impose near the band edges. (In what follows, we will therefore use the terms "tapering function" and "bandpass" interchangeably to describe φ).
Precisely what form the edge tapering takes is a data analysis choice, and as shown in Ref. [43], choosing a good tapering function minimizes leakage of foregrounds in the Fourier plane. Later on, we will extend the work of Ref. [43] to self-consistently incorporate the effect that a tapering function has on error covariances. Aside from the peculiarities of certain tapering functions, working in delay space is simply a change of basis, and as emphasized in Ref. [32], represents no loss of generality. In later sections, we will find that delay space is a particular effi- Sky brightness temperature at angle θ and frequency ν Eq. Visibility response at delay τ of baseline b to sky mode on spatial scale u and spectral scale η. Eq. (17) g(u, η; b, τ ) Same as h(u, η; b, τ ) but integrated over direction on uv plane perpendicular to baseline vector Eq. (30) cient basis to work in, one that makes many of the large matrices required for power spectrum estimation sparse (see Appendix C for details). At this point our data are in a basis specified by baseline vector b and delay mode τ . This basis closely approximates the Fourier basis that the power spectrum inhabits, but the correspondence is not exact and in general must be accounted for. Concretely, suppose we let u be the Fourier dual of θ and η be the Fourier dual of ν. Since the transverse co-moving distance can be used to convert angles on the sky to transverse comoving distance, we have u ∝ k ⊥ (the exact conversion is given in Appendix A, where we define our Fourier conventions). Similarly, the observed frequency of a spectral line can be mapped to a co-moving line-of-sight distance, so η ∝ k .
At a particular frequency, a single baseline b of an interferometer roughly measures the spatial Fourier mode specified by wavenumber u = b/λ [this can be seen by applying Parseval's theorem to Eq. (1) and assuming that the primary beam is wide]. Across different frequencies, however, we see that a single baseline probes different spatial Fourier modes. A Fourier transform of a single baseline's spectrum [Eq. (3)] is thus a Fourier transform along one of the solid lines in Figure 2, where we show the chromaticity of various baselines. As pointed out in Ref. [32], with short baselines it is an excellent approximation to say that η ∼ τ . In this paper, we do not make this "delay approximation", as we consider baselines of TABLE II. Dictionary of vectors and matrices. The quantities shown here are grouped into three categories: those that exist in the vector space of the visibility measurements (indexed, for example, by baseline and delay), those that exist in the vector space of bandpowers (indexed by Fourier wavenumbers), and those that bridge the two spaces. In the column giving the length/dimensions, N bl denotes the number of baselines, Nν the number of frequency channels, and N bands the number of bins in Fourier space (i.e. the number of bandpowers). To properly account for the mapping between the "true" Fourier coordinates (u, η) to our visibilities expressed in the (b, τ ) basis, we can express the sky temperature I(θ, ν) in terms of its "true" Fourier transform, I(u, η): Inserting Eqs. (2) and (4) into Eq. (3) gives The integral over ν can be simplified because A is a much broader function (with a characteristic width of cb/θ 0 ) than γ (which has a characteristic width of B chan ). For a typical EoR experiment, one might have ν 0 = 150 MHz and B chan = 50 kHz. Even with a rather conservative θ 0 = 1 radian, A is wider than γ for any baseline shorter than ∼ 3000λ. With compact arrays giving the highest sensitivity [17,47], very few interferometers that are optimized for an EoR measurement have any relevant sensitivity on baselines this long. As a result, A can be factored out of 2 Note that the distinction between τ modes and η modes implies that our tapering functions differ slightly from those examined in Ref. [43]. In this paper, the tapering functions are applied to the per-baseline Fourier transform, rather than the Fourier transform along the "true" frequency axis of Figure 2.
the integral and evaluated at ν = ν, leaving a Fourier transform over the ν variable. Defining γ to be the Fourier transform of γ, this leaves With our approximations, then, the effect of having frequency channels with non-zero width is to envelope our response in η: limited spectral resolution makes the array less sensitive to high η modes (i.e., rapidly fluctuating spectral modes).
To further simplify our expression, it is useful to orient our u ≡ (u, v) axes so that the u axis is in the same direction as the baseline vector b. If we further assume that the antenna's footprint on the uv plane is separable, i.e.
then our expression becomes To proceed beyond this point requires specific forms for A b and γ. However, it is instructive to examine various limits. For an instrument with short baselines and/or narrow fields of view satisfying bθ 0 c/B band , A b is a slowly varying function that can be factored out of the integral, yielding where we used Eq. (7) to recombine A b and A b⊥ . This gives the "usual" description of an interferometric measurement: each baseline samples a portion of the sky in uvη space, defined by the antenna's footprint on the uv plane and the Fourier transform of the bandpass shape in η, enveloped by the Fourier transform of a frequency channel's profile. On the other hand, for an instrument with long baselines and/or wide fields of view satisfying bθ 0 c/B band , it is φ that is broad compared to A b , and Eq. (8) is well-approximated by 3 We see that in this limit, the bandpass shape φ and the primary beam A b swap roles: the bandpass acts as the convolution kernel on the uv plane, while the primary beam acts as the convolution kernel in the η-direction.
When bθ ∼ c/B band , the full expression given by Eq. (8) interpolates between the two extremes given by Eqs. (9) and (10); in general, the bandpass shape and the primary beam both play a part in determining the uv plane and ηdirection convolution kernels. This point was emphasized in Ref. [32]. Returning to the long baseline limit of Eq. (10), we can see our first glimpse of the foreground wedge. Suppose one were dealing with flat-spectrum foregrounds, where I(u, η) = I 0 (u)δ(η), so that there is no signal beyond η = 0. With such a sky, the measurement becomes pro- . Suppose we now make use of the delay approximation [32], where the quantity | V (b, τ )| 2 can be treated as an estimate of the power spectrum P (u, η) at u ≈ ν 0 b/c and η ≈ τ . The result is Thus, even flat spectrum sources (which would naively only have power at η = 0) gives a non-zero measured power spectrum at higher η. If the primary beam is zero beyond some argument value (which is always true in some sense, since A is identically zero below the horizon), then the power extends only to a finite region on the uηplane. For example, if we define θ 0 to be the angle at which A drops to zero (we are free to do so, since θ 0 is simply a characteristic scale that we have so far kept general as simply "some characteristic scale"), Eq. (11) predicts that the line should be a sharp boundary between zero and non-zero power. Switching from angular/spectral Fourier coordinates u and η to comoving spatial wavenumbers k ⊥ and k using the relations in Appendix A gives where H 0 is the Hubble parameter, Ω m is the normalized matter density, and Ω Λ is the normalized dark energy density. This is precisely the "usual" formula given for the edge of the foreground wedge (e.g., in Refs. [15,32,41,42]). As a function of η (or k ), the wedge has a profile given by the square of the primary beam profile, scaled by u in such a way that power is found at higher η for higher u.
In presenting a preview of the wedge, we have made a number of assumptions that will be relaxed in the following sections. For example, we will no longer approximate τ modes as η modes. Nor will we assume that each baseline cleanly samples just a single value of u. In fact, our use of the delay approximation here was somewhat inappropriate-while one may always take a delay transform, we have seen (e.g., from Figure 2) that the delay approximation (i.e., assuming τ ∼ η) ought to work well only if baselines are short. By using Eq. (10) instead of Eq. (9), however, we are expressly working in the longbaseline limit. It is thus important to emphasize that the preview shown here is included only to build intuition, and a much more rigorous treatment will be presented in later sections. There, we will also generalize to a sky with an arbitrary power spectrum. While there will be minor alterations to details of the foreground wedge, the basic picture will remain intact.
In the strictest sense, the derivation that we have just presented is nothing new. We have simply re-derived a number of results (e.g. the existence of the wedge) that are already known in the literature. However, our rederivations have been part of an analytic formalism, confirming a number of numerical results while bringing together many known (but previously separate) features in a unified framework. In the next section, we will extend this framework to include a fully covariant description of the wedge, including correlated errors in such a way that extends the quadratic power spectrum estimation techniques of Ref. [36] to include wedge physics. Setting up such a framework allows one to systematically examine the statistical properties of various power spectrum estimators in light of the wedge, a study that we perform in Paper II.

III. QUADRATIC ESTIMATOR FORMALISM
Until now, our focus has been on the measurement of a visibility, which is linear in temperature. The power spectrum, however, is a quantity that depends quadratically on temperature. In this section we very briefly review the mathematical machinery that makes possible the fully covariant description of the EoR window we will present in Sections IV and V. The basis of our discussion will be the quadratic estimator formalism, which has a long history in the CMB and galaxy survey literature (e.g., Refs. [49][50][51]), and was explicitly adapted for 21 cm cosmology in Refs. [36,37].
The central quantity that encodes our instrument (and therefore the statistical properties of our estimators) is the data covariance matrix C. To form the covariance matrix, imagine that our input data (organized by baseline and delay mode) is serialized into a data vector x, i.e., The covariance matrix is then given by C ≡ xx † . Importantly, we will keep track of all off-diagonal elements in the matrix, so that all correlations between different baselines and different delay (or spectral) modes are taken into account. Knowledge of these correlations will allow us to formulate both a covariant description of the foreground wedge and the tools to fight its contaminating influence.
Although we omitted noise contributions in the previous section for notational cleanliness, they of course contribute to the variance captured by C. Assuming that instrumental noise is uncorrelated with sky signals, the noise appears as an additive term to the sky covariance, so that we can write C as where S is the sky signal portion of the covariance. Inserting Eq. (6) into our general definition of the data covariance, we obtain where we have defined and have made use of the definition of the power spectrum P (u, η): As we have written it, our power spectrum is expressed in terms of u and η, instead of the more common combination of k ⊥ and k . Doing so minimizes the number of cosmological quantities in our expressions, since u and η are the more "natural" quantities from the perspective of the instrument. This represents no loss of generality, and indeed, in Section V we will express our results in terms of k ⊥ and k .
To form a practical estimator for the power spectrum, it is necessary to discretize. Assuming that the power spectrum possesses cylindrical symmetry 4 (so that the power depends only on the magnitudes of u and η), one can imagine binning the uvη space into a series of annuli, each specified by a radius |u| = √ u 2 + v 2 and "vertical distance" |η| away from the η = 0 plane. Each annulus can then be represented as a small cell on a |u|η plane (which we will henceforth call the uη plane to conform to convention). 5 As long as these cells are made sufficiently small, the power spectrum can be approximated by a constant bandpower p α in each cell, with α indexing a serialized list of locations on the uη plane. With this approximation, our covariance can be written compactly as where C ,α ≡ ∂C/∂p α is the response of the covariance to αth bandpower, and is given by (20) with V α denoting the annular volume that is binned into the αth uη cell. If desired, the sky contribution of the 4 The true power spectrum of the cosmological signal is of course spherically symmetric (thanks to statistical isotropy), and thus can be further binned. However, systematics such as foregrounds tend to be cylindrically symmetric, since the instrument probes the line-of-sight direction differently than it does the angular directions [52]. Another effect to consider (though it is beyond the scope of this paper) is that of redshift-space distortions, which will break spherical symmetry. The cylindrical power spectrum is therefore a useful diagnostic quantity to compute prior to the formation of a spherical power spectrum. Note that in general, we need not make any assumptions about symmetry in our estimation formalism. If desired, our formalism can be used to estimate P (u, η) without even the cylindrical binning step. The presence of symmetry, while useful for increasing signal-to-noise, is much less important than other assumptions such as the smoothness of foregrounds. (Although see Section VI for a discussion of strategies for dealing with unsmooth foregrounds). 5 We note an unfortunate bit of notation: u is used both to denote the first coordinate on the uv plane as well as the magnitude of the u vector. Unfortunately, both usages are standard.
covariance can be further divided into separate contributions from foregrounds and the cosmological signal: where C fg is the foreground covariance, and p sig α is the αth bandpower of the cosmological signal only. Eq. (21) is more general than Eq. (19), because the latter implicitly assumes that C fg is given by α p fg α C ,α , where p fg α is a set of foreground bandpowers. This assumption holds true only when the foregrounds are describable as a power spectrum.
In the quadratic estimator formalism, the bandpowers are extracted by forming weighted pairwise combinations of the data vector x. In particular, one can form a quadratic estimator p α of the true bandpower p α by computing 6 where E α is an estimator matrix of weights to be used for weighting pairwise products of the data when estimating the αth bandpower. To see how this works, consider (as a toy example) a noiseless cosmological survey with uncorrelated real-space measurements in a three-dimensional volume. Further suppose that one's goal is to measure the unbinned, three-dimensional power spectrum P (k), i.e., p α = P (k α ). If x is expressed in a real-space basis (so that it is simply a serialized list of real-space voxel intensities), a sensible choice for the estimator matrix would be E α ij ∝ e −ikα·(ri−rj ) , where r i and r j are the position vectors of the ith and jth voxels, respectively. We therefore see that in this example, the role of E α is to take a Fourier transform of the data. If one desires estimates of a binned power spectrum (for example, one where statistical isotropy allows the binning of power over shells of constant |k|), the relevant E α for different k α are simply averaged together in each bin. Now suppose instead that our data are expressed in a Fourier basis, so that each element of the data vector x represents the Fourier amplitude at some location in k space. The estimator matrix is then even simpler, and in fact becomes diagonal, with E α ij = δ αi δ ij . In our current application, the data are organized by baseline b and delay τ . As discussed above, b and τ closely approximate the Fourier wavenumbers u and η in some regimes, but the correspondence is not perfect. For our application we would therefore expect E α to be diagonal-dominant, but not be perfectly diagonal. (For an explicit form, see Sections IV and V).
In general, correlated errors and other instrumental effects (such as the ones that we seek to model in this paper) make the estimator matrices more complicated than they were in our pedagogical examples. They will typically involve the C ,α matrices, since those provide the link between the "input" space that the data covariance inhabits, and the "output" space of bandpowers. The detailed form of the family of E α matrices is a choice made by the data analyst, and different choices yield estimators with different statistical properties. One such property is the error covariance Σ αβ ≡ p α p β − p α p β of our estimated bandpowers, which is given by 7 a result that can be derived by direct substitution of Eq. (22) into the definition of the error covariance. The error bars on our bandpower estimates are given by ∆p α ≡ (Σ αα ) 1 2 , but it is important to note that the error covariance contains much more information than just the error bars: off-diagonal elements of the covariance encode correlations between different uη cells of the cylindrical power spectra. With a fully covariant formulation of the errors, it is possible to over-resolve in uη space, evading the commonly-made assumption that uη cells are independent (i.e., have a diagonal covariance matrix) so long as they are more than 1/B band apart in the η direction and separated by more than the width of A in the u direction. While this assumption is prevalent in the 21 cm cosmology literature for reasons of computational simplicity, we will see that it is one that should be avoided. More precisely, we will see in Section V that the errors on the uη plane are not independent at high k ⊥ (including the foreground wedge region), necessitating a full accounting of the entire error covariance matrix, and not just its diagonal elements.
In addition to the error covariance, the quadratic estimator formalism also allows the computation of biases and window functions. Taking the expectation value of Eq. (22), recalling that C ≡ xx † , and inserting Eq. (21) gives where we have defined the window function matrix 7 While we write all vector and matrix quantities in boldface, it is important to note that there are two different vector spaces at work here. The error covariance Σ and the window function matrix W defined later inhabit the "output" vector space indexed by locations on the uη plane, unlike matrices such as N and C, which inhabit the "input" vector space indexed by baseline and delay. Hybrid quantities include C,α and E α , which can either be thought of as a family of matrices in the input vector space, or as rank-3 tensors in a combined space. and the contamination bias From Eq. (24), we see that each row of the window function matrix gives a window function that specifies the linear combination of the true bandpowers that each estimate of a bandpower represents. Typically, the E α matrix is normalized such that each row of W sums to unity, allowing the linear combinations to be interpreted as weighted averages. 8 The contamination bias represents an additive bias to the estimated power spectrum that arises due to residual noise and foregrounds in the data. In practice, cross-correlation techniques-such as forming crosspower spectra between odd and even time samples of data, as was done in Ref. [15], or between different subsets of redundant baselines in an array, as was done in Ref. [16]-allow the noise bias to be eliminated without any explicit bias subtraction. The bias that one needs contend with is therefore solely comprised of the foreground bias: If a perfect foreground model is available, this expected level of this bias can be computed and subtracted from the power spectrum estimate. However, because a detailed knowledge of the low-frequency sky is as-yet unavailable, this subtraction step is often omitted to avoid over-subtractions that destroy cosmological information. Instead, one simply hopes that the bias is small in regions of the uη plane where one wishes to make a power spectrum measurement. In the following sections, the error covariance Σ, the window function matrix W, and the bias b α are the quantities that will provide us with a detailed, covariant picture of the EoR window and the foreground wedge. The bias essentially captures the power spectrum contribution from noise and foreground contaminants, and corresponds to the foreground wedge signatures seen in various simulations in the literature. The window functions provide an alternate view of the wedge: the wedge can be thought of as a leakage of power from low to high η (or equivalently, k ) modes that becomes increasingly pronounced as u (or k ⊥ ) increases. For a wedge to exist, window functions for bandpowers centered at high η and high u must have tails that extend to low η, where foregrounds live. Finally, the error covariance provides an estimate of the error bars throughout the Fourier plane (including the wedge region), as well as a quantification of how the chromatic response of an interferometer can cause error correlations between otherwise uncorrelated uη cells.

IV. A COVARIANT DESCRIPTION OF THE FOREGROUND WEDGE FOR A BASIC ESTIMATOR
In this section, we use the quadratic estimator formalism to derive the foreground wedge and the EoR window for a "basic" estimator of the power spectrum. We will make a number of approximations for the sake of analytical tractability, leaving an exact numerical treatment to Section V. The goal here is to formalize the discussion from Section II to obtain a fully covariant, analytic description of power spectrum properties at high u and low η, where the foreground wedge resides. This will provide a basic picture of the foreground challenges that we face, setting the stage for Paper II, where we look at how these challenges can be mitigated with better estimators.
The relatively simple estimator that we will examine in this paper is specified by the relation where N is the instrumental noise covariance and M α is a normalizing scalar 9 for each bandpower α. This choice gives an estimator that is quite similar to the crude | V (b, τ )| 2 estimator discussed informally in Section II. The principal difference between what we will consider here and our previous estimator is the presence of C ,α .
The role of C ,α is two-fold. Its first purpose is to complete a signal-to-noise weighting of our data: the copies of N −1 downweight noisy data, while C ,α (by virtue of its being the derivative of C) upweights the high signal portions. The second purpose is to map the data from the input baseline b and delay τ space to the output uη space. Recall from Section II that while b and τ approximate u and η, respectively, the correspondence is not perfect. Applying C ,α completes the transition to Fourier space, a fact that will become more apparent when we write down an explicit form for the matrix. Studying the basic estimator given by Eq. (28) is worthwhile because it is approximately equivalent to the methods used in a number of state-of-the-art 21 cm power spectrum pipelines for analyzing observations and simulations [43,44,53]. These pipelines typically use an 9 In this paper, we do not consider the more general possibility of a matrix-based normalization, where instead of a simple multiplicative normalization, one multiplies the unnormalized bandpower estimatesp β by a matrix M to form the normalized bandpowerspα, i.e.,pα = β M αβp β . For details, please see Ref. [15] or Paper II. optimal mapmaking approach [54,55] to first go from visibilities to a gridded uvη data cube of Fourier amplitudes. The complex magnitudes of these amplitudes are then squared and binned to estimate power spectrum bandpowers. (Note that while the mapmaking may be optimal in this case, the subsequent power spectrum estimation is not). In Appendix B, we will prove that in the limit of infinitely fine bins in Fourier space, such pipelines are equivalent to estimating power spectra directly from the visibilities using Eq. (28) and the quadratic estimator formalism. Our numerical results will therefore be roughly representative of those seen in the aforementioned pipelines, but with fewer gridding artifacts because we go straight from visibilities to power spectra.
For analytical tractability here and numerical tractability in later sections, we will use an approximate form for the covariance matrix: and we have once again omitted the instrumental noise contribution to the covariance for simplicity. Superficially, this looks quite similar to Eqs. (16) and (17). However, in this case we have gone beyond simply forming C ≡ xx † from Eq. (8), in that we have performed the integral over A b⊥ . This is not always permissible, and represents a subtle additional approximation: we have assumed that baselines that are similar in length but very different in orientation have a negligible correlation with each other, and that those with similar orientations are correlated as though they were identical in orientation.
In other words, we assume that although two baselines can be completely uncorrelated, partially redundant, or perfectly redundant in the direction of the baseline vector, 10 overlaps between baselines in the transverse direction are treated in a binary fashion, so that the overlap is either zero or perfect. This assumption was inherited from the derivation of Eq. (8), which required a re-orientation of the axes of the uv plane so that the u axis would lie along the direction of the baseline. While this can always be done for a single baseline, the covariance matrix C encodes correlations between different baselines, which may be oriented differently. It is thus strictly speaking incorrect to form a covariance matrix from Eq. (8), and in principle one should use Eq. (6) instead. For the purposes of intuition, however, we may continue with our approximate expression as long as we remember that distant baselines have negligible correlation.
As we previously mentioned, C ,α provides the crucial link between the input data and the output Fourier space. It therefore forms a crucial component of any error statistic. Working in the limit of a continuous (rather than discrete) set of bandpowers, we may differentiate C with respect to P (u α , η α ) to obtain where we used the fact that P (u, η) can be written as du α dη α P (u α , η α )δ(u − u α )δ(η − η α ). Inserting this into Eq. (28) provides a concrete example of the general proof of equivalence given in Appendix B. One sees that each copy of g acts on a noise-weighted copy of the data vector x. Examining Eq. (30) reveals that the action of g is to Fourier transform the delay spectrum back into the frequency domain, apply another copy of the tapering functions, grid the result at the appropriate location on the uv plane, and then Fourier transform in frequency again, before adjusting for frequency channel discretization by an additional weighting in η. This is precisely the procedure that one would follow with a mapmaking algorithm in uvη space [55]. The result is then squared to form a power spectrum. While this particular example may at first sight seem to render the delay basis obsolete (since the first action of g is to transform back to a frequency spectrum), it is important to remember that in a more realistic case, one may be unable approximate the bandpowers as being continuous. For example, at small u and η, bin sizes may be comparable to the values of u and η themselves. Many of the algebraic simplifications used in this section then become inapplicable, necessitating full numerical manipulations of the relevant matrices, which are typically more computationally efficient in the delay basis (as we discuss in Appendix C).
Continuing with our approximation scheme for this section, however, Eq. (32) is particularly convenient for computing our suite of error statistics because it is separable. Taking advantage of this, the window functions for our basic estimator reduce to where we defined the shorthand g α ≡ g(u α , η α ; b, τ ), and in the last step assumed that our delay bins were fine enough to be approximated as being continuous. This form for the window function matrix has a straightforward geometric interpretation: when estimating the αth bandpower, one probes a mixture of the true bandpowers; the amount of the βth band that is included in the estimate of the αth band is given by the overlap of our interferometer's response to the αth and βth bands. We stress, however, that this simple form does not hold when one considers more complicated estimators such as the ones that we will consider in Paper II.
Let us now compute some example window functions. Just as we did in Section II, we can gain some analytic intuition by working in the short and long baseline limits. For short baselines satisfying bθ 0 c/B band (or equivalently, if the time-delay τ of a signal between the antennas of the baseline satisfies τ 1/θ 0 B band ), we may invoke the same approximations that led to Eq. (9), and say that Inserting Eq. (34) into Eq. (33) and evaluating the integral over τ yields where φ 2 signifies the Fourier transform of φ 2 , not the square of φ. This expression is in line with what one might intuitively expect from interferometry: the spatial u-dependence of the window function is controlled by the primary beam (or more precisely, its Fourier transform), while the spectral η-dependence is controlled by the bandpass. On the other hand, with long baselines satisfying bθ 0 c/B band we can use the approximations that led to Eq. (10), obtaining Once again, we may insert this into Eq. (33) to get Now, consider first the η-dependence (rather than the u-dependence) of the window functions, since our principal worry is that smooth, low-η foregrounds might scatter to higher η because of the instrument's chromaticity. If we suppress the u dependence by setting u = u α = u β , we obtain Eq. (38) predicts that with our basic estimator, foregrounds should appear in the now-familiar wedge in uηspace. To see this, consider the additive foreground bias for a hypothetical single-baseline interferometer. Since we are concerned with foregrounds, the most relevant regions of the uη plane will be the low η regions. We may therefore safely ignore the γ 2 terms, since at low η they will be approximately unity anyway. If we imagine that such foregrounds are described by a power spectrum p fg α , the foreground covariance is given by α p fg α C ,α , and the foreground contribution b fg α to the bias in Eq. (26) is given by where we have suppressed the u dependence of the foreground power spectrum for notational cleanliness. Now, if we approximate the foregrounds as being completely comprised of flat spectrum sources, then P fg (η) ∝ δ(η), and our foreground bias becomes This equation provides a general mathematical form for the profile of the foreground wedge, as a function of η, and reduces to previously-derived special cases for the wedge profile in the limit of top-hat primary beamshapes and bandpasses [42]. Had we retained the φ 2 terms, they would have enforced the condition that u ≈ bν 0 /c. We therefore see that foreground contamination ought to leak to higher η for higher values of u, since those are probed only by larger b. Importantly, we emphasize that this foreground wedge feature is basis-independent, in that the final result makes no mention of delays. Our choice in this paper to express visibilities and covariances using a delay basis (rather than, say, a frequency basis, i.e. spectra) is a choice that is computationally convenient (see Appendix C for details), but is fundamentally an arbitrary one. The same is true regarding our indexing of measurements by baseline. [The baseline length b that appears in Eq. (40) is an expression of array configuration rather than basis choice; indeed, Eq. (33) shows that all baseline indices are summed over]. The foreground wedge is purely a function of our instrument's design and the form of our power spectrum estimator.
With multiple baselines, Eq. (38) contains cross-terms between different baselines. Working once again at low η and high u (long baselines), one has While this expression is certainly more complicated than the one we had before, the same basic picture holds: the window functions can be quite broad in the η direction, thus allowing foregrounds to be scattered from low to medium values of η. But in general, the resulting contamination is still enveloped by terms like (A b * A b ) c bθ0 η , which limits the possible contamination at high η.
Let us now turn briefly to the behavior of the window functions as a function of u. With the short baseline limit already provided by Eq. (35), we once again focus on the long-baseline limit. Setting η = η α = η β to isolate the u-dependence, Eq. (37) becomes where A 2 b denotes the Fourier transform of the square of A b , not the square of the Fourier transform. As expected, a baseline b roughly probes a u-scale equal to bν 0 /c. Additionally, the window functions peak at u α = u β , a fact that is enforced by the appearance of the two copies of φ as a product, as well as by the presence The width of the window function in the u direction therefore depends on both the primary beam and the bandpass function. The A 2 b term has a characteristic width of θ −1 0 , while the φ functions have a characteristic width of B band b/c. Now, the window function involves the product of these functions, which means that the width of its central portion will be determined mostly by the narrower of the two contributions. In the long baseline limit that we are working in, θ −1 0 B band b/c by definition, so the A 2 b part is the narrower contribution and determines the width of central portion of the window function. We thus predict that the central width is just θ −1 0 , and does not depend on the u (or equivalently, k ⊥ ) value on which the window function is centered.
However, simply characterizing the width of the central peak is insufficient for our purposes. As emphasized throughout this paper, the large dynamic range that exists between the bright foreground emission and the dim cosmological signal means that it is important to accurately capture the weak, low-level wings of the window functions, away from the central peak. These wings will be controlled by the broader contribution in Eq. (42), namely, the product of the bandpasses. As stated above, the bandpasses have a characteristic width B band b/c, and since baselines of length b probe spatial scales given by u ∼ ν 0 b/c, our window function wings will have a width ∆u of The widths of the window function wings are therefore proportional to u, and grow as u increases. Equivalently, since ∆u ∝ u, the fractional wing width is constant, and the wings will appear to have the same width on a logarithmic u (or k ⊥ ) scale. This is intuitively unsurprising, as longer baselines probe a greater spread of spatial scale due to their greater chromaticity, and the width of this spread is proportional to the baseline length (see Figure   2). Since a given u mode is mostly accessed by baselines of length b ∼ uλ, one is then driven to the conclusion that ∆u ∝ u.

V. A NUMERICAL MODEL OF A BASIC ESTIMATOR
Having made various approximations in the previous section to enable an analytic treatment of the foreground wedge, we will now discard most of these approximations in lieu of an exact numerical treatment of our basic estimator. We will find that the basic picture that we presented above remains unchanged.

A. Instrument and foreground model
The model instrument that we consider in this paper is intended to reflect a typical design for an interferometer optimized for 21 cm power spectrum (rather than one that is intended to function as a general-purpose lowfrequency radio observatory). Maximizing power spectrum sensitivity requires antennas to be placed in a way that yields a large number of short, identical baselines [17,47]. With this in mind, we perform our computations for a square, 20 by 20 array of antennas, with a 14 m spacing between adjacent antennas. Each antenna is assumed to have a Gaussian primary beam, with a fullwidth-half-max (FWHM) of 40.5 • that is approximated as frequency-independent in the B band = 8 MHz band that we consider. The frequency width of each individual spectral channel is set at B chan = 50 kHz. With no loss of qualitative generality, we consider only observations centered around ν 0 = 150 MHz. The formalism presented in this paper applies to all redshifts accessible to a 21 cm interferometer, and none of the "lessons learned" in our analysis are substantially changed by examining a different redshift.
(3)] to be Gaussian, even though previous studies in the literature have argued for more desirable choices such as Blackman-Harris function or a Blackman-Nuttal function [16,43]. Using a Gaussian allows us to compute analytically compute the ν integral in our measurement equation [Eq. (8)], giving: where the characteristic scale of the beam θ 0 is given by the standard deviation of our Gaussian beam, which in our case is θ 0 = FWHM/ √ 8 ln 2 = 17.2 • . Intuitively, one sees that each delay mode of each baseline probes a reasonably localized region in uvη space. One also sees that there exists a complex exponential term that mixes spatial and spectral information, which is to be expected given the chromatic nature of an interferometer's synthesized beam.
Following this, we form the covariance matrix under the same coherency approximation as the one we employed in the previous section: two baselines may have any amount of overlap on the uv plane in the direction parallel to their baseline vector, but are either completely non-overlapping or perfectly overlapping in the direction perpendicular to the baseline vector. This allows the integral over v (defined to be the direction on the uv plane perpendicular to a pair of correlated baselines) to be evaluated analytically. The sky signal portion of our covariance matrix is then where α i ≡ 2πθ 0 B band b i /c and similarly for α j . We take the frequency channel response γ to be a Gaussian, which makes its Fourier transformγ also a Gaussian. Computing C ,α is very similar to computing S. Since instrumental noise is random and does not depend on the power spectrum, we have C ,α = S ,α . To find C ,α , then, we simply need to evaluate the integrals in Eq. (45), but with u and η integration limits chosen to match to the band in question, rather than being −∞ and +∞.
Having described our instrument and how it manifests itself in the covariance of our measurements, the final ingredient that we require for our numerical calculations is a model for the total power spectrum P (u, η). We model the total power spectrum as the sum of the cosmological power spectrum and a foreground power spectrum. For the cosmological power spectrum, we use the spherically symmetric power spectrum provided in Ref. [56] and assume statistical isotropy to compute the cylindrical power spectrum needed for our covariance. As for the foregrounds, we consider a relatively simple twocomponent power spectrum model P fg : where A is an overall normalization, C diff is the angular power spectrum of the diffuse Galactic emission, and ν diff c is its frequency coherence length. The corresponding quantities for point sources are given by C ps =2πu and ν ps c . The η dependence of this parametric form is motivated by the mathematical results of Ref. [26]. In that work, empirically-motivated models of foreground spectra were put through an eigenmode analysis, and it was found that the resulting set of eigenmodes were essentially Fourier modes in frequency (i.e. η modes). The eigenvalue spectra were well-fit by a linear exponential in η, with a coherence frequency of 64.8 MHz. For simplicity, we adopt this value for both ν ps c and ν diff c .
To model the angular structure of the diffuse Galactic emission, we use the Global Sky Model (GSM) software [57] to generate a model of the sky at 150 MHz. We then compute the angular power spectrum of this model, which we find to be well-fit by with a 1 = −1.450, a 2 = 0.1003, b 1 = 0.7666, and b 2 = −2.365. For notational simplicity, we have omitted the overall normalization of C diff , as it can be absorbed into A. We note that while only the high (power-law) portion of the angular power spectrum is typically modeled in foreground studies, it is crucial to include the low behavior as well. To see this, note that Eq. (45) predicts a substantial overlap between the response of the shortest baselines and the u = 0 mode of the power spectrum. This is simply reflecting the fact that the sky is not infinite in extent. Thus, even though the autocorrelation/zero-spacing baseline products are typically discarded from interferometric data, the instrument may still be sensitive to the zero mode of the sky. For the point source contribution to foregrounds, we neglect clustering for simplicity, and therefore take C ps =2πu to be a constant. We fix this constant by assuming that at = 1000, the amplitude of the point source angular power spectrum is roughly a factor of 10 smaller than that of the diffuse emission [58]. We therefore set C ps =2πu = 0.1 C diff =1000 . Using a simple, -independent angular power spectrum for point sources is not a required approximation for our formalism, and this assumption can be easily relaxed. In principle, including clustering would boost the point source power at low modes [59]. In practice, however, the low regime is dominated by the diffuse Galactic emission anyway, and we do not expect that an inclusion of clustering would qualitatively impact our numerical results.
We fix the overall normalization A of our foregrounds by considering the zero-mode of the power spectrum. We require that where I GSM = 433 K is the mean temperature of our GSM foreground template. Finally, we must add the noise contribution to our covariance. The computation of a noise covariance matrix N is rather subtle, given the assumptions that we have made above regarding baselines that overlap in directions perpendicular to their baseline vectors. We first sort the baselines of our array by baseline length into 54 equallyspaced bins. If only one baseline fell into each bin, the noise variance assigned to each bin would be [47] N ii where t is the integration time (taken to be 520 hrs), T sys is the system temperature, and Note that our expression for the noise variance differs from equations that are commonly seen in the literature in two ways. First, our variance is proportional to B band . This is simply due to the fact that we are working in a delay basis rather than a frequency basis. In addition, the beam area used here is the integrated square of the beam profile (rather than just the integral of the beam profile itself), which was shown in Ref. [16] to be the correct beam area to use for this calculation. With our Gaussian primary beam, Ω pp = πθ 2 0 . We assume a skynoise dominated instrument and set T sys = I GSM .
We now adjust for the fact that each baseline length bin contains more than just a single baseline. If there are n(b i ) baselines within a particular bin, we simply divide the noise variance by n(b i ), assuming that a combination of instantaneous redundancy and rotation synthesis allow a large fraction of the baselines to be combined coherently, prior to forming a power spectrum. (In practice, this is a somewhat optimistic assumption, given that baselines that are dissimilar in orientation can only be combined statistically, and not coherently). Further assuming that the noise covariance matrix is diagonal, our final form for N is where the Kronecker delta function ensures that instrumental noise is uncorrelated both between baselines and between delay bins. Admittedly, the instrument and noise models presented in this section are quite crude and make use of a large number of simplifying assumptions. The assumptions regarding instrumental noise and rotation synthesis, in particular, are quite optimistic. However, the numerical computations that we perform in this paper are not designed to be definitive sensitivity calculations. Rather, the goal is to use our rough model-which captures the essential features of large-amplitude smoothspectrum foregrounds and lower-amplitude broadband instrumental noise-to gain some statistically rigorous intuition for the EoR window. For further technical details regarding our exact implementation (e.g. for information about bin sizes), we refer the reader to Appendix C.

B. Window functions and foreground bias
We now examine the statistical properties of our basic power spectrum estimator. We will find that with our basic estimator, the foregrounds appear in a wedge in k ⊥ k space, consistent with findings in the previous literature.
In Figure 3, we show the expected foreground bias from Eq. (27) in terms of the cylindrical Fourier coordinates k ⊥ and k . This shows the expected level of additive foreground bias in an estimate of the power spectrum. This foreground bias is the quantity that is most directly comparable to previous studies where a single set of simulated foregrounds are propagated through a power spectrum estimation pipeline. We see that our results are consistent with such studies as well as the analytic arguments of the previous section: the foregrounds are mostly sequestered to low k modes at low values of k ⊥ , rising to higher k modes at high k ⊥ in a wedge-like pattern. As predicted by Eq. (13), the edge of the wedge is defined by a line with unit logarithmic slope. Beyond this edge, the contamination drops off sharply to give a clean EoR window. Importantly, we emphasize that the characteristic shapes seen here are not tied to the specifics of our foreground model, beyond the fact that foregrounds are spectrally smooth. Instead, the foreground wedge is tied to the chromatic nature of our instrument and the properties of our basic estimator. While this may be difficult to establish definitively with a simulation, the covariant formalism that we make use of in this paper allows the foreground models to be easily disentangled from the way they interact with the instrument and the power spectrum estimation pipeline. For example, in Figure 4, we to using the more conventional Fourier coordinates k ⊥ and k when displaying our results. For details about this conversion and our Fourier conventions, please see Appendix A.
show some examples of window functions on the k ⊥ k plane. From Eq. (25), we know that these functions depend only on our choice of estimator (through E α ) and the instrument's response (through C ,β ). Thus, any signatures of the wedge that are independent of the foreground model should be apparent in the window functions.
Consider first the leftmost plot from Figure 4, which shows three window functions that are all centered on the k value of 0.25 hMpc −1 , but different k ⊥ values. As one moves to higher k ⊥ , the window functions become increasingly elongated 12 in the k direction. With long tails to low k (where foreground emission naturally resides), this implies a leakage of smooth foregrounds from low to high k . Since this effect is most pronounced at high k ⊥ , the result is precisely a wedge-like structure. To guide the eye, the black lines on each plot show the edge of the wedge as predicted 13 by Eq. (13) with θ 0 = π/2.
Window functions centered at higher k values (central and rightmost plots in Figure 4) also develop elongations 12 This elongation is real and not just an artifiact of our logarithmic k ⊥ k axes. We know this because the set of three windows shown in each plot of Figure 4 are chosen to be centered on the same k . 13 Our use of θ 0 = π/2 to define the edge of the wedge is somewhat arbitrary, given that there is nothing special about θ 0 = π/2 in our flat-sky approximation. However, as emphasized in Ref. [32], in a proper curved-sky treatment it is the natural scale to consider, since the primary beam of an instrument must vanish at the horizon.
as one moves from low to high k ⊥ , although the effect is visually subtle due to our logarithmic plotting. These elongations are slightly less important for smooth foregrounds, as even the elongated tails are not quite long enough to reach the lowest k for window functions that are centered at high k . However, such effects may be important if foregrounds turn out to contain unsmooth (high k ) components (see Section VI for a brief discussion of this). Except for at the lowest k ⊥ values, the width of the window functions in the k ⊥ direction also appears to increase with increasing k ⊥ . Plotted on the logarithmic axes of Figure 4, it is visually clear that at intermediate to high k ⊥ the window functions have a roughly constant logarithmic width, confirming the proportional increase in width with k ⊥ that we predicted in Eq. (43).
Our analytic, covariant treatment of power spectrum statistics allows us to compute window functions to the high dynamic range shown in Figure 4. This is crucial given that foregrounds are expected to be ∼ 10 10 times brighter in power (i.e., in temperature-squared units) compared to the cosmological signal. It is therefore essential to capture the low-level tails of window functions. Conveniently, once the windows have been computed, the foreground bias for a different foreground power spectrum can be easily determined using Eq. (39a).

C. Error bars and error covariance
In addition to the foreground bias, we may quantify the error covariance in estimates of the power spectrum. In Figure 5, we show the square-root of the diagonal elements of the error covariance matrix [Eq. (23)], i.e. the power spectrum error bars. These error bars capture more than just thermal noise errors, and include contributions from the foreground covariance. Indeed, one sees that the foreground wedge appears not just in the bias that we discussed previously, but also in the form of increased error bars. Outside the wedge, the error bars are dominated by thermal noise, and are quite low, in what constitutes the EoR window. Note that the thermal noise contribution did not appear in the bias, since we have assumed that cross-correlations have eliminated the noise bias. Towards the smallest k ⊥ , the errors increase by a small amount due to cosmic variance. This effect is typically negligible unless k is also small, but does play a small role since our model has a rather low noise level. (Again, we emphasize that our goal is not to perform a definitive sensitivity calculation). At the highest k , the errors also rise slightly. This is due to the finite spectral resolution of our instrument, which is self-consistently included in our error bars via the spectral channel profile γ in our formalism.
Beyond just the error bars, our formalism also delivers the off-diagonal elements of the error covariance Σ, which encode the error correlations between different k ⊥ k cells. As emphasized in Ref. [15], these correlations need to be quantified if one wishes to accurately propagate errors from the cylindrical power spectrum P (k ⊥ , k ) to the spherical power spectrum P (k). To capture this information, we consider the correlation matrix, defined as which is essentially a whitened version of the error covariance. Examining Σ instead of Σ allows the correlations rather than the larger errors within the wedge to be the dominant feature. To further aid visualization, we focus on just small portions of the correlation matrix. Whereas the full correlation matrix would relate all k ⊥ k coordinates to all other such coordinates, in Figure 6 we fix k ⊥ at three separate values, and consider the correlations between different k coordinates. Immediately obvious is the fact that the error correlations form qualitatively different structures at different values of k ⊥ . To calibrate our expectations, we include in each plot a semi-transparent square of size ∆k × ∆k , where This is the scale over which errors are expected to be correlated in the k direction. It is derived by making the assumption that for a survey with bandwidth B band , the error correlation scale in η should be roughly 1/B band , and expressing this scale in cosmological Fourier coordinates. In the leftmost plot of Figure 6, where k ⊥ is fixed at the low value of 0.0042 hMpc −1 , we see that our rough expectations are correct. Having chosen our Fourier cell sizes with the survey volume in mind (see Appendix C for details), the cells are seen to be essentially uncorrelated. This simple picture breaks down, however, as we move towards higher k ⊥ , where the effects of the foreground wedge are more pronounced. Figure 6 shows that increasing k ⊥ (which essentially means moving deeper and deeper into the wedge) causes different k cells to become increasingly correlated. To understand why error correlations are to be expected, consider what it would take for the errors to be uncorrelated in k . Suppose one had an achromatic instrument with noise properties that were uncorrelated and uniform (i.e. white) between all frequency channels. Moving into k space by way of a Fourier transform does not induce any error correlations, because uncorrelated white noise has the same statistical properties in all bases. On the other hand, for our measurement we have an inherently chromatic instrument, which makes our noise chromatic. Fourier transforming non-white noise will result in noise correlations, even if the noise was uncorrelated to begin with. Additionally, interferometers are chromatic in a very specific way, with longer baselines more chromatic (as illustrated in Figure  2). Since higher k ⊥ are probed by longer baselines, error correlations should increase with k ⊥ , as seen in Figure 6.
Viewed together, Figures 5 and 6 suggest that the characteristic k correlation scale coincides roughly with the k extent of the wedge. This suggests that the instrumental effects that caused the wedge have also decreased the number of independent measurements, thus decreasing overall signal-to-noise. (Note, however, that the correlated errors persist even at high k . Though somewhat challenging to see with our logarithmic axes and hybrid binning, one sees that even there the off-diagonal k correlations are greater for large k ⊥ modes). That there is a rough matching of scales between the wedge and our error correlation scale is not entirely surprising (though not a given in the a priori sense 14 ), since that window functions and covariances are closely related to one another. Recall that our window functions exhibited long elongations in k as we moved towards higher k ⊥ . With extremely broad window functions, our instrument essentially smoothed over a large number of modes, and it is unsurprising that the errors in nearby bins ended up being positively correlated.
We can quantify the error correlations in more detail by defining an effective number of independent cells N eff , where with N c being the number of Fourier cells that enter into a simple unweighted averaging of Fourier modes. By con- Effective number of independent cells N eff as a function of Nc, the number of Fourier cells included in an averaging of k modes (going from low to high k ). The approximately linear relationship seen for the low k ⊥ curve is indicative of uncorrelated modes. On the other hand, the high k ⊥ curve approaches linearity only when high k modes dominate the average. Its initially slow increase at low k shows that the modes are highly correlated. Thus, for power spectrum sensitivity calculations to be reliable, one must take error correlations into account or risk overestimating an instrument's sensitivity.
struction, if all the modes that one averages over are independent, N eff equals N c , whereas N eff = 1 if the modes are perfectly correlated. In Figure 7, we consider the effect of averaging Fourier modes along the k axis and show N eff as a function of N c , with N c increasing from unity (when only the lowest k mode is included in the average) to N c = 30 (the total number of k cells in our computation). We do this for two different constant k ⊥ slices on the Fourier plane, k ⊥ = 0.0042 hMpc −1 and 0.13 hMpc −1 (corresponding to the leftmost and rightmost plots of Figure 6, respectively). For low k ⊥ , one sees the linear relationship N eff ≈ N c regardless of how many k cells are included in the averaging. For high k ⊥ , however, N eff increases only very slowly at first, as one averages together the highly-correlated modes within the wedge. The increase in N eff is linear only at higher k , where our hybrid binning becomes logarithmic, and each cell encompasses a greater extent in k , eventually exceeding the error correlation length. In our formalism, the correlated errors within each of these larger cells have already been self-consistently averaged over, giving inde-pendent cells that contribute linearly to N eff . Regardless of the specific of one's power spectrum estimation formalism, it is crucial to take into account error correlations if averaging together bins that are narrower than the error correlation length, or if the bins are wider than this, to ensure that the implicit averaging performed within the bin is done correctly [rather than relying on possibly incorrect a priori assumptions such as that implied by Eq. (53)].
For the particular computational set-up in this paper, one sees about a factor of 2 reduction in N eff for the highest k ⊥ . Going to even higher k ⊥ causes even greater reductions in sensitivity (compared to simple expectations). The correlations discussed here are particularly important for experiments proposing to make measurements deep within the wedge. As discussed above, the extent of the wedge in k provides the characteristic error correlation length between different k cells. It therefore follows that errors are highly correlated whenever one chooses to work within the wedge. Though such measurements may be well-motivated by the fact that the cosmological-signal-to-thermal-noise ratio is largest at low k, previous studies that established this have typically assumed the correlation length given by Eq. (53). The resulting signal-to-noise ratios may therefore have been over-estimated.
We emphasize that the error correlations seen in this section exist whether or not one's power spectrum estimation pipeline includes a direct subtraction of modeled foregrounds from the input data. To see this, suppose a direct foreground subtraction scheme reduces the foreground covariance by some constant multiplicative constant, so that C fg → εC fg where 0 < ε ≤ 1. If one assumes that thermal noise is negligible compared to foreground residuals at low k after a long time-integration, the result of our reduced C fg will be a corresponding decrease in the amplitude of the final bias and the error bars. However, the plots of error correlation in Figure 6 will remain unchanged, since the correlation is insensitive to an overall scaling. The number of independent modes will therefore still decrease in the manner discussed above, reducing sensitivity. Since this loss of sensitivity is the most pronounced at high k ⊥ , it is particularly important to take into account for arrays that make use of long baselines, such as LOFAR or GMRT.

VI. UNSMOOTH FOREGROUNDS?
In many prior works on 21 cm cosmology, the assumption of spectrally smooth foregrounds is considered crucial to one's ability to perform foreground subtraction. At various points in this paper, we too have assumed that foregrounds are smooth, and incorporating this assumption into our general framework gave rise to various predictions, such as the existence of the EoR window. However, since the power spectrum estimation framework itself does not require smooth foregrounds, a number of our key results would survive a (hypothetical) discovery of unsmooth foreground sources. We now briefly discuss how the problem of foreground mitigation would change if such an unfortunate discovery were to be made.
Because the window functions encode only the mapping of the true power spectrum to the estimated power spectrum and do not depend on the actual power spectrum, they do not rely on the assumption of smoothness. Irrespective of whether the foregrounds are smooth or not, the window functions accurately describe how foreground power is smeared out on the k ⊥ k plane by our instrument. If the foregrounds are smooth, their influence is limited to the wedge. This makes foreground mitigation easy, as their avoidance requires no more than a simple cut on the Fourier plane. If the foregrounds are not as smooth as expected, the EoR window will be smaller, but its exact size can still be predicted by convolving the (now unsmooth) model of our foregrounds with the same window functions as before. Forthcoming data from various experiments at higher sensitivity will allow further foreground modeling, and-with the help of the window functions-an accurate determination of the extent of the foreground wedge. Encouragingly, recent theoretical calculations have shown that in physically-motivated models of synchrotron emission, foreground spectra tend to be smooth even under the most pessimistic of assumptions [60].
For the sake of argument, however, let us consider a worst-case scenario where foregrounds are discovered to be sufficiently unsmooth for the EoR window to be drastically reduced in size. In such a scenario, a number of strategies can be employed for foreground subtraction. First, foregrounds can be modeled and subtracted to the best of one's ability in the visibility data. Following that, a more sophisticated estimator (one that downweights the data not by N but by the total covariance C to account for uncertainty in the foregrounds) can be used. Finally, the foreground bias can be subtracted from the power spectrum using Eq. (39a), and window function decorrelation techniques can also be used in an attempt to increase the size of the EoR window. We explore a number of these techniques in Paper II [46].

VII. CONCLUSIONS
In any measurement of the redshifted 21 cm power spectrum, foreground contamination is a serious concern. Fortunately, observations and various theoretical studies have shown that despite complications arising from the inherently chromatic nature of an interferometric measurement, smooth spectrum foregrounds occupy a characteristic wedge region in cylindrical k ⊥ k Fourier space. The complement of this region is expected to be relatively foreground-free, forming an EoR window where measurements might be made.
While there exists an extensive literature on the topic, previous studies have typically focused on how the fore-ground wedge manifests itself in the mean power spectrum signal. However, the same physical effects that cause the wedge in the power spectrum also affect the associated error statistics, such as the error covariance and the window functions. An examination of some of these statistics was performed in Ref. [42] using Monte Carlo methods. In this paper, we have provided a complementary treatment by deriving a rigorous, fully-covariant mathematical description of the foreground wedge and the EoR window. While our methods require the numerical evaluation of some matrix expressions, they differ from previous work in that they do not require numerical simulations of interferometric measurements, since the underlying framework is largely analytic. This makes it possible to compute error statistics with very high dynamic range, which is crucial since the foregrounds are expected to dwarf both the instrumental noise and cosmological signal.
Our formalism takes advantage of the delay spectrum techniques introduced in Ref. [32] to achieve computational savings, and in fact it is the use of the delay basis that makes our covariant, high dynamic range calculations numerically feasible. However, we re-emphasize that this is merely a choice of basis, and that our results are independent of this choice. This was shown explicitly in Section IV, when we developed a description of the foreground wedge in terms of window functions. Our description decouples the causes of the wedge-which depend only on the chromatic nature of the instrument and the specific form of our power spectrum estimator-from the detailed nature of the foreground emission. Independent of foreground properties, window functions that are centered at high k ⊥ will typically develop long tails towards low k . The wedge then results from the additional assumption that foregrounds are spectrally smooth, so that strong signals from low k are transferred to higher k by the long tails. Once the window functions have been computed, however, our formalism allows such assumptions to be relaxed.
With a fully covariant framework, we are able to track all error correlations in our numerical computations. We find that measurements made at high k ⊥ have highly correlated errors, effectively reducing the number of independent measurements that can be made in that part of Fourier space. This is particularly important for sensitivity forecasts that rely heavily on measurements made within the wedge, since the wedge's extent in Fourier space is roughly on the same scale as that of the error correlations. Previous studies have typically neglected error correlations, assuming that errors are independent as long as the spatial Fourier cells are of the same size as an antenna's uv footprint, and the spectral Fourier cells are on the order of 1/B band . Our work suggests that this is likely to be too optimistic an assumption. At the highest k ⊥ considered in our numerical computations (k ⊥ = 0.13 hMpc −1 ), for example, error correlations reduce the number of independent modes by approximately a factor of 2. This effect will be even more pronounced at even higher k ⊥ , which are probed by experiments with extremely long baselines. Since the chromatic effects that caused the wedge are closely related to those that cause error correlations, it will be crucial in future research to address the question of exactly how far the wedge can be pushed back (or equivalently, how much one can expand the EoR window). In Paper II, we use the formalism of this paper to explore statistical methods for enlarging the EoR window [46].
In this paper, our goal was to provide a rigorous treatment of the wedge. Previous treatments have typically made different simplifying assumptions. These include neglecting partially redundant baselines, approximating delay modes as η modes, making assumptions about baseline length, assuming top-hat primary beams, neglecting binning artifacts, or assuming that errors are uncorrelated on the Fourier plane. Our framework discards all of these approximations simultaneously, and it is gratifying to see that the basic picture of the EoR window as a naturally foreground-free region of Fourier space remains unchanged. This bodes well for foreground avoidance efforts that aim to detect the EoR by working outside the wedge, making it possible for 21 cm cosmology to open a new window into the high redshift universe using only existing data analysis techniques, with even more transformative results possible with further advances that expand the EoR window.
where c is the speed of light, z is the redshift of observation, H 0 is the Hubble parameter, Ω m is the normalized matter density, and Ω Λ is the normalized dark energy density. The line-of-sight direction is more subtle. When forming a power spectrum, one typically includes only data from a relatively narrow range in redshift. Otherwise, cosmological evolution invalidates the assumption of a translation-invariant temperature field, which is needed in the definition of a power spectrum. What matters, then, is not the mapping between frequency and the total comoving line-of-sight distance, but rather, the local relation between differences in frequency ∆ν and differences in distance ∆r at the redshift of observation: where ν 21 ≡ 1420 MHz is the rest frequency of the 21 cm line. Having established this mapping we may (in what is perhaps an abuse of notation) recenter our coordinates so that ∆ν → ν and ∆r → r . Such a re-centering introduces a constant phase shift in our Fourier transforms, which has no bearing on quadratic statistics such as the power spectrum.
Making the identification between T (r ⊥ , r ) and I(θ, ν) using our mappings, we may compare Eqs. (A1) and (A4) to conclude that This allows the power spectra defined under the different Fourier conventions to be related to one another: In this paper, we will plot all numerical results in terms of cosmological Fourier coordinates k ⊥ and k , even though we perform all our computations in terms of u and η. We use WMAP9 cosmological parameters: Ω m = 0.28, Ω Λ = 0.72, and H 0 = 69.7 km/s Mpc [1].
thereby proving that the basic power spectrum estimator that we examined in Section IV is equivalent to one where one forms a uvη space dirty map from visibilities, and then squares the result. We stress that the proof that we have just presented is basis-independent, in the sense that our data vector need not be indexed by baseline and delay. For example, one may choose to deal with frequency spectra rather than delay spectra, in which case the data vector would be indexed by baseline and frequency channel. The resulting h α vectors would no longer be given by Eq. (17), but the proof shown here would be unchanged.
Crucially, the proof shown here assumed an infinitelyfinely discretized Fourier space. In practice, this will only be a good approximation on small scales (large Fourier wavenumbers), where the difference in wavenumber ∆k between neighboring discretized bins is small compared to the magnitudes of the wavenumbers k themselves. In Paper II we will formalize this assumption by importing the Feldman-Kaiser-Peacock approximation that is commonly used in galaxy surveys.
Importantly, we emphasize that while squaring a normalized uvη dirty map is a perfectly reasonable way to estimate a power spectrum, it is by no means optimal. Indeed, Eqs. (B4) and (B7) are provably non-optimal, and better estimators are explored in Paper II. Instead, the estimators considered here (and therefore in Section IV) are intended to be representative of simple, "first pass" methods [43,44,53], and their statistical properties provide basic pictures of the challenges that one faces.

Appendix C: Technical computational details
In this Appendix, we provide further details pertaining to the numerical computations described in Section V.

Binning
In the quadratic estimator formalism used in this paper, there are two sets of discretizations: a discretization of the input data and a discretization of the output power spectra.
For all the computations performed here, we have chosen to express the input data in a basis parameterized by baselines and delay. In other words, each component of the input data vector x corresponds to a different baseline and delay pair. For computational tractability, we bin baselines together into 50 linear bins of width 5 m, with the first bin centered about 10 m and the last bin centered about 255 m. (There are in fact some slightly longer baselines in our array. However, we discard them in our analysis for reasons of numerical stability). The delay bins are again linear, and range from −200 µs to 198.75 µs in 320 bins. This gives delay increments of 0.125 µs, equal to the natural bin size of 1/B band . Note that this equality is by no means a requirement. If computational resources are not a concern, it may be preferable to slightly over-resolve. This is perfectly legitimate despite the fact that the bins will no longer be independent, since the inclusion of the channel profile γ and our tracking of full covariance information allows the nonindependence to be self-consistently captured.
Given that a baseline of length b roughly probes modes with u ∼ bν 0 /c, and that a delay bin τ roughly probes η ∼ τ , a logical way to discretize the output uη space would be to use linear bins that matched the input bins (but with u bins scaled by an appropriate factor of ν 0 /c). Such a scheme would be the most appropriate for matching the specifications of the instrument. However, the cosmological power spectrum is expected to evolve on logarithmic k scales. Thus, a linear binning is computationally wasteful at high k, where a large number of bins are used to resolve a power spectrum that does not evolve very much. On the other hand, a logarithmic binning scheme that is appropriate at high k will tend to be computationally wasteful at low k, where one would be over-resolving the instrumental response. As a compromise, we use a hybrid binning scheme that is roughly linear at low k and roughly logarithmic at high k. In this scheme, the (n + 1)th boundary of the u bins u n+1 is given by u n+1 = 1.036u n + 2.5. (C1) Similarly, the (n + 1)th boundary of the η bins η n+1 is given by η n+1 = 1.095η n + 0.125 µs.
At low u and low η the additive terms dominate, yielding bin boundaries that are spaced in an approximately linear fashion well-suited to the instrumental specifications. At high u and η the multiplicative terms dominate, giving logarithmic bins that are a good fit for theoretical expectations. For both u and η, we use 30 bins, giving a total of 900 uη bandpowers. The bottom edge of the lowest u bin is at u = 3, while the bottom edge of the lowest η bin is at η = 0.12 µs.

Sparseness and computational shortcuts
The methods and computations presented in this paper are basis-independent. By this, we mean that while our final goal is to estimate a power spectrum (and its associated error statistics) on the k ⊥ k plane, our input data may be expressed in any basis that we find convenient. We now elaborate on our reasons for working in a baseline/delay basis.
As an example, consider the evaluation of C. With the Gaussian beams and tapering functions used in Section V, the covariance C is given by Eq. (45). For parts of the matrix corresponding to short baselines, one can see by inspection that the matrix will be diagonal-dominant, with a large number of off-diagonal elements that are close to zero. The C ,α matrices are even more sparse, since many of the diagonal elements (those that do not satisfy u ≈ ν 0 b i /c or η ≈ τ i ) will also be zero. In our computations, we skip the evaluation of matrix elements that are expected to be small. This represents significant savings in computation time, given that with our binning scheme each matrix measures 16, 000 × 16, 000, and each element requires numerically integrating a twodimensional integral given by Eq. (45). Moreover, with 900 bandpowers, this process must be repeated 900 times for each of the C ,α matrices.
In our implementation, we set off-diagonal matrix elements of C to zero if the integrand is suppressed by at least 10 −12 relative to the relevant diagonal elements everywhere over the integration volume. For the C ,α matrices, we apply the additional constraint that a diagonal element is to be skipped if the integrand is attenuated by 10 −12 or more compared to a different diagonal element that satisfies u α ≈ ν 0 b/c and η α ≈ τ , where u α and η α are the uη values corresponding to the αth band.
The sparseness that we have described here is a direct result of our using a baseline/delay basis. In contrast, parameterizing the spectral information in a frequency basis results in substantially denser matrices, since the data are highly correlated between frequencies [42]. The delay transform roughly isolates spectral information by η mode. This isolation is imperfect in the long baseline limit, as we saw in Eq. (10). The matrices are therefore still dense for elements corresponding to long baselines, but the sparseness that is available with short baselines provides enough savings to enable the full propagation of covariant information without resorting to Monte Carlo methods, which can sometimes be slow to converge to the dynamic range displayed in Section V. Finally, we note that in an application of our methods to real measurements, the delay transform operates individually on each baseline [32], and therefore can be applied with negligible computational cost to the input data.