© The Authors 2025

1. Introduction

One of the most comprehensive probes of the astrophysics and cosmology of the early Universe is the redshifted 21 cm signal from neutral hydrogen (HI) (Madau et al. 1997; Shaver et al. 1999). The fluctuations in the brightness temperature distribution of the 21 cm signal against the cosmic microwave background radiation at high redshifts (e.g., z > 6) can be probed by low-frequency radio interferometers. However, none of the currently operating instruments have the thermal noise sensitivity to detect the faint signal in tomographic images in reasonable observing time (i.e., thousands of hours). Hence, they attempt to detect the 21 cm signal statistically and quantify it in terms of a power spectrum. Several telescopes such as the GMRT¹ (Paciga et al. 2013), MWA² (Ewall-Wice et al. 2016; Barry et al. 2019; Li et al. 2019; Trott et al. 2020; Yoshiura et al. 2021), LOFAR³ (Patil et al. 2017; Gehlot et al. 2019, 2020; Mertens et al. 2020), PAPER⁴ (Kolopanis et al. 2019), HERA⁵ (Abdurashidova et al. 2022; Adams et al. 2023), OVRO-LWA⁶ (Eastwood et al. 2019; Garsden et al. 2021), and NenuFAR⁷ (Munshi et al. 2024) have already set upper limits on the 21 cm signal power spectrum at cosmic dawn (CD) and during the subsequent epoch of reionization (EoR). The future SKA⁸ (Dewdney et al. 2009; Koopmans et al. 2015) is expected to have the sensitivity to generate tomographic maps of the 21 cm signal in the early Universe. A detection or strong upper limits allow us to constrain the underlying astrophysical models that cause the fluctuations of the 21 cm signal during these epochs (e.g., Furlanetto et al. 2006; Pritchard & Loeb 2012). However, the orders of magnitude brighter foregrounds and instrumental systematics are the two most important obstacles to the detection of the power spectrum of the 21 cm signal.

A crucial diagnostic in any 21 cm cosmology analysis is the cylindrical power spectrum in which the signal brightness temperature cube is averaged within cylindrical shells in the Fourier space, with the axis of the cylinders pointing along the line of sight (LoS). The cylindrical power spectrum is an estimator of the variance in the observed data as a function of the spatial mode along the LoS (k_∥) and the spatial mode on the plane perpendicular to the LoS (k_⊥). The intrinsic chromatic nature of an interferometer causes the spatial scale probed by a certain baseline to vary with frequency. This leads to a mixture of spatial modes with spectral modes, a phenomenon commonly known as mode mixing (Liu & Tegmark 2011). The instrumental mode mixing combined with the smooth spectral nature of foregrounds results in them occupying a region in the cylindrical power spectrum known as the foreground wedge (Datta et al. 2010; Vedantham et al. 2012; Morales et al. 2012; Liu et al. 2014a,b; Murray & Trott 2018). The two main approaches used in the mitigation of foregrounds are subtraction, where the smooth spectral nature of the foregrounds is used to model and subtract them from the visibilities, and avoidance, where the power spectrum is constructed only from the k modes that fall outside the foreground wedge and are thus assumed to be largely free of foregrounds.

The maximum extent of the foregrounds in the cylindrical power spectrum is described by a straight line in the (k_⊥, k_∥) space, commonly known as the horizon line (Vedantham et al. 2012; Morales et al. 2012; Trott et al. 2012; Thyagarajan et al. 2013). However, for cylindrical power spectra estimated by phase-tracking instruments, as opposed to drift-scan approaches, the equations describing such a line are derived under the assumption that sky sources lie on the tangent plane of the sky at the phase center. This assumption leads to the complex visibility function sampled by an interferometer for an arbitrary phase center taking the form of a two-dimensional Fourier transform. As a result, the conventional horizon line equations are only valid for instruments with a small field of view (FoV) without primary beam sidelobes, where such a flat-sky approximation holds. For wide-field instruments, sky sources far away from the phase center can strongly impact the power spectrum, which these equations do not describe well.

In this paper we derive the horizon line equation in the cylindrical power spectrum for a pointing observation without taking a flat-sky approximation. These equations describe the actual maximum extent of foregrounds in cylindrical power spectra estimated using wide-field instruments. This is important because it helps optimize our instruments, observations, and processing for signal extraction. It is also crucial for diagnoses, by helping trace the source of power observed in the power spectrum to the emitter. We find that the predictions of these equations differ significantly from those of the conventional horizon line equations. In addition to the horizon line equation, we also derive an equation that describes the extent of contamination due to a specific source in the sky on the cylindrical power spectrum, also without a flat-sky approximation. Again, we find that the equation deviates from what we obtain from an equation derived using the flat-sky approximation. We confirm the validity of both sets of equations through realistic simulations for a range of observation parameters. We present the Python library pslines, which can be used to generate and plot the horizon and source lines derived in this analysis.

The paper is organized as follows. Section 2 introduces the problem and illustrates the impact of pointing observations on signal delays. In Sect. 3 we derive the analytical expressions for the maximum delays in the visibilities. In Sect. 4 we use these expressions to derive equations for the horizon line in the cylindrical power spectrum without imposing a flat-sky approximation. In Sect. 5 we derive an expression for the modes in the cylindrical power spectrum most strongly affected by a specific point source. Section 6 describes the simulations performed to validate the horizon and source line equations. In Sect. 7 we discuss the limitations of the derived equations and their potential applications in 21 cm cosmology. Section 8 summarizes the main results of this paper.

2. Pointing observations and delays

Power spectrum analyses performed in 21 cm cosmology can be divided into two major groups: the delay spectrum approach and the reconstructed power spectrum approach (Morales et al. 2019). In the delay spectrum approach (Parsons et al. 2012), the individual visibilities from pairs of interferometric elements are directly Fourier transformed along the frequency axis to obtain the visibilities as a function of the signal delay (η) between the pair of elements. The corresponding power spectrum is closely related to the cosmological power spectrum and contains all statistical information about the 21 cm field under the assumption of Gaussianity. This method is usually employed for drift scan telescopes where the data is not phased to a direction fixed to the sky. In this case, the maximum delay (|η|^max) is given by

$$\begin{array}{c}\hfill {|\eta |}^{\max }=|\mathit{b}|/c,\end{array}$$(1)

where b is a baseline vector in units of meters, with the maximum delay corresponding to a source at the physical horizon. Therefore Eq. (1) represents the horizon line equation for the delay power spectrum for a drift scan observation.

In the reconstructed or sky-plane power spectrum approach, the power spectrum is constructed from the gridded visibility data in the uvν space. This requires choosing a phase center fixed to a direction in the sky which then moves with respect to the zenith during an observation. The raw cross-correlation data measured by the instrument needs to be phased in the direction of the phase center. This is commonly referred to as geometric delay correction, and this results in the phase center having a phase of zero across all stations. We consider a Cartesian coordinate system with axes ABC, where the C axis points along the direction of the phase center. Let $\widehat{\mathit{s}}=(l,m,n)$ be a unit vector in the direction of a particular source in the sky and b_p = (U_p, V_p, W_p) be the position vector of the p-th interferometric element. The electric field due to a unit-flux unpolarized source at the p-th element and frequency ν is

$$\begin{array}{c}\hfill e_p(\nu )=e^{-2\pi i({\mathit{b}}_p·\widehat{\mathit{s}})\nu /c}.\end{array}$$

The cross-correlated raw data from the p-th and q-th elements are

$$\begin{array}{cc}\hfill V_{\mathit{pq}}(\nu )=e_p(\nu )e_q^{\ast }(\nu )& =e^{-2\pi i({\mathit{b}}_p-{\mathit{b}}_q)·\widehat{\mathit{s}}\nu /c}\hfill \\ \hfill & =e^{-2\pi i(U_{\mathit{pq}}l+V_{\mathit{pq}}m+W_{\mathit{pq}}n)\nu /c},\hfill \end{array}$$

where b_pq = b_p − b_q is the baseline vector, having the components (U_pq, V_pq, W_pq) in this coordinate system. The phase center $\widehat{\mathit{p}}$ has the components (0, 0, 1) and the phase in the direction of the phase center is e^{−2πiW_pqν/c}. Phasing the data amounts to applying the inverse of this phase to the data, so that the phase center has a phase of zero for all baselines. After phasing, the visibility data is

$$\begin{array}{c}\hfill V_{\mathit{pq}}(\nu )=e^{-2\pi i({\mathit{b}}_{\mathit{pq}}·(\widehat{\mathit{s}}-\widehat{\mathit{p}}))\nu /c}=e^{-2\pi i(U_{\mathit{pq}}l+V_{\mathit{pq}}m+W_{\mathit{pq}}(n-1))\nu /c}.\end{array}$$(2)

This equation gives the visibility at a baseline b_pq due to a unit-flux source $\widehat{\mathit{s}}$ for a phase center $\widehat{\mathit{p}}$, assuming noiseless isotropic receivers with unit gains.

We now examine how phasing the data to a certain phase center affects the horizon line in the delay power spectrum. We hereby drop the pq suffix. Following Eq. (2) for a unit-flux unpolarized source at $\widehat{\mathit{s}}$, the delay spectrum is given by

$$\begin{array}{cc}\hfill V(\eta )& ={\displaystyle {\int }_{\nu }}{e}^{-2\pi i(Ul+Vm+W(n-1))\nu /c}{e}^{2\pi i\eta \nu }\mathrm{d}\nu \hfill \\ \hfill & =\delta (\eta -{\eta }_0),\hfill \\ \hfill \mathrm{where}{\eta }_0& =\mathit{b}·(\widehat{\mathit{s}}-\widehat{\mathit{p}})/c=(Ul+Vm+W(n-1))/c.\hfill \end{array}$$(3)

Here δ(η − η₀) is a Dirac delta function⁹ with a peak at η = η₀. Because $\widehat{\mathit{s}}$ is a unit vector,

$$\begin{array}{cc}& -\frac{|\mathit{b}|}{c}\le \frac{\mathit{b}·\widehat{\mathit{s}}}{c}\le \frac{|\mathit{b}|}{c}\hfill \\ \hfill \mathrm{or},& \frac{-|\mathit{b}|-W}{c}\le \frac{Ul+Vm+W(n-1)}{c}\le \frac{|\mathit{b}|-W}{c}\hfill \\ \hfill \mathrm{or},& \frac{-|\mathit{b}|-W}{c}\le {\eta }_0\le \frac{|\mathit{b}|-W}{c}.\hfill \end{array}$$(4)

Thus, the maximum (horizon) delay is not simply |b|/c. It depends on the W coordinate for the particular baseline. The range of delays remains the same but gets shifted through geometric phasing due to the artificial injection of an additional delay.

To illustrate the impact of phasing the data on the delays, we constructed the delay power spectrum for visibilities from a 1 h full-sky simulation of a NenuFAR observation of the north celestial pole (NCP). The details of the sky model and simulation parameters are described in Sect. 6. In Fig. 1, we show the delay power spectra before and after phasing the data to the NCP. While plotting, the delay spectra of all baselines with a negative W component have been flipped along the delay axis. This effectively makes sure that all baselines have a positive component toward the north (i.e., the vector joining the p-th element to the q-th element has a positive component toward the north). As a result, before phasing, positive delays correspond to sources toward the north of the zenith and negative delays correspond to sources toward the south of the zenith. The dotted lines correspond to Eq. (1), while the solid lines indicate the horizon line from Eq. (4). We see that the two equations correctly describe the horizon limit, before and after phasing the data. There is a slight spillover due to spectral leakage coming from the Blackman-Harris filter which was applied to the data along frequency before the Fourier transformation. This shift in the delays due to phasing can be intuitively understood as a shift of the coordinate system such that the phase center has a delay of zero. Since the NCP is toward the north of the zenith, phasing the data to the NCP decreases the horizon toward the north of the phase center, and increases the horizon limit toward its south. However, it should be noted that this is not a simple rotation of the coordinate system, and the delay due to sources π/2 away from the phase center is no longer confined to |b|/c after phasing (discussed in Appendix A). Normally we plot only the absolute values of the delay, and in that case, the negative delays will spill over the traditional horizon line.

Fig. 1. Effect of geometric phasing on the delay power spectrum. The two panels show the delay power spectra before (left) and after (right) phasing to the NCP for a full-sky simulation with NenuFAR. The dotted lines correspond to the standard horizon limit of Eq. (1), and the solid lines correspond to the extended horizon due to phasing from Eq. (4).

It should be noted that this plot is possible for a long observation (not a snapshot) only when the phase center is the NCP or the south celestial pole (SCP) where the W coordinates are independent of time. Additionally, the lines can be plotted as long as the data is not binned in the baseline direction, as the W values are dependent on the baseline directions. This simple example illustrates how the classic horizon delay lines are modified for interferometers that phase-track a fixed point in the sky. In the next section, we derive a more complete description of maximum delays for phase tracking instruments.

3. Maximum delays

Evaluating the maximum delay from Eq. (4) requires the knowledge of the individual W coordinates for the baseline under consideration. Thus, Eq. (4) cannot describe a situation where multiple baselines oriented in different directions are averaged together, as is the case for the cylindrical power spectrum. Additionally, for any phase center other than the NCP and SCP, the W coordinates, and consequently the horizon lines, will shift with time. In an observation, we typically have a large number of baselines oriented in different directions and a large number of sources distributed throughout the sky. In this section, we derive general equations for the horizon delay as a function of baseline length similar to Eq. (1), but accounting for the phasing to the phase center.

3.1. Analytical treatment

Here, we assume a three-dimensional coordinate system XYZ fixed to the Earth with the Z axis pointing toward the NCP, the X axis pointing due east, and the Y axis pointing in the direction Z × X. Consider another three-dimensional coordinate system ENU (east-north-up) where U points to the zenith, N points to the north along the plane perpendicular to the zenith, and E points to the east. The two coordinate systems are illustrated in Fig. 2. The parameters H, δ and H₀, δ₀ are the hour angle, declination (Dec) of the source ($\widehat{\mathit{s}}$) and the phase center ($\widehat{\mathit{p}}$), respectively, and ϕ is the latitude of the location of the telescope. Table 1 summarizes the notation used in the paper. The coordinates of the baseline (b), source ($\widehat{\mathit{s}}$), zenith ($\widehat{\mathit{z}}$), and the phase center ($\widehat{\mathit{p}}$) in these systems are

Fig. 2. Coordinate systems used in the calculations. The origin of the coordinate system is placed at the location of the array for illustration. Since the sources are at an infinite distance, the exact location of the coordinate center is irrelevant.

$$\begin{array}{cc}\hfill \mathit{b}& =(B_E,B_N,B_U)[\mathrm{ENU}]\hfill \\ \hfill & =(B_E,B_{N}sin\varphi -B_{U}cos\varphi ,B_{N}cos\varphi +B_{U}sin\varphi )[\mathrm{XYZ}],\hfill \\ \hfill \widehat{\mathit{s}}& =(-cos\delta sinH,-cos\delta cosH,sin\delta )[\mathrm{XYZ}],\hfill \\ \hfill \widehat{\mathit{z}}& =(0,-cos\varphi ,sin\varphi ){\text{[XYZ]}}^{10},\hfill \\ \hfill \widehat{\mathit{p}}& =(-cos{\delta }_{0}sin{H}_0,-cos{\delta }_{0}cos{H}_0,sin{\delta }_0)\text{[XYZ]}.\hfill \end{array}$$

We now assume that¹⁰ all baselines lie on the plane perpendicular to the zenith (or the Earth’s tangent plane, considering the Earth to be a sphere), so B_U = 0. This is a fair assumption for EoR or CD experiments where the baselines used to construct the power spectrum are small relative to the curvature of the Earth. The components can now be written in terms of the baseline length |b| and the angle θ it makes with the E axis. In this analysis we assume that the array has baselines in all directions and θ can take any value between 0 and 2π. This is a fair assumption for imaging interferometers with dense cores used in 21 cm cosmology. The coordinates of the baseline are then given by

$$\begin{array}{c}\hfill \mathit{b}=|\mathit{b}|(cos\theta ,sin\theta sin\varphi ,sin\theta cos\varphi )[\mathrm{XYZ}].\end{array}$$

The corresponding expression for the delay due to a source (η₀ from Eq. (3)) is given by

$$\begin{array}{cc}\hfill {\eta }_0=& \frac{|\mathit{b}|}{c}(cos\theta cos{\delta }_{0}sin{H}_0-cos\theta cos\delta sinH\hfill \\ \hfill & +sin\theta sin\varphi cos{\delta }_{0}cos{H}_0-sin\theta sin\varphi cos\delta cosH\hfill \\ \hfill & +sin\theta cos\varphi sin\delta -sin\theta cos\varphi sin{\delta }_0).\hfill \end{array}$$(5)

The condition that both $\widehat{\mathit{s}}$ and $\widehat{\mathit{p}}$ are above the physical horizon takes the form

$$\begin{array}{cc}\hfill \widehat{\mathit{s}}·\widehat{\mathit{z}}\ge 0& \mathrm{or},cos\delta cosHcos\varphi +sin\delta sin\varphi \ge 0;\hfill \\ \hfill \widehat{\mathit{p}}·\widehat{\mathit{z}}\ge 0& \mathrm{or},cos{\delta }_{0}cos{H}_{0}cos\varphi +sin{\delta }_{0}sin\varphi \ge 0.\hfill \end{array}$$(6)

We would like to obtain the extrema of Eq. (5) for a given δ₀, subject to the horizon condition (Eq. (6)). Instead of solving the full three-dimensional equation ∇η₀(H, δ, H₀, θ) = 0, we first reduce the problem to the plane containing the zenith vector and the Earth’s axis (the meridian plane). It should be noted that for any phase center other than the NCP and SCP, this reduction of the problem cannot be done without loss of generality, since H₀ affects the delay values and the extrema do not necessarily lie on the meridian plane. We show in Sect. 3.2 that this reduction is valid only in the regime of |δ₀ + ϕ|> π/2. The situation for the NCP is described in Appendix A. Nevertheless, to gain insight into the problem, we investigate the situation of a north-south baseline with the phase center and source both lying on the meridian plane. To get an extreme case of maximum delay, we choose H₀ and H so that the phase center and the source can be located the furthest apart. This amounts to setting H = 0, H₀ = π, and θ = π/2 in Eq. (5). The equation simplifies to

$$\begin{array}{cc}\hfill {\eta }_0(\varphi ,\delta ,{\delta }_0)=\frac{|\mathit{b}|}{c}& (-sin\varphi cos{\delta }_0-sin\varphi cos\delta +cos\varphi sin\delta \hfill \\ \hfill & -cos\varphi sin{\delta }_0).\hfill \end{array}$$

The declination of the source responsible for the maximum absolute delay is analytically obtained by setting $\frac{\partial {\eta }_0}{\partial \delta }=0$. This gives us the following:

$$\begin{array}{cc}\hfill \delta {|}_{\max }& =\{\begin{array}{cc}\varphi -\pi /2\hfill & \mathrm{if}\varphi \ge 0\hfill \\ \varphi +\pi /2\hfill & \mathrm{if}\varphi \le 0\hfill \end{array}\hfill \\ \hfill \mathrm{and}|{\eta }_0{|}^{\max }& =\{\begin{array}{cc}|\mathit{b}|(1+sin({\delta }_0+\varphi ))/c\hfill & \mathrm{if}\varphi \ge 0\hfill \\ |\mathit{b}|(1-sin({\delta }_0+\varphi ))/c\hfill & \mathrm{if}\varphi \le 0.\hfill \end{array}\hfill \end{array}$$(7)

The maximum delays occur due to a source at the physical horizon¹¹ and hence it is not necessary to perform the maximization subject to Eq. (6). These equations reduce to the ones obtained for NCP phasing at δ₀ = π/2 and ϕ ≥ 0 (Eq. (A.3) in Appendix A). In the next section, we show that this equation is only valid for a subset of all possible situations since the assumption that the maximum delays will occur when the source and phase center lie on the meridian plane is not always correct.

3.2. Full synthesis

We consider a full synthesis when H₀ takes all values between 0 and 2π (i.e., a 24 h synthesis) during the observation due to the rotation of the Earth. At each such time, there can be a source anywhere in the sky and any possible baseline can be responsible for the maximum delay. The expression we want to maximize is $\mathit{b}·(\widehat{\mathit{s}}-\widehat{\mathit{p}})/c$. Now $|\widehat{\mathit{s}}-\widehat{\mathit{p}}|$ is maximum when $\widehat{\mathit{s}}$ and $\widehat{\mathit{p}}$ point in the opposite directions, and in that case $(\widehat{\mathit{s}}-\widehat{\mathit{p}})$ has a magnitude of 2. Additionally, $\mathit{b}·(\widehat{\mathit{s}}-\widehat{\mathit{p}})$ is maximum when b is in the same direction as $(\widehat{\mathit{s}}-\widehat{\mathit{p}})$. Let us only consider the situation when the phase center is above the horizon for at least a part of the full synthesis, where |δ₀ − ϕ|≤π/2. Now we consider two cases:

• When |δ₀ + ϕ|≤π/2 (i.e., when the phase center is non-circumpolar), at some point during a full synthesis, $\widehat{\mathit{p}}$ will go below the horizon. Whenever H₀ assumes a value such that the phase center crosses the horizon, H and δ get optimized to make $\widehat{\mathit{s}}$ point along the opposite direction of $\widehat{\mathit{p}}$. Finally θ is optimized to have b pointing along $\widehat{\mathit{p}}$ or $-\widehat{\mathit{p}}$. Since the baselines are constrained to lie on the Earth’s tangent plane, this situation can occur when $\widehat{\mathit{s}}$, $\widehat{\mathit{p}}$ and b all lie on the physical horizon, irrespective of whether we impose the horizon condition of Eq. (6). Thus, the maximum delay is always 2|b|/c.

• When |δ₀ + ϕ|> π/2 (i.e., when the phase center is circumpolar), $\widehat{\mathit{p}}$ will always be above the horizon and can never cross it. Here the maximum delay occurs when $\widehat{\mathit{p}}$ is closest to the horizon, which by definition lies on the meridian plane. Thus, H₀ = π and the assumption taken to derive Eq. (7) is satisfied and Eq. (7) is valid in this regime.

In order to account for both regimes of |δ₀ + ϕ|≤π/2 and |δ₀ + ϕ|> π/2, Eq. (7) is modified to have the form:

$$\begin{array}{c}\hfill |{\eta }_0{|}^{\max }=\frac{|\mathit{b}|}{c}(1+sin(max\{|\varphi +{\delta }_0|,\pi /2\}))\mathrm{for}|{\delta }_0-\varphi |\le \pi /2.\end{array}$$(8)

We verify Eq. (8) by numerically maximizing Eq. (5) over H, δ, H₀, and θ, subject to Eq. (6) for a given ϕ and δ₀. This was done for a baseline length of 100 m, by evaluating Eq. (5) on a discrete grid in H, δ, H₀, θ space, masking out regions where Eq. (6) is not valid and finding the extrema. To decrease numerical errors due to finite grid size at some δ₀ values, the optimized parameter values were used as initial values for the Nelder-Mead simplex algorithm (Gao & Han 2012) to obtain the final extrema. We repeated this process for a series of phase center declinations (δ₀), for the latitudes ϕ = −26.82°, 47.38°, and 52.91° (corresponding to SKA, NenuFAR, and LOFAR respectively). In the left panels of Fig. 3, we see that Eq. (7) (blue dashed line) agrees with the numerically maximized delays only for |δ₀ + ϕ|> π/2. For all other values of δ₀, the numerically maximized delay is 2|b|/c. Eq. (8) (red solid line) agrees with the numerically maximized delays in both regimes |δ₀ + ϕ|> π/2 and |δ₀ + ϕ|≤π/2. It should be noted that even if Eq. (6) is not imposed, we get the same result, thus confirming the claim that the maximization results should be unaffected by whether or not we impose the condition that the source and phase center lie above the physical horizon.

Fig. 3. Maximum delay in visibilities measured by a 100 m baseline, for a given phase center declination, assuming a full synthesis. Left: maximum delay as a function of declination of the phase center δ0, obtained by numerical maximization of Eq. (5) over H, δ, H0, θ (black dots), the maximum delays predicted by the partial equation (Eq. (7): blue dashed line) and those predicted by the full equation (Eq. (8): red solid line). Right: values of H, δ, and H0 at the maximum. The vertical dotted lines separate regimes where |δ0 + ϕ|≤π/2 and |δ0 + ϕ|> π/2.

In the right panels of Fig. 3, we see that in the regime |δ₀ + ϕ|≤π/2, H and H₀ are separated by π and the plot of δ against δ₀ has a slope of −1. This shows that the source vector and the phase center vectors are in diametrically opposite directions as explained previously. Also, in this regime, H and H₀ will give the same delay for two values separated by π. These π jumps between consecutive points have been removed while plotting for visual clarity, along with any 2π jumps. In the regime, |δ₀ + ϕ|> π/2, the source and phase center are no longer in opposite directions, and δ|_max agrees with the values predicted by Eq. (7). Also, H₀ assumes a value of π (except at δ₀ = π/2 or −π/2 where the value of H₀ is irrelevant). Again, this confirms that the assumption made in deriving Eq. (7) is correct in this regime.

3.3. Snapshot at a given LST

In the previous section, we investigated the situation of a full synthesis where the phase center moves with respect to the array and derived an expression for the maximum delay over this 24 h period. In this section, we instead consider a snapshot observation at a given LST, where the phase center is fixed with respect to the array. The previous analysis shows that in both regimes (|δ₀ + ϕ|≤π/2 and |δ₀ + ϕ|> π/2), the maximum delay occurs due to a source on the physical horizon, lying on the plane containing $\widehat{\mathit{z}}$ and $\widehat{\mathit{p}}$, and in the opposite direction of the phase center. Using this insight, we can derive a general equation for the horizon delay for an arbitrary H₀. At a given local sidereal time (LST), the maximum additional delay due to phasing is equal to the W coordinate of a baseline lying in this plane. This additional delay is $\mathit{b}·\widehat{\mathit{p}}=|\mathit{b}|cos{a}_0$, where a₀ is the altitude of the phase center. Thus, the equation for the maximum delay is given by:

$$\begin{array}{c}\hfill |{\eta }_0{|}^{\max }=\frac{|\mathit{b}|}{c}(1+cos{a}_0).\end{array}$$

For a given LST (denoted by t) and right ascension (RA) of the phase center (α₀), this reduces to:

$$\begin{array}{c}\hfill |{\eta }_0{|}^{\max }=\frac{|\mathit{b}|}{c}(1+\sqrt{1-(sin\varphi sin{\delta }_0+cos\varphi cos{\delta }_{0}cos(t-{\alpha }_0){)}^2}).\end{array}$$(9)

We verify this equation by maximizing Eq. (5) over H, δ, and θ at a particular H₀, δ₀, both with and without imposing the horizon condition of Eq. (6). We repeated this for multiple latitudes to verify Eq. (9) for a range of combinations of ϕ, H₀, and δ₀. Figure 4 shows the maximum delay as a function of the phase center direction and latitude. The left column shows the values obtained by maximizing Eq. (5) at each H₀, δ₀. The right column shows the expected maximum delays obtained using Eq. (9). We find that the equation can predict the maximum delays, which confirms the result of Sect. 3.2 and the premise in deriving Eq. (9) that the maximum occurs due to a source on the physical horizon, lying on the plane containing zenith and the phase center and in the direction opposite to the phase center. Also, in the regime of |δ₀ + ϕ|≤π/2, the maximum values of the delay occur along the edge of the striped regions, which indicate the physical horizon. Thus, we affirm the result from Sect. 3.2 that in the regime of |δ₀ + ϕ|≤π/2, the maximum delay occurs when H₀ and δ₀ are such that the phase center lies on the horizon. In the regime of |δ₀ + ϕ|> π/2, we see that the maximum value of the delay (2|b|/c = 6.66 × 10⁻⁷ s for |b| = 100 m) is not reached at any H₀, since the phase center can never be on the horizon for this declination range. Also, the maximum occurs at H₀ = π, which agrees with the result from Sect. 3.2 and confirms again that the assumption made in deriving Eq. (7) is valid in this regime.

Fig. 4. Maximum delay in visibilities measured by a 100 m baseline for a given phase center location in the sky (H0, δ0). Left: Maximum delays, obtained by numerical maximization of Eq. (5) over H, δ, and θ. Right: Maximum delays expected from Eq. (9). The vertical dotted lines separate regimes where |δ0 + ϕ|≤π/2 and |δ0 + ϕ|> π/2. The slanted stripes indicate the regions where the phase center is below the physical horizon.

We now have two equations that give us the maximum delays for a given phase center, one assuming full synthesis (Eq. (8)) and the other for a snapshot at a particular LST (Eq. (9)). Both equations have been derived assuming that the sky is filled with sources, the array has baselines in all directions, and the baselines lie on the plane perpendicular to the zenith.

4. Horizon wedge in the cylindrical power spectrum

In the reconstructed power spectrum approach, the visibilities are phased to the phase center direction and then gridded in the uv plane before Fourier transformation along frequency. The gridded data is then used to construct the cylindrical power spectrum. In this section, we derive the equations governing the horizon line in the cylindrical power spectrum.

We work in the same ABC coordinate system as Sect. 2 with the C axis pointing in the direction of the phase center. Rewriting Eq. (2) after dropping the pq suffix, we get:

$$\begin{array}{c}\hfill V(\nu )=e^{-2\pi i\mathit{b}·(\widehat{\mathit{s}}-\widehat{\mathit{p}})\nu /c}=e^{-2\pi i(Ul+Vm+W(n-1))\nu /c}.\end{array}$$(10)

Under the flat-sky approximation l, m ≪ 1. Hence, $n=\sqrt{1-l^2-m^2}\approx 1$, and Eq. (10) becomes:

$$\begin{array}{c}\hfill V(\nu )\approx e^{-2\pi i(Ul+Vm)\nu /c}=e^{-2\pi i{\mathit{b}}^{\prime }·{\mathit{s}}^{\prime }\nu /c},\end{array}$$(11)

where b′=(U, V) and s′=(l, m) are both two-dimensional vectors lying in the tangent plane of the sky perpendicular to $\widehat{\mathit{p}}$. Here s′ is not a unit vector. For a source located at an angle ψ from the phase center $\widehat{\mathit{p}}$, |s′| = sin ψ, and for large ψ, the flat-sky assumption is broken.

To construct the cylindrical power spectrum, the visibilities are placed on the uv plane, where u = Uν/c and v = Vν/c. In the case of an ideal interferometer with full uv coverage, for each point in the uv plane, there is always a baseline crossing that point for each frequency channel. In this case, there will be no fluctuation in the visibility amplitude along frequency at a given uv point, because the entire frequency dependence has been absorbed in the uv term. Thus the entire foreground power is confined to zero delay and there is no foreground wedge. However, no interferometer is ideal, since the uv coverage is sparse and the measurements are sampled at a discrete set of frequencies. Thus, a uv grid with a finite-sized uv cell is chosen and the observed visibilities are placed on this grid using a convolution kernel. Therefore, the gridded visibility at a certain uv cell has contribution from visibilities lying within the area of support of the convolution function (e.g., Offringa et al. 2019). So even in the sky-plane power spectrum, we end up with a horizon wedge.

Since we do not want to impose the flat-sky approximation, we must work directly with Eq. (10). Let us consider a baseline for which we reach the extreme case where it stays in the same uv cell for all frequencies¹². This extreme case corresponds to the situation of the maximum delay due to mode mixing, whereas if multiple baselines cross the uv cell at different frequencies, the mode-mixing effect would be lower. Suppose the baseline crosses the point u,v at the central frequency ν₀ (u = Uν₀/c, v = Vν₀/c). This visibility in such a uv cell is now Fourier transformed along frequency:

$$\begin{array}{cc}\hfill V(\eta )& ={\displaystyle {\int }_{\nu }}{e}^{-2\pi i\mathit{b}·(\widehat{\mathit{s}}-\widehat{\mathit{p}})\nu /c}{e}^{2\pi i\eta \nu }\mathrm{d}\nu \hfill \\ \hfill & =\delta (\eta -{\eta }_0),\mathrm{where}{\eta }_0=\mathit{b}·(\widehat{\mathit{s}}-\widehat{\mathit{p}})/c.\hfill \end{array}$$(12)

We note that η₀ here has the same expression as the signal delay after phasing from Eq. (3). In Sect. 3.1, we have derived the maximum values of η₀ for two cases: full synthesis in Sect. 3.2 and at a given LST in Sect. 3.3. We consider these two cases separately for the horizon line equation.

4.1. Snapshot at a given LST

We first consider the case of a particular LST. Inserting the expression for maximum delay¹³ from Eq. (9) in Eq. (12) gives ${\eta }_0^{\max }=|\mathit{b}|(1+cos{a}_0)/c$. Now, from Sect. 3.3, we know that the maximum delay occurs due to a source lying on the physical horizon, in the opposite direction of the phase center, on the plane containing $\widehat{\mathit{z}}$ and $\widehat{\mathit{p}}$. The baseline which will have the maximum delay due to this source will be oriented such that it lies on this same plane. The length of the projection of such a baseline on the plane perpendicular to the phase center is given by |b′| = |b|sin a₀ where a₀ is the altitude of the phase center at that LST. The maximum delay is then given by:

$$\begin{array}{c}\hfill {\eta }_0^{\max }=\frac{|{\mathit{b}}^{\prime }|}{csin(a_0)}(1+cos{a}_0)=\frac{\sqrt{u^2+v^2}}{{\nu }_0}\left(\frac{1+cos{a}_0}{sin{a}_0}\right).\end{array}$$

The parameters $\sqrt{u^2+v^2}$ and η are related to the cosmological coordinates as (Morales & Hewitt 2004; McQuinn et al. 2006):

$$\begin{array}{c}\hfill \sqrt{u^2+v^2}=k_{\perp }\frac{D_M(z)}{2\pi };\eta =k_{\Vert }\frac{c{(1+z)}^2}{2\pi {\mathcal{H}}_0{\nu }_{21}E(z)},\end{array}$$(13)

where D_M(z) is the conversion factor from angular units to comoving distance units, ℋ₀ is the Hubble constant, ν₂₁ = 1420 MHz and E(z) is the dimensionless Hubble parameter. Converting to cosmological coordinates, we get:

$$\begin{array}{cc}& k_{\Vert }^{\mathrm{LST}}=k_{\perp }\frac{D_M(z){\mathcal{H}}_{0}E(z)}{c(1+z)}\left(\frac{1+cos{a}_0}{sin{a}_0}\right),\hfill \\ \hfill & \mathrm{where}{a}_0={sin}^{-1}(sin\varphi sin{\delta }_0+cos\varphi cos{\delta }_{0}cos(t-{\alpha }_0)).\hfill \end{array}$$(14)

4.2. Full synthesis

For full synthesis, we do a similar treatment on Eq. (8). The horizon delay is then:

$$\begin{array}{c}\hfill {\eta }_0^{\max }=\frac{|{\mathit{b}}^{\prime }|}{csin{a}_0}(1+sin(max\{|\varphi +{\delta }_0|,\pi /2\})).\end{array}$$

From Sect. 3.2, we know that for |δ₀ + ϕ|≤π/2, the maximum delay occurs when the phase center is at the horizon. Thus, a₀ = 0 and the horizon line has a slope of infinity. For |δ₀ + ϕ|> π/2, the maximum delay occurs when the phase center crosses the meridian plane closest to the horizon (H₀ = π). Thus, sin a₀ = sin ϕ sin δ₀ + cos ϕ cos δ₀ cos H₀ = −cos(δ₀ + ϕ). Inserting these into the expression for maximum delay and converting to cosmological coordinates using Eq. (13), the horizon line equation takes the form:

$$\begin{array}{c}\hfill k_{\Vert }^{\mathrm{syn}}=\{\begin{array}{cc}{k}_{\perp }\frac{D_M(z){\mathcal{H}}_{0}E(z)}{c(1+z)}\left(\frac{1+sin(|{\delta }_0+\varphi |)}{-cos({\delta }_0+\varphi )}\right)\hfill & \mathrm{if}|{\delta }_0+\varphi |>\pi /2\hfill \\ \mathrm{Undefined}\hfill & \mathrm{if}|{\delta }_0+\varphi |\le \pi /2.\hfill \end{array}\end{array}$$(15)

It should be noted that a full synthesis observation is seldom performed over a single day. However, for target fields close to the celestial poles, and for deep integrations spanning several hundred hours over multiple months, it is likely that the full 24 h of LST will be covered. The power spectrum estimated from the combined data will have foreground power reaching up to the full synthesis line. Equation (15) thus indicates the maximum possible foreground wedge extent for a chosen phase center, without considering additional details of observation parameters.

We now have two equations: Eq. (14) and Eq. (15), that describe the horizon line in the cylindrical power spectrum and have been derived without imposing a flat-sky approximation. It should be noted that the cylindrical power spectrum constructed in Cartesian coordinates is not the correct basis for describing curved sky effects, and this can be achieved by using a spherical harmonic spatial basis. However, if the cylindrical power spectrum is still used, these derived equations correctly describe the extent of full-sky foregrounds in it. This is the case for most tracking interferometers where the main lobe of the primary beam is relatively narrow and the 21 cm signal of interest is dominated by its footprint, while the bright foregrounds picked up by the primary beam sidelobes can still contaminate the power spectrum up to the horizon lines given by these equations. In Sect. 6.1, we use full-sky simulations to validate both equations.

4.3. Special case of zenith phasing

Starting with Eq. (11), which assumes a flat-sky approximation and considering an extreme case of a baseline that remains in the same uv cell for all frequency channels as in the previous section, we get ${\eta }_0=\mathit{b}\prime ·\mathit{s}\prime /c$. Given a certain uv length, the delay η₀ due to that source s′ has the maximum value for a uv cell such that b′ points in the direction of s′. Thus, the maximum delay due to this source is given by η₀ = |b′||s′|/c and the wedge equation under a flat-sky approximation is given by:

$$\begin{array}{c}\hfill k_{\Vert }=k_{\perp }\frac{D_M(z){\mathcal{H}}_{0}E(z)}{c(1+z)}|{\mathit{s}}^{\prime }|=k_{\perp }\frac{D_M(z){\mathcal{H}}_{0}E(z)}{c(1+z)}sin\psi .\end{array}$$(16)

Setting ψ = π/2 for the horizon (maximum k_∥), it follows that:

$$\begin{array}{c}\hfill k_{\Vert }^{\mathrm{flat}}=k_{\perp }\frac{D_M(z){\mathcal{H}}_{0}E(z)}{c(1+z)}.\end{array}$$(17)

This equation has conventionally been used to define the horizon lines in cylindrical power spectra of 21 cm cosmology. However, setting ψ = π/2 in Eq. (16) violates the flat-sky approximation that was used to derive the equation in the first place. As a result, the equation incorrectly suggests that the maximum k_∥ value occurs at ψ = π/2, which does not correspond to the physical horizon, except when the phase center points exactly toward the zenith. For zenith phasing (a snapshot drift scan observation) with an array lying on the Earth’s tangent plane, all W coordinates are zero. So even without imposing the flat-sky approximation, Eq. (10) reduces to Eq. (11). We note that for the special case of a zenith phasing, a₀ = π/2, and the LST-dependent horizon line derived in this analysis (Eq. (14)) reduces to the conventional horizon line equation (Eq. (17)).

5. Source wedge in the cylindrical power spectrum

The horizon line is useful to understand the maximum extent of foreground contamination in the data. However, in many situations, the sky might be dominated by a few very strong sources far from the phase center. It is thus very useful to be able to predict the impact of a specific source on the cylindrical power spectrum. In this section, we derive an expression for the Fourier modes most strongly affected by a given source in the cylindrical power spectrum.

5.1. Effect of sampling and gridding

Under the flat-sky approximation, the delay due to a source lying in the plane perpendicular to the LoS (s′) is given by η₀ = b′⋅s′/c (Sect. 4.3). For the horizon delay, we used the maximum value of η₀, which occurs when b′ and s′ point in the same direction. However, in this section, the aim is not to determine the maximum delay caused by a source, but rather to identify the delay at which the source has the maximum power. This is important for analyzing the cause of artifacts in the cylindrical power spectrum. If β is the angle between b′ and s′, all delays between zero (for β = π/2) and |b′||s′|/c (for β = 0) have the same power. Thus, the signature of a single source at an angle ψ from the phase center is a wedge with equal power along lines of k_∥ ∝ k_⊥, with the maximum slope given by Eq. (16).

However, this simple situation is only applicable to an interferometer with perfect uv coverage. An imperfect uv coverage gives rise to a point spread function (PSF) that is not a delta function. If the uv tracks are circular, the PSF sidelobes also have a circular symmetry. Given a certain FoV within which we make the image, the PSF sidelobes of sources far from the phase center run through the field. For sources sufficiently far away compared to the FoV, the signature of the sidelobes in the uv plane approaches a straight line of higher source power along the direction of the source, that is, perpendicular to the sidelobe ripples. Thus, instead of uniform power in the uv plane, we end up with higher power in the uv plane in the direction of the source (as illustrated in the middle and right panels of Fig. 5).

Fig. 5. Signature of a source on the uv plane. The left panel is a schematic diagram illustrating the quasi-sinusoidal behavior of the phase of visibilities due to a source along a circular uv track. The middle and right panels show the visibility amplitude in the uv plane for simulations of Cyg A with the phase center at the NCP, with FoVs of 8° and 32°, respectively.

In the uv domain, we can explain this as follows. The visibility due to a source is a complex exponential in the uv plane and an interferometer samples this visibility along uv tracks. These sampled visibilities are gridded on the uv plane with a convolution kernel whose width is decided by the target FoV. Now the visibility (both real and imaginary parts) along a given uv track has slower fluctuations in the direction of the source and faster fluctuations in the perpendicular direction. The left panel of Fig. 5 is a schematic diagram of the phase of the visibilities due to a source at a certain angular distance from the phase center located at the bottom right corner of the image, after sampling at a circular uv track. We highlight this track and the uv cells not sampled are shown with faded colors. Here we can see a quasi-sinusoidal behavior of the sampled visibility as a function of the angle along the uv track, with the slowest variations in the direction of the source. Thus, a convolution kernel of a given width suppresses the amplitude of the sampled visibility more in the direction perpendicular to the source direction (Offringa et al. 2012).

The middle and right panels of Fig. 5 show the visibility amplitude in the uv plane for a 12 h simulated observation of Cygnus A (Cyg A) with the phase center at the NCP. We see that the uv plane has higher power in the direction of Cyg A (lying toward the bottom right of the image). The gridded uv plane is constructed using only the data up to a maximum FoV from the phase center. The visibility amplitudes in the uv plane for FoVs of 8° (middle panel) and 32° (right panel) show that the signature of the source in the uv plane is wider for a larger FoV. This is because a smaller FoV leads to a wider convolution kernel, which suppresses fluctuations as we move slightly away from the source direction. For a large number of uv tracks, this results in a sharper feature in the uv plane compared to that of a wider FoV. This effect is also discussed in Appendix B, where we illustrate the effect of sampling and gridding analytically, for a circular uv track. It should be noted that the line of highest power is oriented exactly along the source direction only for circular uv tracks because the real and imaginary parts of the visibility along a uv track have the slowest fluctuations exactly in the source direction only if the uv track is circular. For elliptical uv tracks, for example, the slowest fluctuations will not occur exactly along the source direction. This results in a line of highest power that is not aligned along the source direction, with the extent of the deviation being higher for ellipses with higher eccentricity. Additionally, the functional form of the gridding kernel affects the direction of maximum power in the uv plane. Both these effects are discussed in Appendix B.

5.2. Delay due to a source

In the previous section, we discussed how the sampling and gridding performed in the process of measurement by an interferometer and construction of gridded visibility cubes results in a higher amplitude of the gridded visibilities along the source direction in the uv plane. We now derive an expression for the delay at which this higher amplitude occurs. Working in the same XYZ coordinate system as Sect. 3.1, the delay due to a source located at H, δ is given by Eq. (5). In the absence of the exact uv coverage, throughout this section and in the entire derivation of source line equations, we make the simplifying assumption that the uv tracks are circular. Thus, the highest power is due to a baseline pointing in the direction of the source from the phase center (i.e., along $\widehat{\mathit{s}}-\widehat{\mathit{p}}$)¹⁴. Since we assume that the baselines lie on the tangent plane of the Earth, the expression for $\mathit{b}·(\widehat{\mathit{s}}-\widehat{\mathit{p}})$ for such a baseline is given by the baseline length multiplied by the length of the projection of $(\widehat{\mathit{s}}-\widehat{\mathit{p}})$ on the plane perpendicular to $\widehat{\mathit{z}}$. The corresponding delay (${\eta }_0=\mathit{b}·(\widehat{\mathit{s}}-\widehat{\mathit{p}})/c$) is given by:

$$\begin{array}{cc}\hfill {\eta }_0=& \frac{|\mathit{b}|}{c}|(\widehat{\mathit{s}}-\widehat{\mathit{p}})-((\widehat{\mathit{s}}-\widehat{\mathit{p}})·\widehat{\mathit{z}})\widehat{\mathit{z}}|\hfill \\ \hfill =& \frac{|\mathit{b}|}{c}({(cos{\delta }_{0}sin{H}_0-cos\delta sinH)}^2+(sin\varphi (cos{\delta }_{0}cos{H}_0\hfill \\ \hfill & -cos\delta cosH)+cos\varphi (sin\delta -sin{\delta }_0){)}^2{)}^{1/2}.\hfill \end{array}$$(18)

5.3. Source line in the cylindrical power spectrum

We next derive the equation describing the line in the cylindrical power spectrum due to a given source. From Eq. (12), the delay in a uv cell without imposing the flat-sky approximation is $\mathit{b}·(\widehat{\mathit{s}}-\widehat{\mathit{p}})/c$. Inserting the expression for the delay at which a specific source $\widehat{\mathit{s}}$ contributes the maximum power (from Eq. (18)), we get:

$$\begin{array}{c}\hfill {\eta }_0^{\mathrm{source}}=\frac{|\mathit{b}|}{c}|(\widehat{\mathit{s}}-\widehat{\mathit{p}})-((\widehat{\mathit{s}}-\widehat{\mathit{p}})·\widehat{\mathit{z}})\widehat{\mathit{z}}|\equiv \frac{|\mathit{b}||\mathit{a}|}{c},\end{array}$$(19)

where $\mathit{a}=(\widehat{\mathit{s}}-\widehat{\mathit{p}})-((\widehat{\mathit{s}}-\widehat{\mathit{p}})·\widehat{\mathit{z}})\widehat{\mathit{z}}$. It should be noted that a points in the direction of the baseline responsible for the maximum power. Now, the projected length of the baseline b on the plane perpendicular to the phase center is given by |b′| = |b|sin γ, where γ is the angle between b and $\widehat{\mathit{p}}$. Since b points along $\widehat{\mathit{a}}=\mathit{a}/|\mathit{a}|$, the relation between |b′| and |b| becomes:

$$\begin{array}{c}\hfill |\mathit{b}\prime |=|\mathit{b}|\sqrt{1-(\frac{\mathit{a}}{|\mathit{a}|}·\widehat{\mathit{p}}{)}^2}.\end{array}$$(20)

Combining Eqs. (19) and (20), we get:

$$\begin{array}{c}\hfill {\eta }_0^{\mathrm{source}}=\frac{\sqrt{u^2+v^2}}{{\nu }_0}\left(\frac{{|\mathit{a}|}^2}{\sqrt{{|\mathit{a}|}^2-{(\mathit{a}·\widehat{\mathit{p}})}^2}}\right).\end{array}$$(21)

Finally, converting to cosmological coordinates using Eq. (13) and expanding the expressions for the vectors in terms of H, δ, H₀, δ₀, and ϕ, the source line equation in the cylindrical power spectrum takes the form:

$$\begin{array}{cc}\hfill k_{\Vert }^{\mathrm{source}}& =k_{\perp }\frac{D_M(z){\mathcal{H}}_{0}E(z)}{c(1+z)}\hfill \\ \hfill \times & \frac{\begin{array}{cc}& {(sinHcos\delta -sin{H}_{0}cos{\delta }_0)}^2+(sin\varphi cos\delta cosH\hfill \\ \hfill & -sin\varphi cos{\delta }_{0}cos{H}_0-sin\delta cos\varphi +sin{\delta }_0{cos\varphi )}^2\hfill \end{array}}{\sqrt{\begin{array}{cc}& {(sinHcos\delta -sin{H}_{0}cos{\delta }_0)}^2-((sinHcos\delta -sin{H}_{0}cos{\delta }_0)\hfill \\ \hfill & sin{H}_{0}cos{\delta }_0+(sin\varphi cos{\delta }_{0}cos{H}_0-sin{\delta }_{0}cos\varphi )(sin\varphi cos\delta cosH\hfill \\ \hfill & -sin\varphi cos{\delta }_{0}cos{H}_0-sin\delta cos\varphi +sin{\delta }_0{cos\varphi ))}^2+(sin\varphi cos\delta \hfill \\ \hfill & cosH-sin\varphi cos{\delta }_{0}cos{H}_0-sin\delta cos\varphi +sin{\delta }_0{cos\varphi )}^2\hfill \end{array}}}.\hfill \end{array}$$(22)

By substituting H = t − α and H₀ = t − α₀, we get the expression of the line due to a given source in the cylindrical power spectrum as a function of LST, for a particular phase center. In Sect. 6.2, we verify the validity of this equation against representative simulations for a range of observation parameters. We note that in Eq. (22), the slope of the line in the (k_⊥, k_∥) space is a function of LST. Hence, given an observation duration, the line due to a source will assume different slopes at different times. This will lead to a wedge-like region in the cylindrical power spectrum most affected by a source during the observation. This will be illustrated in Sect. 6.2.

5.4. Special case of zenith phasing

Under a flat-sky approximation, the expression of the delay in a uv cell due to a source s′ is b′⋅s′/c (Sect. 4.3). We showed in Sect. 5.1 that the maximum power due to a source occurs on uv cells that lie along the direction of s′. Thus η₀^source = |b′||s′|/c and the signature of the highest power due to a source under the flat-sky approximation is given by:

$$\begin{array}{c}\hfill k_{\Vert }^{\mathrm{source}}=k_{\perp }\frac{D_M(z){\mathcal{H}}_{0}E(z)}{c(1+z)}sin\psi .\end{array}$$(23)

However, as in the case of the horizon line equation with flat-sky approximation (Eq. (17)), this equation breaks down for sources far away from the phase center where the condition ψ → 0 does not hold.

As discussed in Sect. 4.3, the special case of zenith phasing is equivalent to imposing a flat-sky approximation. When the phase center is at the zenith, $\widehat{\mathit{z}}=\widehat{\mathit{p}}$. If ψ is the angle between $\widehat{\mathit{s}}$ and $\widehat{\mathit{p}}$, $\mathit{a}=\widehat{\mathit{s}}-\widehat{\mathit{p}}cos\psi$. It follows that |a|=sin ψ and $\mathit{a}·\widehat{\mathit{p}}=0$. Inserting these expressions in the source line equation derived in this analysis (Eq. (21)) and converting to cosmological coordinates, we get back the conventional source line equation (Eq. (23)).

6. Simulations

In Sect. 4, we derived two equations (Eqs. (14) and (15)) which describe the horizon line in the cylindrical power spectrum at a given LST and for full synthesis. In Sect. 5, we derived an equation that describes the source line in the cylindrical power spectrum at a given LST (Eq. (22)). These equations are dependent on several parameters, such as the latitude of the telescope, the coordinates of the phase center, and the time of observation. We now verify the validity of these equations over these parameters through a set of representative simulations for a range of observation parameters and a set of latitudes and longitudes of current and upcoming radio telescopes.

6.1. Horizon line

In this section, we verify both Eqs. (14) and (15) against full-sky simulations. The sky model contains 1000 sources uniformly distributed in the sky with the flux (S) of the sources described by a flux density distribution dN/dS ∝ S^−1.54 (Franzen et al. 2016) between 1 and 10 Jy, with each component having a ν^−0.8 frequency dependence. This model is sufficient to assess the effect of sources from all over the sky down to the horizon for any location on the Earth, at any LST. The sky model used in the simulations is shown in the top panel of Fig. 6. We constructed a fiducial array that (locally) has a random distribution of 79 stations in a 2 km diameter region, with the vertical positions having a jitter of ±2 m. This array was then placed at the geographical locations of NenuFAR, LOFAR, and SKA. The array configuration and the baseline coverage (for NCP phasing with the array placed at the location of NenuFAR) are shown in the bottom panels of Fig. 6.

Fig. 6. Inputs to the simulations used for validating the derived equations. Top: sky model used in the simulations performed for validating the horizon line equations. Bottom left: Local array configuration. Bottom right: example baseline coverage for a 1 h simulation with NCP phasing with the array placed at the location of NenuFAR.

The simulations were performed using OSKAR (Dulwich 2020). Since the aim is to test the derived equations that describe the maximum extent of foregrounds without considering primary beam effects, we disabled station beam effects and used “station type = Isotropic beam” for all simulations. In these simulations, the visibility of a source is not predicted when it goes below the physical horizon, as would be expected in an actual observation. The frequency range was chosen to be 61–72 MHz with 177 channels and the time resolution was chosen to be 4 s¹⁵. We performed 1 h simulations centered at six LSTs each separated by 4 h (and hence spanning a day). We repeated the simulations for three phase centers and chose the phase center RA to be 1 h such that the phase center crosses the local meridian at t = 1 h and thus at the maximum altitude. It should be noted that since a visibility and its complex conjugate are sampled at u,v and −u,−v respectively, 12 hours is sufficient to yield a complete uv coverage. However, the situations separated by 12 h are not identical since the phase center can be at different elevations, thus requiring different geometric phasing. The settings used for the simulations for validating the horizon lines are summarized in the top section of Table 2.

The power spectra were constructed using the power spectrum pipeline pspipe¹⁶, which is used for estimating the power spectra for the NenuFAR and LOFAR 21 cm cosmology projects. pspipe uses WSClean (Offringa et al. 2014) to construct gridded dirty image cubes. A Hann filter with a FoV of 8 deg was applied in the spatial direction and a Blackman-Harris filter along the frequency direction of the dirty image cube before Fourier transformation into the uvη space. Finally, the modulus squared of the gridded data cube in the Fourier space was used to construct the cylindrical power spectrum.

Figure 7 shows the cylindrical power spectra for the simulations corresponding to NenuFAR, LOFAR, and SKA for a phase center declination of 70° (−70° for SKA). The results of the remaining simulations corresponding to the celestial pole and zenith transit phase centers are shown in Fig. E.1. For each situation, both the LST-dependent line (Eq. (14): solid black line) and the full synthesis line (Eq. (15): red dashed line) are shown. The LST-dependent horizon line (Eq. (14)) will in general move through the (k_⊥, k_∥) space slightly during the 1 h observation. So the horizon k_∥ is evaluated at a series of 100 equally spaced time points within the 1 h observation duration and the line corresponding to the maximum k_∥ is shown. We see that both the LST-dependent line and the full synthesis line perfectly predict the full sky horizon limit in all situations in both Figs. 7 and E.1. For each row, the full synthesis line indicates the maximum possible slope of the LST-dependent horizon line, that is, the horizon line if we were to perform a full synthesis. The altitude of the phase center (a₀) at the LST at which the simulation is centered is mentioned on each plot. We see that the horizon line is steeper when the phase center goes closer to the horizon, as predicted by Eq. (14). The horizon extent is lower for the LOFAR simulations compared to the NenuFAR simulations. This is because LOFAR is at a higher latitude, and the same phase centers are closer to the zenith compared to NenuFAR. The conventional horizon line (Eq. (17): black dotted line) fails to correctly describe the full sky foreground wedge horizon except in situations where the phase center is very close to the zenith (bottom left panel for each location in Fig. E.1). In these situations, the LST-dependent horizon line equation (Eq. (14)) reduces to the conventional horizon line equation (Eq. (17)), as described in Sect. 4.3. Also, the conventional equations predict that the horizon line is independent of time, the coordinates of the phase center, and the latitude of the telescope. However, in Figs. 7 and E.1, we can see that this is not the case for full-sky simulations with an arbitrary phase center, where the extent of the foreground wedge is highly dependent on these parameters.

Fig. 7. Cylindrical power spectra for full-sky simulations with the array placed at the locations of NenuFAR, LOFAR, and SKA corresponding to a phase center declination of 70° (−70° for SKA). The different columns correspond to different LSTs at which the 1 h simulations are centered. The dotted line is the conventional horizon line. The black solid line is the LST-dependent horizon line (Eq. (14)). The red dashed line is the horizon line for full synthesis (Eq. (15)). a0 is the altitude of the phase center at the LST at which the simulation is centered. The results of the remaining simulations for the celestial pole and zenith transit phase centers are shown in Fig. E.1.

6.2. Source line

In this section, we verify Eq. (22) against simulations of a single source for a range of observation parameters. In the simulations, we used the same random array at the locations of NenuFAR, LOFAR, and SKA. The simulation and power spectrum generation were performed in the same manner as described in Sect. 6.1, except here the sky model consists of a single delta function at the position of a chosen source. Also, here we performed simulations for six LST ranges: three of them being 1 h observations at different LSTs separated by 4 h, and the other three being observations of different durations with the same start LST. This allows us to study how a longer observation duration results in a range of modes in the power spectrum getting affected by a strong source. The simulations were repeated for two sources, two phase centers, and three array locations corresponding to NenuFAR, LOFAR, and SKA. We chose Cyg A and Cassiopeia A (Cas A) for NenuFAR and LOFAR because they have the strongest impact on the data of these telescopes. For SKA, we chose Virgo A (Vir A) and Centaurus A (Cen A) and a phase center RA of 13 h to ensure that they remain above the horizon for all the observation durations considered. The settings for which the simulations for validating the source line were repeated are summarized in the bottom section of Table 2.

Figure 8 shows the cylindrical power spectra for the simulations used for validating the source line equations corresponding to Cyg A (for NenuFAR and LOFAR) and Vir A (for SKA). The remaining simulations of Cas A (for NenuFAR and LOFAR) and Cen A (for SKA) are shown in Fig. E.2. In each situation, the source line is evaluated from Eq. (22) at 100 equally spaced time samples and the maximum and minimum source lines for the observation duration are shown (in solid black lines). The dotted black lines indicate the conventional source line derived with a flat-sky approximation (Eq. (23)). We see that in all situations, the new source lines provide a better depiction of the range of modes most strongly affected by a source in the cylindrical power spectrum. The conventional source lines only depend on the angular separation between the source and the phase center, and hence do not move with LST (i.e., along a row in Figs. 8 and E.2). On the other hand, the new equations capture the trend of the movement of the source lines in the (k_⊥, k_∥) space. For longer observations of 4, 8, and 12 h (the three right-most columns in Figs. 8 and E.2) the new equations can describe the range of modes affected by a source during a synthesis observation.

Fig. 8. Cylindrical power spectra for Cyg A simulations at the locations of NenuFAR and LOFAR, and Vir A simulations at the location of SKA. For each location, the two rows correspond to two phase center declinations (δ0) of the celestial pole (upper row) and zenith transit (lower row). The different columns correspond to the different LST ranges for which the simulations were repeated. The dotted line is the conventional source line. The black solid lines correspond to the maximum and minimum source lines for the observation duration, obtained from the source line equation derived in this analysis (Eq. (22)). The remaining Cas A and Cen A simulation results are shown in Fig. E.2.

An important distinction between the horizon line and the source line equations is that the horizon line equations (both Eqs. (14) and (15)) have the maximum delay in the numerator (from Eqs. (9) and (8)) divided by the projection of the baseline responsible for the maximum delay on the uv plane, in the denominator. Since the baseline responsible for the maximum delay lies on the plane containing the zenith and the phase center, this is also the baseline which has the smallest projection on the uv plane. This means that the numerator is the maximum possible delay and the denominator is the minimum possible projection, thus resulting in the maximum possible extent of power in the cylindrical power spectrum. However, for the source line equation (Eq. (22)), though the numerator corresponds to the maximum delay due to the source (see Appendix C), the denominator does not necessarily correspond to the minimum possible projection of a baseline on the uv plane. Thus, it is not surprising to have modes above the maximum source line that are affected by the power due to the source in some of the cases. Also, note that the power does not go to zero below the minimum source lines, since baselines oriented in other directions also contribute power, even though it is attenuated due to sampling and gridding (as described in Sect. 5.1).

In some situations, we see that the peak power extends beyond the source lines. This is possible because, in the derivation of the source line equation, we take the assumption that the maximum power in the uv plane is in the direction of the source. However, the actual signature of a source in the uv plane does not always have the highest power exactly in the source direction, and the exact functional form of the amplitude of the sampled gridded visibilities will depend on the shape of the uv tracks and also on the nature of the convolution function (see Sect. 5.1 and Appendix B). Therefore the power due to the source slightly away from the source direction can end up slightly away from the predicted source line direction in the cylindrical power spectrum. This effect will be more pronounced for phase centers very close to the physical horizon. This is when the denominator of Eq. (22) can have values very close to zero, and slight changes in the numerator will result in large variations of the k_∥ predicted by the source line equation. In Figs. 8 and E.2, we indeed see that the peak of the source line deviates from the predicted lines only for SKA simulations with the phase center at the SCP, which is very close to the physical horizon. It should be noted that phase centers that reach very close to the horizon are rarely chosen for 21 cm cosmology analyses since the sensitivity to the 21 cm signal decreases significantly due to attenuation by the dipole beam of receiving antennas. In Appendix D, we describe simulations that were performed for specific situations where the uv tracks run along the source direction and illustrate its impact on the uv plane and the power spectrum.

7. Discussion

The revised and new horizon and source line equations that we have developed accurately describe the full sky signature of foregrounds in the cylindrical power spectrum. In this section, we discuss the main limitations of the equations, their key characteristics, and their potential practical applications in 21 cm cosmology analyses.

7.1. Limitations

The following are the main limitations of the horizon and source line equations derived in this analysis.

Impact of the primary beam on horizon lines. The horizon line equations (Eqs. (15) and (14)) describe the maximum possible extent of foregrounds in the cylindrical power spectrum. However, we have not considered the impact of the primary beam of the interferometric elements in the derivation, which has a damping effect on the power far from the phase center and will decrease the extent of contamination at high k_∥, particularly so for instruments with a narrow primary beam extent. So while these lines mark the maximum boundary of the foreground wedge, the foreground power that reaches this threshold could be much lower than the thermal noise level.

Impact of uv coverage on source lines. The source delay line equations have been derived assuming circular uv tracks. However, this assumption breaks down for phase centers far away from the NCP or the SCP, particularly for latitudes close to the equator where the uv tracks can deviate considerably from being circular. We have seen in the simulations that the source line equations still describe these situations in general. However, in specific situations where the uv tracks deviate significantly from being circular, the signature of the source might not be well described by the derived equations. This is investigated in Appendix D for a specific situation where the uv coverage has radial structures in the direction of the source. The maximum and minimum of the source line equations shown in Figs. 8 and E.2 do not necessarily indicate the region in which the power will be confined in the parameter space. This is because, in the derivation of the equation, we assume that the maximum power in the uv plane occurs exactly along a line in the direction of the source. In reality, this is not always the case even for circular uv tracks. This can result in power spilling slightly beyond the minimum and maximum source lines for the observation duration, particularly for phase centers close to the physical horizon. This is discussed in more detail at the end of Sect. 6.2.

Non-ideal basis. The horizon and source line equations derived in this analysis describe the impact of foregrounds on the cylindrical power spectrum constructed from gridded data in the uvη domain. However, this Cartesian Fourier domain is not the ideal basis for wide-field instruments. Instead, a spherical harmonic power spectrum, with a spherical spatial Fourier basis, is more appropriate here (e.g., Liu et al. 2016; Ghosh et al. 2018). A cylindrical power spectrum is usually employed by 21 cm cosmology experiments because the uvν domain is a natural basis for interferometers. It should be noted that even for instruments with a narrow primary beam FoV, where the cylindrical power spectrum is appropriate for describing the 21 cm signal, foreground sources far from the phase center might still be visible in the primary beam sidelobes. As a result, curved sky effects still need to be taken into account, either by using the equations derived in this paper or by switching to a spherical spatial basis.

7.2. Features and practical applications

The equations derived in this analysis have several important features that give us new insights into the nature of foreground contamination in the power spectrum. Additionally, there are several situations where the equations derived in this analysis will significantly aid our objective of detecting the 21 cm signal.

Dependence of the foreground wedge on observation parameters. The conventional horizon and source line equations suggest that the foreground wedge and the signature of strong sources in the cylindrical power spectrum are independent of time, latitude, and the location of the phase center. We show that this is not the case if we consider foregrounds across the full sky, and the signature of the foregrounds depends on several parameters. The horizon line can even be vertical in the extreme example of the phase center crossing the physical horizon (Sect. 4.2), and in such cases, the foreground wedge covers the entire (k_⊥, k_∥) space. For the source line, the dependence on time results in a wedge-like signature in the cylindrical power spectrum for a synthesis observation. The derived equations (Eqs. (14), (15) and (22)) are presented as analytical functions of the observation parameters like the RA, Dec of the phase center, LST, and the latitude of the interferometer, making them directly applicable to actual observations.

Generalized version of conventional equations. The conventional equations for both the horizon line (Eq. (17)) and the source line (Eq. (23)) are special cases of the equations derived in this paper. Both the LST-dependent horizon line equation and the source line equation reduce to these conventional equations for zenith phasing (snapshot drift scan observations) where it does not matter whether we impose a flat-sky approximation. Therefore, the equations derived in this analysis serve as a more general version of the conventional equations.

Impact on foreground avoidance approach. Most radio interferometers aiming for a statistical detection of the 21 cm power spectrum are phased arrays. These are often also wide-field instruments that are sensitive to sources close to the physical horizon. Thus, the cylindrical power spectra estimated from their data will have foreground power well beyond the conventional horizon line extent. Foreground avoidance analyses involve the selection of an “EoR window” that is assumed to be free of foregrounds. The updated horizon line equations derived in this paper describe the actual maximum extent of foregrounds and will allow a more prudent choice of the “EoR window”.

Impact on foreground subtraction approach. In foreground subtraction approaches with Gaussian process regression (Mertens et al. 2018, 2024; Acharya et al. 2024), it is necessary to provide priors on the expected frequency coherence scale of foregrounds, to separate them from the 21 cm signal. These new horizon line equations will allow us to provide a more accurate prior on the frequency coherence scale as a function of k_⊥.

Planning observations to minimize the wedge. The horizon line equations suggest that the extent of the foreground wedge is closely related to the location of the phase center. This understanding will help in planning observations so that the extent of the horizon wedge is minimized. Planning observations such that the phase center is always above a certain altitude will decrease the extent of the foreground wedge at the cost of lower integration speeds due to limited time spent observing the field. The ideal trade-off will depend on the foreground mitigation strategy used, the modes of interest, and the thermal noise sensitivity of the instrument compared to expected 21 cm signals.

Distinguish between foregrounds and systematics. The source line equations will allow us to more accurately diagnose the impact of bright sources on the cylindrical power spectrum and correctly distinguish between features in the power spectrum caused by instrumental systematics and those caused by bright sources. Features that lie beyond the conventional horizon lines and were previously, maybe erroneously, assumed to be caused by instrumental or data processing effects could instead be caused by strong foreground sources far from the phase center.

Impact of bright sources on the uv plane. The understanding of the origin of the signature of sources in the gridded uv plane can help pinpoint the range of modes most strongly affected by a bright source. This can help in systematically flagging these sources in gridded data accurately. It could also be utilized to subtract the power due to foregrounds utilizing the intrinsic frequency smoothness of foregrounds. Additionally, it can serve as a diagnostic to distinguish between power caused by foregrounds and that caused by instrumental effects or radio frequency interference.

Planning observations to minimize the impact of strong sources. We have seen that in certain situations the source lines in the cylindrical power spectrum can be extremely steep. The maximum steepness is caused by sources lying on the physical horizon, in the direction opposite to the phase center. This steepness poses considerable challenges since it restricts access to the smallest k modes or largest scales where we have the highest sensitivity for foreground avoidance approaches. For foreground subtraction approaches, steep source lines make blind separation of foregrounds from the 21 cm signal difficult since the foregrounds now exhibit strong spectral fluctuations and are more difficult to separate from the 21 cm signal. The source line equations derived in this analysis can help plan observation phase centers or LST ranges in a way that the bright sources do not result in steep source lines which can extend to high k_∥ values.

8. Summary

In this paper, we derive and generalize the equations that describe the full-sky foreground wedge in the cylindrical power spectrum for 21 cm cosmology analyses. These equations are derived without taking a flat-sky approximation and accurately describe the signature of foregrounds in the cylindrical power spectrum for an arbitrary phase center. We derive two kinds of equations, corresponding to the horizon and source lines. Horizon lines define the maximum extent of foregrounds in the cylindrical power spectrum, whereas source lines describe the signature of a specific source. Both sets of equations are tested against full simulations for a fiducial interferometric array placed at the locations of NenuFAR, LOFAR, and SKA, across a range of observational parameters. We find that the foreground wedge can extend well beyond what the conventional horizon line equations suggest, and the equations derived in this paper perfectly describe this extended horizon limit. The updated source line equations significantly improve the characterization of the impact of a specific strong source on the cylindrical power spectrum. These new equations are crucial for correctly describing the signature of foregrounds in 21 cm cosmology analyses employing either foreground avoidance or foreground subtraction strategies. Additionally, the equations will allow a more informed choice of the phase center and LST ranges of observations, particularly in the context of wide-field instruments such as the NenuFAR, LOFAR, MWA, and the future SKA, that employ the reconstructed power spectrum approach.

We have made a Python library pslines¹⁷ openly available. This generic library can generate and plot the horizon and source lines using the new equations presented in this paper, for any tracking interferometer.

Acknowledgments

SM, LVEK, SAB, JKC, BKG, SG, and CH acknowledge the financial support from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (Grant agreement No. 884760, “CoDEX”). FGM acknowledges the support of the PSL Fellowship. EC would like to acknowledge support from the Centre for Data Science and Systems Complexity (DSSC), Faculty of Science and Engineering at the University of Groningen.

References

Appendix A: Maximum delays for NCP phasing

The NCP is the main target field for the LOFAR and NenuFAR 21 cm projects. Here we analytically derive the expression for the maximum delay for NCP phasing, by maximizing the expression for delay (Eq. (5)) subject to the horizon condition (Eq. (6)). We work in the XYZ coordinate system described in Sect. 3.1 and set δ₀ = 90° and ϕ ≥ 0 for NCP phasing. To make the problem analytically tractable, we reduce it to the meridian plane. It should be noted that this can only be done without loss of generality when the phase center is at the NCP or SCP since these are the only two points in the sky that are fixed relative to the Earth. The baseline vector and source responsible for the maximum delay lie on the meridian plane. We choose a baseline pointing toward the north and a source toward the south to get the extreme case of maximum absolute delay. This amounts to setting H = 0° and θ = 90° in Eqs. (5) and (6). Thus, we need to maximize the equation:

$$\begin{array}{c}\hfill {\eta }_0(\varphi ,\delta )=\frac{|\mathit{b}|}{c}(-sin\varphi cos\delta +cos\varphi sin\delta -cos\varphi ),\end{array}$$(A.1)

subject to the horizon condition:

$$\begin{array}{c}\hfill \delta \ge \varphi -\frac{\pi }{2}.\end{array}$$(A.2)

The maximum absolute delay is now analytically obtained by setting $\frac{\partial {\eta }_0}{\partial \delta }=0$. This gives us:

$$\begin{array}{c}\hfill \delta {|}_{\max }=\varphi -\frac{\pi }{2}{|{\eta }_0|}^{\max }=\frac{|\mathit{b}|}{c}(1+cos\varphi ).\end{array}$$(A.3)

This extremum also satisfies the horizon condition of Eq. (A.2) and the source lies exactly on the physical horizon, as was seen for the case of an arbitrary phase center (Sect. 3). To ensure that this reduction of the problem to the meridian plane can indeed be made without loss of generality, Eq. (5) was numerically maximized subject to Eq. (6) over H, δ and θ for δ₀ = 90°, |b| = 100 and ϕ = 47.38°. The maximum occurs at (H, δ, θ)≈(0, −42.62° ,270° ) with a value of 5.59 × 10⁻⁷ s while the minimum occurs at (H, δ, θ)≈(0, −42.62° ,90° ) with a value of −5.59 × 10⁻⁷ s. Both situations occur at a declination of ϕ − π/2 and have a delay value of (|b|/c)(1 + cos ϕ). The maximum occurs when the baseline points along the south and the minimum occurs when the baseline points along the north, as is expected for a source located toward the south. These results match with the predictions of Eq. (A.3) (irrespective of whether Eq. (A.2) is imposed). This shows that for the specific case of NCP phasing, the problem can be reduced to the meridian plane without loss of generality.

For drift scan observations, the horizon delay is given by |b|/c which occurs due to sources 90° away from the zenith. However, once the data is phased to the NCP (or any other location), we do not need to go 90° away from the phase center to get a delay equal to |b|/c. Setting |η₀|=|b|/c in Eq. (A.1), we get δ = sin⁻¹(cos ϕ − 1)+ϕ = 28.54° for ϕ = 47.38° (NenuFAR). So a source located more than 61.46° away from the phase center already will have a delay greater than |b|/c. For the horizon line in the power spectrum, we need to set |η₀|/sin ϕ = |b|/c (following Sect. 4) where sin ϕ is the projection of b in the plane perpendicular to $\widehat{\mathit{p}}$. The angle ζ away from the phase center beyond which the horizon line is steeper than the conventional horizon line is then given by:

$$\begin{array}{c}\hfill \zeta =\pi /2-\delta =\pi /2-{sin}^{-1}(cos\varphi -sin\varphi )-\varphi .\end{array}$$(A.4)

Inserting ϕ = 47.38° for NenuFAR, we get ζ = 45.99°. Therefore, the primary beam must be able to attenuate all power beyond an angle less than 45.99°, otherwise the foreground wedge will be steeper than that predicted by the conventional horizon line equation (Eq. (17)). Since Eq. (A.4) is a monotonically increasing function of ϕ, this effect is more pronounced for lower latitudes, where sources slightly away from the phase center can result in a foreground wedge steeper than the conventional horizon extent.

Appendix B: Sampling and gridding on a circular uv track

In Sect. 5.1, we find that sampling visibilities along uv tracks and gridding them using a convolution kernel results in higher power along the source direction, particularly for sources far away from the phase center. Here we investigate this effect analytically, for a circular track in the uv plane with a radius u₀, on which the interferometer samples the visibility. For a unit unpolarized source s′=(l, m), under the flat-sky approximation, the visibility function is given by:

$$\begin{array}{c}\hfill V(u,v)=e^{-2\pi i(\mathit{u}·{\mathit{s}}^{\prime })}=e^{-2\pi i(ul+vm)}\text{where}\mathit{u}=(u,v).\end{array}$$

The sampled visibility on the circular track is given by:

$$\begin{array}{cc}\hfill V_{\mathrm{track}}(u,v)& =e^{-2\pi i(ul+vm)}S(u,v),\hfill \\ \hfill \text{where}& S(u,v)=\{\begin{array}{cc}1\hfill & \text{if}\sqrt{u^2+v^2}=u_0\hfill \\ 0\hfill & \text{otherwise}.\hfill \end{array}\hfill \end{array}$$

To reduce the problem to one dimension, we can write down the expression of V_track in terms of β (the angle between u and s′).

$$\begin{array}{cc}\hfill V_{\mathrm{track}}(\beta )& =e^{-2\pi i|\mathit{u}||{\mathit{s}}^{\prime }|cos\beta }=e^{-2\pi i{u}_{0}sin\psi cos\beta }.\hfill \end{array}$$(B.1)

These sampled visibilities are finally distributed on the uv plane using a convolution kernel. We assume, for simplicity, that a two-dimensional Gaussian is used as the convolution kernel. If f is the FoV or the width of the Gaussian in the image space, the corresponding gridding kernel in the uv space is given by:

$$\begin{array}{c}\hfill G(u,v)=e^{-2\pi f^2(u^2+v^2)},\end{array}$$

and the convolved visibility is:

$$\begin{array}{cc}\hfill V_{\mathrm{conv}}(u,v)& ={\displaystyle {\int }_{u^{\prime }{v}^{\prime }}}{V}_{\mathrm{track}}(u^{\prime },v^{\prime })G(u-u^{\prime },v-v^{\prime })d{u}^{\prime }d{v}^{\prime }\hfill \\ \hfill & ={\displaystyle {\int }_{u^{\prime }{v}^{\prime }}}{V}_{\mathrm{track}}(u^{\prime },v^{\prime })e^{-2\pi f^2{|\mathit{u}-{\mathit{u}}^{\mathbf{\prime }}|}^2}d{u}^{\prime }d{v}^{\prime }.\hfill \end{array}$$

Now $V_{\mathrm{track}}(u\prime ,v\prime )$ is zero everywhere except where $\sqrt{u^2+v^2}=u_0$. On such a track, u and u′ form two sides of an isosceles triangle with the length of the equal sides being u₀ and the angle between the sides being β − β′. Thus, the length of the third side (|u − u′|) is $u_0\sqrt{2-2cos(\beta -\beta \prime )}$. Now since the convolution does not get any contribution from visibilities outside the circular track, the expression for the convolved visibility along the track becomes:

$$\begin{array}{cc}\hfill V_{\mathrm{conv}}(\beta )& ={\displaystyle {\int }_{{\beta }^{\prime }}}{V}_{\mathrm{track}}({\beta }^{\prime })e^{-4\pi f^2{u}_0^2(1-cos(\beta -{\beta }^{\prime }))}{u}_{0}d{\beta }^{\prime }\hfill \\ \hfill & =V_{\mathrm{track}}(\beta )\otimes G(\beta ),\text{where}G(\beta )=u_0{e}^{-4\pi f^2{u}_0^2(1-cos\beta )}.\hfill \end{array}$$

Here u₀dβ′ is an infinitesimal element along the track and ‘⊗’ denotes a convolution over β. Finally, inserting the expression of V_track(β) from Eq. (B.1), we get:

$$\begin{array}{cc}\hfill V_{\mathrm{conv}}(\beta )& =cos(2\pi u_{0}sin\psi cos\beta )\otimes G(\beta )\hfill \\ \hfill & +isin(2\pi u_{0}sin\psi cos\beta )\otimes G(\beta ).\hfill \end{array}$$(B.2)

Figure B.1 illustrates the effect of convolving V_track(β) with G(β) for u₀ = 10, ψ = 60° and f = 15°. Both the real and imaginary parts of V_track(β) are quasi-sinusoids that have fast oscillations near β = −π/2 and β = π/2 (top left panel). These fast oscillations are suppressed by convolution with the gridding kernel (top right panel) in both real and imaginary parts of V_conv(β) (bottom left panel). As a result, the amplitude of V_conv(β) is lowered near β = −π/2 and β = π/2 (bottom right panel). The bottom right panel also shows the effect of the FoV on the amplitude of the convolved visibilities. Larger FoVs (narrower gridding kernels) filter out less power and result in a wider signature. This is also seen in the middle and right panels of Fig. 5. Changing the radius of the uv track (u₀) does not change the width of the peak of |V_conv(β)| (not shown in the figure). This is because, even though the fluctuations of |V_track(β)| are faster in a longer uv track, G(β) is also narrower (since the same G(uv) is now evaluated on a longer uv track), thus exactly compensating for the effect. This means that the signature of a source in the uv plane is a sector with a higher power, and the FoV decides its opening angle Δβ. This is seen in the middle and right panels of Fig. 5. Also, the nature of the convolution kernel affects the exact shape of |V_conv(β)|. Thus, for non-Gaussian convolution kernels, the power need not be the highest at exactly β = 0. It should be noted that here we have assumed a single uv track. For multiple uv tracks, Eq. (B.2) is not exactly valid since the convolved visibility at a point on a given uv track will have a contribution from other uv tracks. Also, we have assumed the uv track is circular. Thus, the track runs perpendicular to the oscillations exactly at β = nπ (for integer n). At these points, both the real and imaginary parts of the visibility have the slowest fluctuations along the track. As a result, the convolution does not decrease the amplitude of the visibilities at β = nπ, leading to higher power in the source direction. For a non-circular track, the slowest oscillations along the track do not necessarily occur exactly at β = nπ. Thus, the signature of the source on the uv plane will be shifted accordingly, and cannot be predicted without knowing the exact uv coverage of the observation.

Fig. B.1. Impact of sampling visibilities due to a source 60° away from the phase center (ψ = 60°) on a circular uv track of radius 10 (u0 = 10) followed by convolution with a gridding kernel corresponding to a FoV of 15° (f = 15°). The different panels show the visibilities along the uv track (top left, Eq. (B.1)), the convolution kernel (top right), the sampled visibilities convolved with the kernel (bottom left, Eq. (B.2)) and the absolute value of the convolved visibilities (bottom right). The amplitudes of the convolved visibilities for a few other FoVs are overplotted in the bottom right panel as blue dashed lines.

Appendix C: Source delay validation

Given a specific source, the delay corresponding to the maximum power in the power spectrum is given by Eq. (18). The baseline for which we get this delay lies along the projection of $\widehat{\mathit{s}}-\widehat{\mathit{p}}$ on the Earth’s tangent plane. Incidentally, this is also the baseline orientation which results in the maximum delay, since it makes the smallest angle with the $\widehat{\mathit{s}}-\widehat{\mathit{p}}$ vector. As a result, this expression of delay not only corresponds to the situation of maximum power but also the maximum delay. We verify this by maximizing Eq. (5) over θ for Cyg A and Vir A, corresponding to t = 18 h.

Fig. C.1 shows how the numerically maximized delays due to Cyg A and Vir A compare against those predicted by Eq. (18), for all possible phase center positions. We see that the delays match for both sources and for all three latitudes.

Fig. C.1. Comparison between the delay due to a source obtained by numerical maximization of Eq. (5), against that expected from Eq. (18). The delays due to Cyg A (left two columns) and Vir A (right two columns) for three different latitudes (different rows) are shown. The columns titled “Numerical” correspond to delay values obtained by numerically maximizing Eq. (5) over θ. The columns titled “Equation” correspond to the delay values for the baseline responsible for maximum power as predicted by Eq. (18). The slanted stripes indicate regions where the phase center is below the physical horizon.

Appendix D: Source lines for specific cases

In Sect. 5.1 we find that sampling on uv tracks and gridding with a convolution kernel results in the visibilities due to a source having a specific signature in the uv plane. In the derivation of the source line equation (Eq. (22)), we assumed that the maximum power in the uv plane occurs along the direction of the source, which is exactly correct only for circular uv tracks. For telescopes located at latitudes close to the equator, the uv tracks can look very different, particularly when the phase center moves far away from the zenith. Here, we investigate the impact of non-circular uv tracks on the source lines, in the uv plane and the cylindrical power spectrum. We performed simulations with the fiducial array described in Sect. 6 placed at the location of SKA, with a zenith transit phase center (RA = 1 h, Dec = −26.82°). Three 1 h simulations centered at LSTs of 21 h, 1 h, and 5 h were performed, so that the phase center crosses the zenith at the midpoint of the simulation centered at LST = 1 h. The uv coverages for the three simulations are shown in the top row of Fig. D.1. We can see that the first and third simulations result in radially outward uv tracks along the bottom right and bottom left, respectively. The simulations were repeated for three different sources. The first source chosen is Fornax A (α = 3.4 h, δ = −37.21°), which lies in the bottom left corner of the image plane. We then repeated the simulations for a source placed at the same declination, but at α = 22.6 h, so that it lies in the bottom right corner of the image plane. Finally, we repeated the simulations for a source located toward the bottom of the image plane (α = 1 h, δ = −60°). The visibility amplitudes in the uv plane for all simulations are shown in Fig. D.1. We see that the source located at the bottom left corner results in higher power approximately along this direction for the first two LSTs, but for the third LST, the line of higher power is in a completely different direction. This happens because for a source on the bottom left corner, for the third LST, most of the uv tracks run along the direction of the source. As a result, there are fast fluctuations of the visibility due to the source along these uv tracks, which are suppressed by the convolution kernel used in gridding (as described in Sect. 5.1). The opposite scenario occurs for the simulations with the source on the bottom right since now the first LST simulation has uv tracks running in the direction of the source. For the simulations with the source at the bottom, the source line in the uv plane behaves as expected for all three LSTs, since the uv tracks do not run along the source direction for any of the three LSTs. To see the impact of the uv tracks running in the direction of the source on the cylindrical power spectrum, we filtered out the uv cells which lie more than 40° away from the source direction and constructed the power spectrum. Fig. D.1 shows that when the source is at the bottom left or the bottom right of the image plane, the power due to the source is reduced whenever we have uv tracks running along the source direction. When the source is at the bottom, this reduction of power does not occur, as we would expect when the uv tracks do not run along the source direction. The source line equation derived in this analysis (Eq. (22): black solid line) and the conventional source line equation with a flat-sky approximation (Eq. (23): black dotted line) are also shown in the plots. We see that the conventional source line equation can accurately predict the signature of the source only when the phase center is at the zenith (as described in Sect. 4.3), but the new equation can describe the signature in other situations as well.

Fig. D.1. Impact of uv coverage on source lines. The top row shows the baseline coverages for the three LST ranges for which the simulations are repeated. The other rows show the visibility amplitude in the uv plane and the cylindrical power spectra, for a source located at the bottom left, bottom right, and bottom of the image.

Appendix E: Additional figures

Fig. E.1. Cylindrical power spectra for the remaining full-sky simulations that were not shown in Fig. 7. For each location, the two rows correspond to the celestial pole (top) and zenith transit (bottom) phase center declinations. In the lowermost row, for the simulations centered at LST = 9h, 13h, and 17h, the phase center is below the horizon for part of the observation, and hence the visibilities are not predicted. The full synthesis line is not shown in this case since it is vertical (from Eq. (15)).

Fig. E.2. Cylindrical power spectra for the remaining single source simulations that were not shown in Fig. 8. Here the power spectra correspond to Cas A simulations at the locations of NenuFAR and LOFAR, and Cen A simulations at the location of SKA.

All Tables

All Figures

	Fig. 1. Effect of geometric phasing on the delay power spectrum. The two panels show the delay power spectra before (left) and after (right) phasing to the NCP for a full-sky simulation with NenuFAR. The dotted lines correspond to the standard horizon limit of Eq. (1), and the solid lines correspond to the extended horizon due to phasing from Eq. (4).
In the text

	Fig. 2. Coordinate systems used in the calculations. The origin of the coordinate system is placed at the location of the array for illustration. Since the sources are at an infinite distance, the exact location of the coordinate center is irrelevant.
In the text

Fig. 3. Maximum delay in visibilities measured by a 100 m baseline, for a given phase center declination, assuming a full synthesis. Left: maximum delay as a function of declination of the phase center δ0, obtained by numerical maximization of Eq. (5) over H, δ, H0, θ (black dots), the maximum delays predicted by the partial equation (Eq. (7): blue dashed line) and those predicted by the full equation (Eq. (8): red solid line). Right: values of H, δ, and H0 at the maximum. The vertical dotted lines separate regimes where |δ0 + ϕ|≤π/2 and |δ0 + ϕ|> π/2.

In the text

Fig. 4. Maximum delay in visibilities measured by a 100 m baseline for a given phase center location in the sky (H0, δ0). Left: Maximum delays, obtained by numerical maximization of Eq. (5) over H, δ, and θ. Right: Maximum delays expected from Eq. (9). The vertical dotted lines separate regimes where |δ0 + ϕ|≤π/2 and |δ0 + ϕ|> π/2. The slanted stripes indicate the regions where the phase center is below the physical horizon.

In the text

	Fig. 5. Signature of a source on the uv plane. The left panel is a schematic diagram illustrating the quasi-sinusoidal behavior of the phase of visibilities due to a source along a circular uv track. The middle and right panels show the visibility amplitude in the uv plane for simulations of Cyg A with the phase center at the NCP, with FoVs of 8° and 32°, respectively.
In the text

	Fig. 6. Inputs to the simulations used for validating the derived equations. Top: sky model used in the simulations performed for validating the horizon line equations. Bottom left: Local array configuration. Bottom right: example baseline coverage for a 1 h simulation with NCP phasing with the array placed at the location of NenuFAR.
In the text

Fig. 7. Cylindrical power spectra for full-sky simulations with the array placed at the locations of NenuFAR, LOFAR, and SKA corresponding to a phase center declination of 70° (−70° for SKA). The different columns correspond to different LSTs at which the 1 h simulations are centered. The dotted line is the conventional horizon line. The black solid line is the LST-dependent horizon line (Eq. (14)). The red dashed line is the horizon line for full synthesis (Eq. (15)). a0 is the altitude of the phase center at the LST at which the simulation is centered. The results of the remaining simulations for the celestial pole and zenith transit phase centers are shown in Fig. E.1.

In the text

Fig. 8. Cylindrical power spectra for Cyg A simulations at the locations of NenuFAR and LOFAR, and Vir A simulations at the location of SKA. For each location, the two rows correspond to two phase center declinations (δ0) of the celestial pole (upper row) and zenith transit (lower row). The different columns correspond to the different LST ranges for which the simulations were repeated. The dotted line is the conventional source line. The black solid lines correspond to the maximum and minimum source lines for the observation duration, obtained from the source line equation derived in this analysis (Eq. (22)). The remaining Cas A and Cen A simulation results are shown in Fig. E.2.

In the text

Fig. B.1. Impact of sampling visibilities due to a source 60° away from the phase center (ψ = 60°) on a circular uv track of radius 10 (u0 = 10) followed by convolution with a gridding kernel corresponding to a FoV of 15° (f = 15°). The different panels show the visibilities along the uv track (top left, Eq. (B.1)), the convolution kernel (top right), the sampled visibilities convolved with the kernel (bottom left, Eq. (B.2)) and the absolute value of the convolved visibilities (bottom right). The amplitudes of the convolved visibilities for a few other FoVs are overplotted in the bottom right panel as blue dashed lines.

In the text

Fig. C.1. Comparison between the delay due to a source obtained by numerical maximization of Eq. (5), against that expected from Eq. (18). The delays due to Cyg A (left two columns) and Vir A (right two columns) for three different latitudes (different rows) are shown. The columns titled “Numerical” correspond to delay values obtained by numerically maximizing Eq. (5) over θ. The columns titled “Equation” correspond to the delay values for the baseline responsible for maximum power as predicted by Eq. (18). The slanted stripes indicate regions where the phase center is below the physical horizon.

In the text

	Fig. D.1. Impact of uv coverage on source lines. The top row shows the baseline coverages for the three LST ranges for which the simulations are repeated. The other rows show the visibility amplitude in the uv plane and the cylindrical power spectra, for a source located at the bottom left, bottom right, and bottom of the image.
In the text

Fig. E.1. Cylindrical power spectra for the remaining full-sky simulations that were not shown in Fig. 7. For each location, the two rows correspond to the celestial pole (top) and zenith transit (bottom) phase center declinations. In the lowermost row, for the simulations centered at LST = 9h, 13h, and 17h, the phase center is below the horizon for part of the observation, and hence the visibilities are not predicted. The full synthesis line is not shown in this case since it is vertical (from Eq. (15)).

In the text

	Fig. E.2. Cylindrical power spectra for the remaining single source simulations that were not shown in Fig. 8. Here the power spectra correspond to Cas A simulations at the locations of NenuFAR and LOFAR, and Cen A simulations at the location of SKA.
In the text