Mitigating Internal Instrument Coupling for 21 cm Cosmology. I. Temporal and Spectral Modeling in Simulations

A reflection in the signal chain of an antenna inserts a copy of the original signal with a time lag. An example of this is a reflection at the end of an analog cable, where the signal travels back up the cable, reflects again at the start of the cable, and travels back down and is transmitted through the system along with the original signal. The time lag, or delay, the reflected signal has acquired is two times the cable length divided by the speed of light in the cable. The reflected signal also acquires an amplitude suppression, meaning that it is generally only a fraction of the input signal, but even a small fraction of the foreground signal in the data can dwarf the expected EoR signal and therefore needs to be accounted for. If v₁ is antenna 1's voltage spectrum without a signal chain reflection, then the presence of a reflection can be encapsulated as

$$\begin{eqnarray}&&{v}_{1}^{{\prime} }(\nu ,t)={v}_{1}(\nu ,t)+{\epsilon }_{11}(\nu ){v}_{1}(\nu ,t),\end{eqnarray} \tag{ 2 }$$

where ${v}_{1}^{{\prime} }$ is the voltage spectrum of antenna 1 with the reflection component, and is a coupling coefficient describing the reflection in antenna 1's signal chain (denoted as 11 because it is coupling the signal with itself). The coupling coefficient can be broken into three constituent parameters as

$$\begin{eqnarray}&&{\epsilon }_{11}(\nu )={A}_{11}{e}^{2\pi i{\tau }_{11}\nu +i{\phi }_{11}},\end{eqnarray} \tag{ 3 }$$

where A is the amplitude, τ is the delay offset (the total time it takes to be reflected), and ϕ is the phase offset the reflected signal may have acquired relative to the original signal. In Equation (2) we have assumed time and frequency stability of the reflection parameters, although in practice these parameters will have some variation with time and frequency.

If we insert the corrupted voltage spectra into the visibility equation (Equation (1)), we get

$$\begin{eqnarray}&&{V}_{12}^{{\prime} }={v}_{1}{v}_{2}^{* }+{\epsilon }_{11}{v}_{1}{v}_{2}^{* }+{v}_{1}{\epsilon }_{22}^{* }{v}_{2}^{* }+{\epsilon }_{11}{v}_{1}{\epsilon }_{22}^{* }{v}_{2}^{* }.\end{eqnarray} \tag{ 4 }$$

We can see that in addition to the original cross-correlation term (${v}_{1}{v}_{2}^{* }$) we now also have copies of it at positive and negative delay offsets that are suppressed in amplitude by a factor of A₁₁ and A₂₂, respectively. The time behavior of a reflection mimics that of the original data, in that it shows the same temporal oscillation (i.e., fringing) as the foregrounds, and thus also appears at the same fringe-rate modes as the foregrounds (e.g., right panel of Figure 2). The conjugation of ₂₂ means that the reflected signal from antenna 2 appears at negative delays in V₁₂, while the reflected signal from antenna 1 appears at positive delays.

The resultant autocorrelation visibility can also be computed, and is given by

$$\begin{eqnarray}&&{V}_{11}^{{\prime} }={v}_{1}{v}_{1}^{* }+{\epsilon }_{11}{v}_{1}{v}_{1}^{* }+{v}_{1}{\epsilon }_{11}^{* }{v}_{1}^{* }+| {\epsilon }_{11}{| }^{2}{v}_{1}{v}_{1}^{* },\end{eqnarray} \tag{ 5 }$$

where we see that the first-order reflections show up at $\pm {\tau }_{11}$, while the second-order reflection appears at τ = 0 ns, due to the conjugation of the coupling coefficient with itself. If we assume that the first-order reflections are at sufficiently high delay and neglect the second-order term, we can approximate the visibility at τ = 0 ns as ${\widetilde{V}}_{11}^{{\prime} }(\tau =0)\approx ({v}_{1}{v}_{1}^{* })(\tau =0)$. Similarly, near the reflection delay of τ₁₁ the autocorrelation visibility simplifies to

$$\begin{eqnarray}&&{\widetilde{V}}_{11}^{{\prime} }(\tau ={\tau }_{11})\approx {\epsilon }_{11}{v}_{1}{v}_{1}^{* }.\end{eqnarray} \tag{ 6 }$$

This means that one can estimate the reflection coefficient amplitude in delay space as

$$\begin{eqnarray}&&{A}_{11}=\displaystyle \frac{\left|{\widetilde{V}}_{11}^{{\prime} }(\tau =\pm {\tau }_{11})\right|}{\left|{\widetilde{V}}_{11}^{{\prime} }(\tau =0)\right|},\end{eqnarray} \tag{ 7 }$$

which will be useful when modeling reflection systematics.

If one can estimate their parameters from the data, reflections can be removed via standard (direction-independent) antenna-based calibration. In this paradigm, the raw voltage spectrum of antenna 1 corrupted by the instrument is related to its true value as

$$\begin{eqnarray}&&{v}_{1}^{\mathrm{raw}}={v}_{1}{g}_{1},\end{eqnarray} \tag{ 8 }$$

which when inserted into the visibility equation yields the standard antenna-based calibration equation,

$$\begin{eqnarray}&&{V}_{12}^{\mathrm{raw}}={V}_{12}{g}_{1}{g}_{2}^{* }=\langle {v}_{1}{v}_{2}^{* }\rangle {g}_{1}{g}_{2}^{* }.\end{eqnarray} \tag{ 9 }$$

The g term is called the antenna gain, and accounts for amplitude and phase errors introduced by the various stages of the signal chain from the feed all the way to the correlator. Note that this form of the calibration equation does not account for polarization leakage induced by cross-feed coupling, which is generally a higher order correction (Hamaker et al. 1996; Sault et al. 1996). By rearranging Equation (2) as

$$\begin{eqnarray}&&{v}_{1}^{{\prime} }={v}_{1}(1+{\epsilon }_{11})={v}_{1}{g}_{1},\end{eqnarray} \tag{ 10 }$$

we can see that signal chain reflections can be completely encompassed in this gain term, and hence corrected for by dividing the corrupted data by a gain constructed from an estimate of the reflection coefficient.

We now turn our attention to another systematic we refer to as antenna cross coupling, which acts to couple one antenna's voltage stream with another antenna's voltage stream before reaching the correlation stage. Note that our model for cross coupling is different than "capacitive crosstalk" created by the electric field of two parallel signal chains interacting with each other within cabling, receivers, and analog-to-digital conversion units, which is a common systematic for radio interferometers (Parsons & Backer 2009; Zheng et al. 2014; Ali et al. 2015; Patil et al. 2017; Cheng et al. 2018). There are well-established hardware solutions for suppressing crosstalk, such as phase switching (Chaudhari et al. 2017). However, if residual crosstalk remains, or if phase switching is not implemented in the system, we need to model and remove it for robust EoR measurements.

Our cross coupling systematic model simply states that, before correlation, one antenna's voltage is added to another antenna's voltage with a coupling coefficient that determines the amplitude and relative delay with which the voltage is added. For the purposes of our description, we additionally assume that this coupling coefficient can be decomposed into the same three parameters as before, technically making it a form of reflection systematic.⁵ While this model may indeed be capable of describing certain some forms of capacitive crosstalk, we do not expect all forms of capacitive crosstalk to necessarily fall within the bounds of these assumptions.

To write down how this affects the interferometric visibility, we can start by writing the corrupted antenna voltages as

$$\begin{eqnarray}\begin{array}{rcl}{v}_{1}^{{\prime} } & = & {v}_{1}+{\epsilon }_{21}{v}_{2}\\ {v}_{2}^{{\prime} } & = & {v}_{2}+{\epsilon }_{12}{v}_{1},\end{array}\end{eqnarray} \tag{ 11 }$$

where ₂₁ describes the voltage coupling of antenna 2 into antenna 1 and vice versa for ₁₂. Substituting these equations into Equation (1), we get

$$\begin{eqnarray}&&{V}_{12}^{{\prime} }={v}_{1}{v}_{2}^{* }+{v}_{1}{\epsilon }_{12}^{* }{v}_{1}^{* }+{\epsilon }_{21}{v}_{2}{v}_{2}^{* }+{\epsilon }_{21}{v}_{2}{\epsilon }_{12}^{* }{v}_{1}^{* }.\end{eqnarray} \tag{ 12 }$$

We can see that the cross-correlation visibility now contains the autocorrelation visibility terms ${v}_{1}{v}_{1}^{* }$ and ${v}_{2}{v}_{2}^{* }$ at the first-order level, which are purely real quantities and thus have identically zero phase. In the complex plane, the cross-correlation term ${v}_{1}{v}_{2}^{* }$ winds around the origin as a function of time because its phase varies temporally. Because the autocorrelation has no phase, the act of the cross-coupling terms is to introduce an additive bias to the data with an arbitrary phase set by the coupling coefficient itself. Assuming the coupling coefficient is slowly variable (if not completely stable), the first-order systematic terms in Equation (12) only change in amplitude over time set by the natural variation in the amplitude of the autocorrelation, (e.g., ${v}_{1}{v}_{1}^{* }$). This variation is generally fairly slow on timescales of a beam crossing, which for HERA is roughly 1 hr. This leads us to two critical insights about the behavior of the cross-coupling terms: (1) their time variability is slow, thus occupying low-fringe rate modes (e.g., see right of Figure 2) and (2) they have a time-stable phase determined solely by the phase of the coupling coefficient.

In a more generalized case, the cross coupling between antenna 1 and 2 may have an angular dependence on the sky. Take for example the case of mutual coupling (or feed-to-feed reflections), where part of the radiation incident on antenna 1's feed is reflected and received by antenna 2's feed. This behavior will be highly angular dependent due to the nontrivial electromagnetic properties of the feed itself. Nonetheless, we can reason that the systematic phenomenology will be similar to as before. We can think of this angular dependence as a windowing function on the primary beam of the underlying autocorrelation, meaning that the first-order terms in Equation (12) will be proportional to only a fraction of ${v}_{1}{v}_{1}^{* }$, such that they have a smaller amplitude. It may also mean that these terms will have a slightly faster time dependence in the data, as the "effective beam" created by the angular windowing function is smaller on the sky than the total primary beam, and thus leads to a faster "effective beam crossing time."

We can also compute the effects of cross coupling on the measured autocorrelation visibility, V₁₁, which yields

$$\begin{eqnarray}&&{V}_{11}^{{\prime} }={v}_{1}{v}_{1}^{* }+{v}_{1}{\epsilon }_{21}^{* }{v}_{2}^{* }+{\epsilon }_{21}{v}_{2}{v}_{1}^{* }+| {\epsilon }_{21}{| }^{2}{v}_{2}{v}_{2}^{* }.\end{eqnarray} \tag{ 13 }$$

In this case, we find that the cross correlation is inserted into the measured autocorrelation at the first-order level with a delay offset of τ₂₁. These terms are likely many order of magnitudes below the peak autocorrelation visibility amplitude, given that the cross-correlation visibilities are generally a few orders of magnitude below the autocorrelation inherently, which is further compounded by the amplitude suppression from ₂₁.

We can see simply from Equation (12) that the corruption of ${V}_{12}^{{\prime} }$ by cross coupling cannot be factorized into antenna-based gains, based simply on the presence of the ₁₂-like terms, which are baseline dependent. Removal of cross-coupling terms in the data must therefore be done on a per-baseline basis by constructing a model of the systematic in each visibility and then subtracting it.

To summarize, reflections along a single antenna's signal chain produces a duplicate of the signal with suppressed amplitude and some delay offset. This is true for both the cross-correlation and autocorrelation visibility products. Example mechanisms include cable reflections and dish-to-feed reflections within the confines of a single antenna. Reflections in the cross-correlation visibility have the same time structure as the un-reflected visibility, meaning reflected foreground signal occupies the same fringe-rate modes as un-reflected foreground signal, but is shifted to high delays (e.g., Figure 2). Reflections can be removed from the raw data by creating a model of the reflections and incorporating them into the per-antenna calibration gains.

Another systematic we describe is created by antenna-to-antenna cross coupling, which mixes the voltage signals between the antennas. This has the effect of introducing a copy of the autocorrelation visibility into the measured cross-correlation visibility at positive and negative delay offsets, and similarly introduces copies of the cross-correlation visibility into the measured autocorrelation visibility. In the measured cross-correlation visibility, the first-order coupling terms are slowly time variable and occupy low-fringe-rate modes centered at f = 0 Hz. Cross-coupling terms cannot be removed via antenna-based calibration, and must be modeled and subtracted at the per-baseline level.

Modeling signal chain reflections can in theory be done simultaneously with standard gain calibration because reflections factor as an antenna-based effect (see Section 2). However, there are reasons why we might be wary of using standard calibration techniques for deriving reflection parameters. Principally, standard bandpass calibration typically operates on the $\sim {N}_{\mathrm{ant}}^{2}$ number of cross-correlation visibilities, and generally allows each frequency channel's gain to be solved independently from other channels. Frequency-dependent calibration errors will therefore set a fundamental floor to the precision with which reflections can be calibrated via standard antenna-based bandpass techniques (Barry et al. 2016; Ewall-Wice et al. 2017; Orosz et al. 2019). Furthermore, the dynamic range of the signal to noise is considerably higher in the autocorrelation than in the cross-correlation visibility, as it is a measurement of the total power received by an antenna: in many cases a signal chain reflection cannot even be seen above the noise floor in a cross-correlation visibility but is highly apparent in the autocorrelation visibility. This latter point is important, because it implies that reflection parameters estimated from the autocorrelation visibilities will have a signal-to-noise ratio (S/N) that is drastically higher than the S/N of the cross-correlation visibilities, meaning that, when calibrating out systematics in the cross-correlation visibilities, the S/N of the derived reflection parameters will never be a limiting factor. Of course other real-world factors can limit the precision of the derived reflection parameters such as nontrivial frequency evolution, which is discussed in more detail in the context of HERA in Kern et al. (2019).

Reflection parameters must be estimated to high precision in order to get even a modest suppression of their systematic power in the visibilities. For example, Ewall-Wice et al. (2016) employed a reflection fitting algorithm on MWA autocorrelation visibilities by fitting sinusoids in the frequency domain and was able to suppress reflection systematics by a couple orders of magnitude in the power spectrum, although their end-result band powers were still systematic limited at some k modes. Similarly, Beardsley et al. (2016) explored reflection calibration on MWA data as an extension to their restricted polynomial gain calibration scheme and also achieved a couple of orders of magnitude of suppression in the 2D power spectrum, although their more deeply integrated power spectra show its reemergence.

The algorithm we present here also operates on the autocorrelation visibility, but we choose to model the reflection in the delay domain rather than the frequency domain. Recall from Section 2 that in the autocorrelation visibility, V₁₁, a signal chain reflection appears as a shifted copy of the original visibility at symmetric positive and negative delay offsets. The Fourier transformed autocorrelation visibility, $| {\widetilde{V}}_{11}|$, is intrinsically quite peaky in delay space, meaning that a reflection is essentially a narrow spike appearing at its corresponding reflection delay. The reflection amplitude is then estimated via Equation (7), and the reflection phase is estimated by transforming back to frequency space and aligning reflection templates while varying their phase until a squared error metric is minimized. This yields an initial estimate of the reflection parameters, but it needs to be refined in order to suppress the systematic to high dynamic range. In order to refine our parameter estimates, we set up a nonlinear optimization system that perturbs the initial guesses, applies the calibration to the data in frequency space, transforms to Fourier space, estimates the amplitude of the residual reflection bump, and repeats until it is minimized or a stopping threshold is reached. Solving this with an iterative minimization technique allows us to estimate the reflection parameters with sufficient accuracy to suppress the reflection in our systematic simulations by eight orders of magnitude in the power spectrum.

To summarize, our approach for estimating the reflection parameters from the data takes the following steps:

• 1.  
  Zero-pad the autocorrelation in frequency space and apply a windowing function before Fourier transforming to delay space to minimize sidelobe power.
• 2.  
  Fit for the peak of the reflection bump in $| \widetilde{V}|$ via quadratic interpolation of its nearest neighbors.
• 3.  
  The estimated reflection delay, τ, is equal to τ_peak.
• 4.  
  The estimated reflection amplitude, A, is the ratio $| \widetilde{V}(\tau \,={\tau }_{\mathrm{peak}})| /| \widetilde{V}(\tau =0)|$ (Equation (7)).
• 5.  
  Set all modes of $\widetilde{V}$ to zero except the modes nearest τ_peak and Fourier transform back to frequency space to get V_filt.
• 6.  
  The estimated reflection phase, ϕ, is found by minimizing $| {V}_{\mathrm{filt}}-{{Ae}}^{2\pi i\nu \tau +i\phi }{| }^{2}$ while varying ϕ from 0 to 2π.
• 7.  
  Set up a nonlinear optimization that perturbs the initial reflection parameter estimates, applies the calibration and transforms to Fourier space until the residual near the original reflection bump is minimized or a stopping threshold is reached.

We demonstrate this algorithm on a simulated HERA autocorrelation visibility corrupted by a cable reflection with a (frequency-independent) amplitude of 10⁻² at a delay of 600 ns (Figure 3). The natural delay resolution of HERA data is 10 ns, which is much too coarse to achieve precision estimates of the reflection delay. Zero-padding the data by a factor of three gets us to a delay resolution of 3 ns, but this is still not precise enough for accurate reflection delay estimates. By employing quadratic interpolation on the spectral peak, we can recover the input cable delay to roughly ±0.1 ns (left of Figure 3). In this idealized, noise-free simulation, the initial reflection calibration estimates the reflection delay to within ±0.1 ns of its true value, and its phase to within ±0.01 radians, yet we only see systematic power suppression of two orders of magnitude in the visibilities. This is representative of the precision needed to achieve even modest systematic suppression. With the refined reflection calibration, however, we find we can achieve reflection systematic suppression of up to four orders of magnitude in the visibility, and recover the reflection parameters to within 1 part in 10⁶ of their true value.

Zoom In Zoom Out Reset image size

Fiducial EoR levels are expected to be roughly 10⁻⁵ times the peak cross-correlation foreground power in the visibility at ${k}_{\parallel }\sim 0.1\ h$ Mpc⁻¹, and is generally thought to be even weaker at higher k_∥ (Mesinger et al. 2011). If reflection systematics have inherent amplitudes of around 10⁻³, then a few orders of magnitude of further suppression will push them below expected EoR levels at low k. In practice, nonideal effects like frequency evolution of the reflection parameters will limit the precision of reflection calibration, which has been observed in real instruments (Ewall-Wice et al. 2016). Indeed, inclusion of such effects in our systematic simulation will likely degrade our algorithm performance. However, if frequency evolution is a limiting factor for reflection calibration, one simple strategy is to split the full bandwidth into multiple subbands and perform reflection calibration independently on each of them, with the caveat that non-negligible frequency evolution within each subband may need to be mitigated in other ways. One can also do this more self-consistently by estimating the reflection parameters and their frequency dependence across the full band, however, we defer this to future work.

Antenna cross-coupling systematics are a baseline-dependent effect and as such must be modeled and subtracted for each cross-correlation visibility independently. In other works, cross coupling has been modeled as a phase-stable term in the data that can be removed by applying a finite-impulse response (FIR) high-pass filter to the data (Parsons et al. 2010; Ali et al. 2015; M. Kolopanis et al. 2019, in preparation). The algorithm described in this work is conceptually quite similar in that we take advantage of cross-coupling's slow time variability to model it, but is different in its methodology. A comparison of cross-coupling subtraction techniques is done in Appendix B. Note that the method presented here does not assume that the instrumental bandpass has been calibrated out: these two steps are in principle interchangeable. In practice, it helps if at least the antenna cable delays are calibrated out such that the main foreground lobe shows up at the expected delays of τ ≈ 0 ns, but again this is not strictly necessary.

Our semiempirical approach starts with the cross-correlation visibility Fourier transformed across frequency, ${\widetilde{V}}_{12}$, such that it is in the time and delay domain. If we think of the visibility as a 2D rectangular matrix, we can use singular value decomposition (SVD) to decompose the matrix as

$$\begin{eqnarray}&&\widetilde{{\boldsymbol{V}}}={{\boldsymbol{TSD}}}^{\dagger },\end{eqnarray} \tag{ 14 }$$

where ${\boldsymbol{T}}$ is a unitary matrix containing basis vectors (also referred to as eigenmodes) across time, ${\boldsymbol{D}}$ is a unitary matrix containing basis vectors across delay, and ${\boldsymbol{S}}$ is a diagonal matrix containing the weight (or their singular values) of each mode in the data. There may be components of our data matrix, $\widetilde{V}$, that are inherently low rank, like a slowly time-variable systematic for example. Thermal noise in the visibility, on the other hand, occupies the full rank of the data matrix. SVD can help us model and pull out the low-rank components of the matrix, thus providing an approach for systematic removal on a per-baseline basis.

Our SVD-based systematic removal algorithm operating on an individual visibility takes the following steps:

• 1.  
  Fourier transform the visibility waterfall to delay space.
• 2.  
  Apply a rectangular band-stop window across delay to down-weight foregrounds at low delays.
• 3.  
  Decompose the visibility via SVD.
• 4.  
  Choose the first N modes to describe the systematic and truncate the rest.
• 5.  
  Take the outer product of the remaining ${\boldsymbol{T}}$ and ${\boldsymbol{D}}$ modes to form N data-shaped templates.
• 6.  
  Multiply each template with their corresponding singular value in ${\boldsymbol{S}}$ and sum them to generate the full-time and delay-dependent systematic model.
• 7.  
  Fourier transform the systematic model from delay space back to frequency space and subtract it from the data.

A demonstration of this process on simulated visibilities is shown in Figure 4. In this example, the simulation contains foregrounds and cross-coupling systematics, but is free of both thermal noise and EoR components. We start with a simulated visibility spanning roughly 9 hr in LST with 1000 time bins and spanning a bandwidth from 120–180 MHz with 256 frequency bins. By Fourier transforming it to the delay domain (a), we see that foregrounds are confined to low delays, while cross-coupling systematics span a wide range of delays at positive and negative delay offsets. Applying a band-stop windowing function across the delay before taking the SVD (hatched region) assigns zero weight to delay modes dominated by foreground signal at $| \tau | \lt 300$ ns. The result of the SVD shows significant isolation of the information content of the systematic in the visibility into the first eigenmode (c). The first time eigenmode (b) indeed shows it to be slowly time variable, as we would expect for a cross-coupling systematic (see Section 2). The outer product of the first time and delay eigenmodes multiplied with their singular value yields a systematic template with the shape of our original data matrix (d). Taking this single template as our systematic model (equivalent to setting N = 1) and subtracting it from the data yields the systematic-subtracted data (f).

Zoom In Zoom Out Reset image size

We see in (d) that a cross-coupling model with N = 1 provides good subtraction of the systematic, but is not enough to completely remove it: based on Figure 4 (c) we can see that it provides effectively three orders of magnitude of suppression, which, depending on the inherent amplitude of the systematic, may or may not be enough to push it below fiducial EoR levels. We can remove more and more of the systematic by increasing the number of SVD modes we incorporate into the systematic model. However, as is the case with empirically based models, this has the side effect of possibly introducing structure from other components of the visibility that we do not want in our systematic model, such as the EoR itself. If the EoR signal was somehow soaked up by our systematic model, then by subtracting the model from the data we are inducing EoR signal loss, which is highly undesirable.

To limit EoR signal loss in the process of systematic subtraction, we can filter the systematic model to reject Fourier modes that we know hold EoR power. In general, if a signal occupies the visibility in the fringe-rate domain with variance given by σ(f)² and we enact a filter on it by multiplying by a weighting function w(f), then the total power of the signal before filtering, P_before, can be related to the total power after filtering, P_after, as

$$\begin{eqnarray}&&\displaystyle \frac{{P}_{\mathrm{before}}}{{P}_{\mathrm{after}}}=\displaystyle \frac{\int {df}\ \sigma {(f)}^{2}}{\int {df}\ \sigma {(f)}^{2}w{(f)}^{2}}.\end{eqnarray} \tag{ 15 }$$

Therefore, if we know statistically how the EoR will populate the visibilities in the fringe-rate domain—in a sense deriving their power spectral density functions (PSDs)—we can construct a Fourier filter that is tailored to reject Fourier modes in the systematic model that we know hold EoR power.

This is closely related to the optimal fringe-rate formalism outlined in Parsons et al. (2016). In Figure 5 we show peak-normalized PSDs of an EoR sky model at 120 MHz for various baselines in the array, which describe the relative amount of signal power occupied by different fringe rates. We derive these via ensemble simulations of the same EoR sky while varying the initial seed (Appendix A). The PSD of a simulated HERA cross-coupling systematic is also plotted (black-dashed curve). While centered at a fringe rate of 0 mHz, the systematic has a tail that extends out to negative and positive fringe rates, the latter being where most of the power from the EoR also lies. Conversely, while most of the EoR power lies at positive fringe rates for many baselines, some of its power also extends to zero and negative fringe rates. Baselines that are longer along the east–west direction have more EoR power pushed to higher fringe rates and thus are more naturally separated from cross-coupling systematics, while baselines that are purely north–south in orientation see an EoR signal that is centered at f ∼ 0 mHz, and almost completely overlaps the cross-coupling systematic.

Zoom In Zoom Out Reset image size

Each curve in Figure 5 integrated out to 99% of their total area yields the domain in fringe-rate space where 99% of the EoR power is contained in the visibility. We tabulate these bounds in Table 1 for a few HERA baselines at ν = 120 MHz. If we low-pass filter any of the visibilities in Figure 5 in fringe-rate space by applying a symmetric top-hat filter with a maximum extent f_max given by these lower bounds, then Equation (15) tells us we will retain 99% of the EoR power in the data after filtering, which for our purposes is an effectively lossless operation given other more dominant sources of error. This result means that our ability to safely remove cross-coupling systematics is baseline-dependent: for baselines with large east–west lengths (e.g., blue curve), we can filter out the vast majority of the systematic without attenuating the EoR. For shorter baselines (e.g., green curve), we may only be able to remove part of the systematic, and for baselines oriented along the north–south direction (e.g., red), we may not be able to remove any cross-coupling systematics, if they exist.

Armed with the ability to filter an arbitrary signal without attenuating its EoR component, we can return to the problem of choosing the appropriate number of eigenmodes to use in describing a cross-coupling systematic in the data. We noted in Figure 4 that by increasing the number of SVD eigenmodes used to describe our systematic model, we might be able to remove more of the systematic from the data. Figure 6 proves this on a simulated visibility, now simulated with an EoR and foreground component (green curve) corrupted by a cross-coupling systematic at high delays (blue curve). We also show SVD-based systematic removal with increasingly more eigenmodes used to describe the systematic (orange curve). We can see that going from N = 1 (a) to N = 5 (b) enables us to subtract more of the systematic from the data. However, at some point we expect low-level eigenmodes to be influenced by other components of the data, such as EoR, foregrounds, or noise, which raises the possibility of removing those components along with the systematic. This is shown in (c), where with N = 15 we have over subtracted the systematic and caused signal loss of EoR power at high delays. One might conclude from this that N = 5 is the "sweet spot" choice for the number of eigenmodes to use, but this choice is conditional on the relative amplitude between EoR and systematic: without knowing the amplitude of EoR in the data a priori, we have no way of knowing the appropriate number of eigenmodes to use that would enable us to subtract the systematic without attenuating EoR. This makes the algorithm as originally described in effect unusable, because it is an operation that is dangerously lossy to EoR signal.

Zoom In Zoom Out Reset image size

The solution to this problem is to apply a low-pass time filter to the systematic model that is tailored to reject fringe-rate modes occupied by the EoR. Specifically, we can apply a filter to the systematic model ${\boldsymbol{T}}$ matrix that only keeps structure below some predefined maximum fringe rate, f_max, such that in the process of subtracting it from the data, all fringe rates $| f| \gt | {f}_{\max }$ are left unaffected. For example, if we could tolerate a maximum of 1% attenuation of EoR power in the process of systematic removal, then f_max would be set at the lower bounds tabulated in Table 1. The result of applying such a filter to the SVD eigenmodes is demonstrated in Figure 6(d), which shows the systematic-subtracted data with N = 15 having first applied a low-pass filter to ${\boldsymbol{T}}$. Although a significant amount of systematic remains, we can now be confident that we have not attenuated the EoR signal in the data, even while using a large number of eigenmodes to describe the systematic.

In this section we have argued and shown via sky and instrument simulations that we can construct a cross-coupling model that removes the vast majority of the systematic while remaining lossless to the EoR signal (for certain baseline orientations). In real data, however, the fidelity of this model will be fundamentally limited by the thermal noise floor of the observation, as is the case for any signal term modeled on a per-baseline basis. If the cross-coupling systematic is truly baseline-dependent and is uncorrelated between baselines, then the residual systematic term will integrate down like thermal noise when we combine visibilities and we would not expect it to re-appear in the integrated power spectra. In HERA, for example, there is evidence that the observed cross-coupling systematics are at least partially uncorrelated between baselines (Kern et al. 2019).

To test signal loss in the context of reflection calibration, we precompute the visibilities for a single diffuse foreground model and 100 independent EoR models, with each simulation spanning 8 hr of LST and 1000 individual time integrations. Because HERA has a beam crossing time of about 1 hr, this yields an effective number of independent foreground + EoR simulations of ∼800. The adopted EoR model is an uncorrelated Gaussian field across angular position and frequency with a variance of 25 mK² (Appendix A).

A single signal loss trial takes the following steps. First we choose a random EoR model from our library of precomputed visibilities and add it to our foreground visibility (V₁). We then make a copy of it and insert a reflection with a delay of 600 ns, a random phase, and a single amplitude across frequency using Equation (4) (V₂). We then model the reflection in the simulated autocorrelation knowing only the approximate delay range at which it appears, and then apply the derived gain solution to the cross-correlation visibilities (V₃). Next we compute power spectra of each data product (P₁, P₂, and P₃). We then repeat this on the order of 100 times, each with a different random EoR model, and then take their average to approximate the ensemble average in Equation (19). We then form the R₂ and R₃ metrics as a function of time and delay, i.e., ${R}_{3}(t,\tau )$, and average over time to collapse them onto a single axis across delay. This entire procedure produces one signal loss trial, which is defined uniquely by the amplitude, A, of the reflection inserted into the visibility.

Figure 7 shows multiple trials for different reflection amplitudes in the range of 10⁻⁶–10⁻¹. The top row shows power spectra of each of the three visibility products as a function of delay for one trial when A = 10^−2.5. Recall that the reflection amplitude is defined with respect to the visibility, meaning that the observed reflection amplitude in the power spectrum is A². The bottom row shows a heatmap of the signal loss R metric as a function of delay (x-axis) and each trial's reflection amplitude (y-axis). The left panels shows R₂ and the right panel shows R₃, highlighting the amount of delay space that is brought down to R ∼ 1 after reflection calibration, with negligible amounts of signal loss (purple-shaded regions). Furthermore, the weak levels of fluctuating residual systematic and signal loss observed at the ∼2% level are within the $1/\sqrt{N}$ sample variance of our finite ensemble average. For context, HERA cable reflection amplitudes are seen at around 10⁻³ (Kern et al. 2019).

Zoom In Zoom Out Reset image size

4.1.1. Multi-reflection Regime

Above we probed for signal loss when performing reflection calibration on a single reflection that was isolated in delay space. Next, we relax this assumption and test how the the algorithm performance and signal loss properties change when we add in more reflections, which is relevant for any instrument with multiple cables, or with cables that have sub-reflections along the length of the cable (Ewall-Wice et al. 2016; Kern et al. 2019). We choose to model the relative amplitude of these reflections as an inverse power law as a function of delay with a nominal reflection amplitude of $A=3\times {10}^{3}$. Our algorithm models and calibrates out each reflection one at a time, starting with the reflection with the largest amplitude. We do not feed the algorithm the position of each reflection (it searches for it automatically within a specified range of delays), but we do assume we know the number of reflection inherent in the data, which controls how many times we iterate the algorithm.

Our first test shown in Figure 8 involves only five reflections inserted across a relatively wide region in delay, such that they can be considered nonoverlapping (or unconfused). In the top panels, we show power spectra of the uncorrupted data (P₃; green), the corrupted data (P₂; blue line) and the reflection calibrated data (P₃; orange-dashed line), along with a line marking the EoR amplitude in the data (gray). We show this for the autocorrelation visibility (left) and a 29 m cross-correlation visibility (right). In the bottom panels we show the signal loss metrics R₂ (blue) and R₃ (orange-dashed line). Recall that reflection calibration builds up a set of gains strictly from the autocorrelation, and then applies those gains to the cross correlations. We find that our algorithm performs exceptionally well in this regime: removing reflections down to the inherent sidelobe floor of the autocorrelation, which is more than enough to bring the systematics to the EoR level in the cross-correlation visibility. We find that while reflection calibration can lead to slight signal loss in the autocorrelation visibility, we find no appreciable levels of signal loss in the cross correlation.

Zoom In Zoom Out Reset image size

Our next test shown in Figure 9 increases the number of reflection modes inserted into the same region, such that they become almost entirely overlapping. In this case we can see that our algorithm fails to perfectly calibrate out the reflections in both the autocorrelation (left) and cross correlation (right) due to the partial confusion. Nonetheless, we still find that while imperfect reflection calibration can lead to slight signal loss in the autocorrelation, the cross correlation is resistant to signal loss. Reflection calibration's resistance to signal loss in the cross correlations is perhaps not surprising, given that reflection calibration operates in an antenna-based space, while EoR and other sky signals live in a baseline-based space. Furthermore, our algorithm solely uses the autocorrelations to estimate the reflection parameters.

Zoom In Zoom Out Reset image size

In total, our findings suggest that (1) our reflection calibration algorithm performs moderately well even in the many-reflection regime, and that (2) even if a broad delay region is contaminated by reflections, we can still in principle use this region for EoR measurements (or upper limits in the case of imperfect systematic removal) after reflection calibration because it does not suffer appreciable levels of signal loss.

In this section, we quantify signal loss for cross-coupling subtraction in a similar manner. Based on our conclusions from Section 3.2, we expect our cross-coupling subtraction to be effectively lossless to EoR by construction, but testing this against ensemble signal loss trials is a good double-check and validation of our arguments. The rough functional form of the simulated cross coupling inserted into the visibilities is informed by cross-coupling systematics observed in the HERA Phase I system (Kern et al. 2019). The simulations used are fundamentally the same as those used in Section 4.1, except simulated with cross-coupling systematics rather than cable reflections. To simulate cross coupling in a cross-correlation visibility, we use Equation (12) to insert ∼25 modes spanning $700\lt | \tau | \lt 900$ ns with a decaying power law as a function of $| \tau |$ for their relative amplitudes, normalized such that the maximum amplitude relative to the peak autocorrelation foreground power equals a predefined amplitude, A. In the process of systematic removal, we model the systematic with N = 20 SVD eigenmodes and apply a low-pass fringe-rate filter on the SVD ${\boldsymbol{T}}$ modes using a Gaussian process (GP) smoothing (see Appendix B) with a maximum fringe rate given by the lower bound in Table 1. We then Fourier transform the data from frequency to delay space using a seven-term Blackman–Harris window, and average across ensemble trials. After forming the R₂ and R₃ metrics as a function of time and delay, we truncate 5% of the time bins on either edge of the time axis before taking their time average to limit the influence of boundary effects in the smoothing process.

The result is shown in Figures 10 and 11, which shows signal loss trials run on a 15 m east–west baseline and a 29 m east–west baseline, respectively. As expected, we see better systematic suppression for the longer baseline, where the EoR signal is inherently more isolated from the systematic in fringe-rate space. We also see that in the case of strong coupling amplitudes we cannot completely suppress the systematic down to EoR levels; however, we can nonetheless suppress the systematic by roughly four orders of magnitude in power for the 15 m baseline and eight orders of magnitude in power for the 29 m baseline. Importantly, we show that the algorithm presented does not significantly attenuate EoR in the data, with residual fluctuations about R = 1 at the 1% level in power. While this result is not surprising given the fact that in Section 3 we constructed our systematic model to explicitly be lossless to EoR, it is still a useful cross-check on our algorithm and its implementation on the data.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We use healvis-simulated HERA observations to calculate the expected PSD functions of an EoR-like signal for HERA visibilities. At any given time, a point source locked to the celestial sphere generates a complex sinusoid in the inteferometric visibility with a time period determined by its position on the sky and the fringe profile of the baseline at hand. Over a short time interval, its Fourier transform across time is nearly a delta function at a fringe rate set by the inverse of its time period. Thus, over short time intervals, points on the sky map to specific fringe rates in the visibility (Parsons et al. 2016), which is also related to the m-mode analysis (Shaw et al. 2014). The analytically derived "optimal fringe-rate filter" in Parsons et al. (2016) is a filter that minimizes the noise component of the visibility-based power spectrum error bars, and is related to the PSD of EoR-like signals in the visibility. However, the analytic derivation hinges crucially on the assumption that the beam crossing time is much longer than the fringe-crossing time for a point source on the sky, which was a fairly valid assumption for a wide-field experiment like PAPER. This is not the case for HERA, which has both shorter baselines with wider fringes and also a more compact primary beam. Therefore, calculating the PSD of an EoR-like sky signal is more easily done numerically. We do this by generating over 100 independent EoR sky models with the same variance but different initial seeds. We perform HERA visibility simulations of each sky using healvis, Fourier transform them from LST into the fringe-rate domain, and then take their absolute value and average across each realization. The square of these profiles is shown in Figure 5, which represents a numerically derived PSD of a theoretical EoR-like signal in HERA visibilities. These profiles allow us to tailor visibility Fourier filters to do things like minimize EoR signal loss in systematic subtraction, or to minimize the thermal noise component in the power spectrum error bars. We tabulate the 99% power bounds of these curves in Table 1 for the few baselines presented in this study. We also plot these bounds and provide a fitting formula as a function of projected east–west baseline length in Figure 14, such that one can extrapolate these data points to other baselines in the HERA array. Error bars on the 99% bounds in Figure 14 are calculated via bootstrap resampling over the independent realizations (Efron & Tibshirani 1994).

Zoom In Zoom Out Reset image size